A new preconditioned Richardson iterative method

Linear systems arise in various scientific and engineering applications, including differential equations, signal processing, and optimization [4, 5, 6, 7, 8, 12]. Therefore, it is crucial to find the efficient solution of a linear system in both theoretical and practical contexts. Iterative methods of solving linear systems have been extensively studied because they are more efficient than direct methods when solving large-scale problems. Furthermore, iterative methods have the advantage of being easily adaptable to different types of linear operators, which makes them versatile and applicable to a broad range of problems. Such operators frequently arise in practical applications, including image processing and machine learning.

In the present work, we propose a new iterative method for solving linear systems with bounded, self-adjoint operators and then present a new, modified version of the process by using the Chebyshev polynomials. The method utilizes a flexible preconditioner that adapts to the structure of the operator and exploits its self-adjointness. The resulting algorithms have favorable convergence properties and require fewer iterations than some existing methods. Also, we present the results of some numerical experiments that demonstrate the efficiency and effectiveness of our methods compared to the existing, state-of-the-art methods.

Consider the linear system

in a Hilbert space $\mathscr{H}$, where $A$ is a bounded and self-adjoint operator. The basic idea of an iterative method is to use an initial guess $x_0$, and then apply a sequence of updates to the solution, refining the solution with each iteration. There are many different iterative methods, each having its advantages and disadvantages, depending on the properties of the equation to be solved.

The most straightforward approach to an iterative solution is to rewrite (1) as a linear, fixed point iteration. If we consider the function $ℒ(x)=Ax+x-y$ on $\mathscr{H}$, then $x^{\ast }$ is a solution of $Ax=y$ if and only if $x^{\ast }$ is a fixed point of $ℒ$. Thus, it seems natural to consider fixed-point theorems.

One of the best and most widely used methods to do this is the Richardson iterative method [12]. The abstract Richardson iterative method has the form

where and $r_k$ is the residual vector $y-A{x}_k$. Note that we can rewrite (2) as

The Conjugate Gradient method is a more sophisticated method that uses information about the residual error to guide the iterative updates [12]. The update for the $k$th iteration is given by

where ${\alpha }_k$ is the step size, $p_k$ is the search direction given by

and

is the residual, with ${\beta }_k$ determined by the Conjugate Gradient method.

However, it is possible to precondition the linear equation (1) by multiplying both sides by an operator $B$ to obtain $BAx=By$, so that the convergence of the iterative methods is improved. In this case, the residual $By-BA{x}_k$ better reflects the actual error. This method is a very effective technique for solving differential equations, integral equations, and related problems [2, 9, 10].

In general, iterative methods are powerful tools for solving large systems of linear equations. By choosing an appropriate iterative method based on the properties of the operator under consideration, it is possible to obtain fast and accurate solutions to a wide range of linear systems.

Based on the properties of the operator $A$ in equation (1), positive constants $c_1$ and $c_2$ exist such that for each $u\in \mathscr{H}$,

Although the positive definiteness of $A$ is essential for methods such as conjugate gradient, our method can disregard this property. This is more evident in the Numerical results section.

It is not difficult to show that the optimal parameters $c_1$ and $c_2$ in equation (3) are equal to ${\lambda }_{\min }(A^2)$ and ${\lambda }_{\max }(A^2)$, respectively. Considering the properties of matrix $A$, there exist constants $c_1$ and $c_2$ that satisfy the aforementioned equation. It is not imperative to utilize the optimized coefficients within this context. Consequently, precise knowledge of the eigenvalues of the matrix $A^2$ is not mandatory; an approximation thereof is sufficient.
Now, based on the Richardson iterative method, we consider the iteration

In this case, the following lemma holds.

Note that the optimal convergence rate in the Richardson iterative method, obtained by letting $\alpha =\frac{2}{{\lambda }_{\min }(A)+{\lambda }_{\max }(A)}$, is $\frac{{\lambda }_{\max }(A)+{\lambda }_{\min }(A)}{{\lambda }_{\max }(A)+{\lambda }_{\min }(A)}$.
In the following theorem, we show that the convergence rate of the iterative method (4) is ${\left(\frac{c_2-c_1}{c_2+c_1}\right)}^2$.

Note: Assume $c_1^{\ast }={\lambda }_{\min }(A^2)$ and $c_2^{\ast }={\lambda }_{\max }(A^2)$. Let $c_1=c_1^{\ast }-\delta$ and $c_2=c_2^{\ast }+\delta$, then $c_1+c_2=c_1^{\ast }+c_2^{\ast }$ and $c_2-c_1=c_2^{\ast }-c_1^{\ast }+2\delta$. Therefore, since ${\left(\frac{c_2^{\ast }-c_1^{\ast }}{c_2^{\ast }+c_1^{\ast }}\right)}^2\le {\left(\frac{c_2-c_1}{c_2+c_1}\right)}^2$, Lemma 1 and so Theorem 2 hold in this case.

To summarize the results obtained so far, we present an algorithm that generates an approximate solution to equation (1) with prescribed accuracy.

In this section, the properties of the Chebyshev polynomials are used to modify the previous algorithm to accelerate the convergence rate.

Considering the iteration (3), let $u_n={\sum }_{k=1}^n{a}_{n_k}{x}_k$, where ${\sum }_{kj=1}^n{a}_{n_k}=1$. In this case, based on the proof of the previous theorem,

By setting $B={\left(I-\frac{2}{c_1+c_2}{A}^2\right)}^2$ and $P_n(x)={\sum }_{k=1}^n{a}_{n_k}{x}^k$, we see that

Since $A$ is invertible and self-adjoint, $A^2$ is a positive definite operator. Also, by the previous lemma, the spectrum of $B$ is a subset of the interval $[-\rho ,\rho ]$, where $\rho ={\left(\frac{c_2-c_1}{c_2+c_1}\right)}^2$. Therefore, in view of (6), the spectral theorem leads to

In the sequel, we aim to minimize ${\Vert x^{\ast }-u_n\Vert }_{\mathscr{H}}$. We have to find

where ${\mathbb{P}}_n:=\{P_n(x):P_n(x)$ is a polynomial of degree $n$ with $P_n(1)=1.\}$.

We investigate the solution to this problem in terms of the Chebyshev polynomials [3]. These polynomials are defined by

and satisfy the recurrence relations

The following lemma presents a minimization property of these polynomials which will be used later.

This lemma shows that the minimization problem (8) can be solved by setting $a=-\rho$ and $b=\rho$. These lead to

which solves (8).

Now, we rewrite $u_k$ by using the Chebyshev polynomials. First of all, combining (9) with the definition of the Chebyshev polynomials (the recurrence relation) for $n\ge 2$, we obtain

Replacing $x$ with $B$ and applying the resulting operator identity to $x^{\ast }-x_0$ yield

By (6) and the fact that $Q_n(x)$ is the solution of the minimization problem (8), the above equation can be recovered as

or equivalently,

Repeated application of the recurrence relations of the Chebyshev polynomials for $n\ge 2$ leads to

and hence

If we set

then according to the properties of the Chebyshev polynomials we obtain

Therefore, we can rewrite (10) as

or equivalently,

This yields

Also, a straightforward computation gives us the following recursive relation for ${\rho }_n$,

Now, based on the recursive relation (12), we design the following algorithm to approximately solve equation (1).

The following theorem investigates the convergence rate of the Algorithm.

It is expected that methods utilizing more spectral information to yield better results compared to those relying solely on matrix-vector products. However, the exact eigenvalues are not necessary for our method, as we assume that the coefficients $c_1$ and $c_2$ exist.

In this section, we present several examples to evaluate the efficiency and performance of our algorithms and compare them with several well-known algorithms in some cases. In addition, we compare our algorithms, namely, Algorithm 1 and Algorithm 2, with each other as well as with the Richardson and Conjugate Gradient (CG) algorithms in some cases.

The reported experiments were performed on a 64-bit 2.4 GHz system using MATLAB version 2010.

In the following two examples, we show the efficiency of our novel algorithms by using $3\times 3$ and $5\times 5$ systems, respectively. Both examples use the tolerance threshold $ϵ=0.001$.

As discussed in Section 2, it is unnecessary to utilize the optimal values of $c_1$ and $c_2$ as the eigenvalues of the matrix $A^2$; a mere approximation is sufficient. Methods such as the Power Method [13] and the Jacobi Method [13] can be employed to approximate the eigenvalues of $A^2$.

To demonstrate the efficacy of our algorithms with approximate values of $c_1$ and $c_2$, we employed several approximated values for $c_1\le {\lambda }_{\min }(A^2)$ and $c_2\ge {\lambda }_{\max }(A^2)$ in Example 5. The resulting data is summarized in Table 1. In this table, the first column corresponds to the optimal values of $c_1$ and $c_2$, while the second and third columns represent the data corresponding to the approximate values.

As shown in Table 1, the further we deviate from the optimal values ${\lambda }_{\min }(A^2)$ and ${\lambda }_{\max }(A^2)$, the number of iterations increases in both our algorithms, but the computational error decreases due to the increased number of iterations. This demonstrates that our algorithms also perform well with approximate values for $c_1$ and $c_2$.

Here, the error function can be defined by $Err(k)=\Vert x_k-x^{\ast }\Vert$ at each step $k$ of the iterative method, where $x_k$ is the $k$th approximation of the solution and $x^{\ast }$ is the exact solution of the system. We indicate the value of this function at the final iteration by $Err$. For the approximate solution obtained from Algorithm 1 in the previous example, this value is $6.4613\times {10}^{-4}$, which is a little less than that of Algorithm 2 with $Err=6.5223\times {10}^{-4}$. Nevertheless, these error function values indicate that the final approximation obtained from each of our two algorithms is accurate up to four decimal places. Fig. 1 shows how Algorithm 1 and Algorithm 2 converge to the exact solution of the system in Example 6.

Although the Richardson iterative method is a light calculation method that is used in many applications of iterative methods in the solution of linear systems, our examples show that in many cases, the number of iterations required to reach the desired solution in Richardson’s method is greater than that of our second method. For instance, in the following example and the special case $ϵ=0.001$, the number of iterations required in Richardson’s method is equal to $30322$ and the final error value is equal to $Err=2.1836\times {10}^{-4}$. But, in our second algorithm, this number is equal to $23896$ with an error value equal to $Err=1.2679\times {10}^{-4}$.

The optimal values for $c_1$ and $c_2$ in the previous example are ${\lambda }_{\min }(D^2)=0.3090$ and ${\lambda }_{\max }(D^2)=13556690$. As can be seen, in this example, there is no need to use optimal values for $c_1$ and $c_2$. Table 2 presents the numerical results of the iterations required to converge to the exact solution of the systems in Example 7.

As shown in Table 2, in this special example, the well-known method CG does not work at all, while our first algorithm works, although slowly!

Also, we can see in Table 2 that the required number of iterations with different values of $ϵ$ in Algorithm 2 is less than that of Richardson’s method.

According to the above example, if we ignore the condition of positive definiteness, then the algorithm of the Conjugate Gradient method fails to work. Moreover, if the matrix of the system has negative eigenvalues, then the Richardson iterative method fails to work either. The following examples illustrate such systems.