A Modification of Minimal Residual Iterative Method to Solve Linear Systems

Abstract

We give a modification of minimal residual iteration (MR), which is 1V-DSMR to solve the linear system 𝐴𝑥=𝑏. By analyzing, we find the modifiable iteration to be a projection technique; moreover, the modification of which gives a better (at least the same) reduction of the residual error than MR. In the end, a numerical example is given to demonstrate the reduction of the residual error between the 1V-DSMR and MR.

1. Introduction

One of the important computational problems in the applied science and engineering is the solution of 𝑛×𝑛 nonsingular linear systems of equations:𝐴𝑥=𝑏, (1.1)where 𝐴∈𝑅^𝑛×𝑛 is a symmetric positive definite matrix (referred to as an SPD matrix), 𝑏∈𝑅^𝑛 is given, and 𝑥∈𝑅^𝑛 is unknown. To solve this problem, usually an iterative method is spurred by demands, which can be found in excellent papers [1, 2]. Most of the existing practical iterative techniques for solving larger linear systems of (1.1) utilize a projection process in one way or another; see, for example, [3–9].

Projection techniques are presented in different forms in many other areas of scientific computing, which can be formulated in abstract Hilbert functional spaces and finite element spaces. Furthermore, projection techniques are the process in which one attempts to solve a set of equations by solving each separate equation by correcting so that it is small in some norm. The idea of projection process is to extract an approximate solution to (1.1) from a subspace of 𝑅^𝑛. Denote 𝒦 and ℒ the search subspace and the constraints subspace, respectively. Let 𝑚 be their dimension and 𝑥₀∈𝑅^𝑛 be an initial guess to the solution of (1.1). A projection method onto the subspace 𝒦 and orthogonal to ℒ is a process which finds an approximation solution 𝑥∈𝑅^𝑛 to (1.1) by imposing the Petrov-Galerkin conditions that 𝑥 belongs to the affine space 𝑥₀+𝒦 and the new residual vector is orthogonal to ℒ, that is,find𝑥∈𝑥₀+𝒦,suchthat𝑏−𝐴𝑥⟂ℒ. (1.2)

From this point of view, the basic iterative methods for solving (1.1), such as Gauss-Seidel Iteration (GS), Steepest Descent Iteration (SD), Minimal Residual Iteration (MR), and Residual Norm Steepest Descent Iteration (RNSD), all can be viewed as a special case of the projection techniques.

In [2], Ujević obtained a new iterative method for solving (1.1), which is considered as a modification of Gauss-Seidel method. In [10], Jing and Huang pointed that this iterative method is also a projection process and named this method as “one-dimensional double successive projection method” (referred to as 1D-DSPM). In the same paper [10], the authors obtained another iterative method, which is named as “two-dimensional double successive projection method” (referred to as 2D-DSPM). The theory indicates that 2D-DSPM gives a better reduction of error than 1D-DSPM.

2. Notations and Preliminaries

In this paper, we will consider the following linear system of equations:𝐴𝑥=𝑏, (2.1)where 𝐴∈𝑅^𝑛×𝑛 is not a symmetric but a positive definite matrix of order 𝑛,𝑏∈𝑅^𝑛 is a given element, and 𝑥 unknown. For linear systems (2.1), we can use the classical minimal residual iteration to solve, which can be found in [11]. Here, we will give a modification of minimal residual iteration to solve linear system (2.1). We call the modification method as one vector double successive MR (abbreviated 1V-DSMR) and compare reduction of the residual error at step 𝑘+1 between the modification iteration and original MR. Hence, we find that the modification iteration gives a better reduction of the residual error than the original MR. ⟨𝑥,𝑦⟩=𝑦^𝑇𝑥 denotes a vector inner product between the vector 𝑥,𝑦∈𝑅^𝑛.

We define the inner products as𝑎=⟨𝐴𝑣₁,𝐴𝑣₁⟩,𝑐=⟨𝐴𝑣₁,𝐴𝑣₂⟩=⟨𝐴𝑣₂,𝐴𝑣₁⟩,𝑑=⟨𝐴𝑣₂,𝐴𝑣₂⟩,𝑝=⟨𝑏−𝐴𝑥_𝑘,𝐴𝑣₁⟩=⟨𝑟_𝑘,𝐴𝑣₁⟩,𝑞=⟨𝑏−𝐴𝑥_𝑘,𝐴𝑣₂⟩=⟨𝑟_𝑘,𝐴𝑣₂⟩. (2.2)(2.3)

In this subsection, we will recall minimal residual iterative method and give some properties of this iteration.

For the linear system (2.1), we can use the following algorithm which is called minimal residual iteration, viewed in [11].

Algorithm 2.1. (1) Choose an initial guess solution 𝑥₀∈𝑅^𝑛 to (2.1), 𝑘∶=0.
(2) Calculate𝑟₀=𝑏−𝐴𝑥₀,𝛼₀=⟨𝐴𝑟₀,𝑟₀⟩⟨𝐴𝑟₀,𝐴𝑟₀⟩. (2.4)
(3) If 𝑟_𝑘=0, then stop; else,
(4) calculate𝑥_𝑘+1=𝑥_𝑘+𝛼_𝑘𝑟_𝑘,𝑟_𝑘+1=𝑏−𝐴𝑥_𝑘+1,𝛼_𝑘+1=⟨𝐴𝑟_𝑘+1,𝑟_𝑘+1⟩⟨𝐴𝑟_𝑘+1,𝐴𝑟_𝑘+1⟩,𝑘∶=𝑘+1. (2.5)
(5) Go to step (3).

The minimal residual iteration can be interpreted with projection techniques. Here we represent the principles of this method in our uniform notation as follows:𝑥_𝑘+1=𝑥_𝑘+𝛼𝑣₁, (2.6)where 𝛼=𝑝/𝑎.

If we choose 𝒦=span{𝑣₁} and ℒ=span{𝐴𝑣₁}, then (1.2) turns to find𝑥_𝑘+1∈𝑥_𝑘+𝒦,suchthat𝑏−𝐴𝑥_𝑘+1⟂ℒ, (2.7)where 𝑥_𝑘+1=𝑥_𝑘+𝛼𝑣₁.

Equation (2.7) can be represented in terms of inner product as⟨𝑏−𝐴𝑥_𝑘+1,𝐴𝑣₁⟩=0, (2.8)which is⟨𝑏−𝐴𝑥_𝑘−𝐴𝛼𝑣₁,𝐴𝑣₁⟩=⟨𝑏−𝐴𝑥_𝑘,𝐴𝑣₁⟩−𝛼⟨𝐴𝑣₁,𝐴𝑣₁⟩=𝑝−𝛼𝑎=0, (2.9)giving rise to 𝛼=𝑝/𝑎, which is the same as in (2.6).

If we choose a special 𝑣₁=𝑟_𝑘, then (2.6) is the minimal residual iteration (MR); up to now, it is clear that MR is a special case of projection methods.

For the MR, we have the following property.

Lemma 2.2. Let {𝑥_𝑘} and {𝑟_𝑘} be generated by Algorithm 2.1, then we have‖𝑟_𝑘+1‖²=‖𝑟_𝑘‖²−𝑝²𝑎, (2.10)where 𝑎,𝑝 are defined in (2.2) and (2.3), respectively.

Proof. Using (2.6), we obtain𝑟_𝑘+1=𝑏−𝐴𝑥_𝑘+1=𝑏−𝐴𝑥_𝑘−𝛼𝐴𝑣₁=𝑟_𝑘−𝛼𝐴𝑣₁, (2.11)where 𝑣₁=𝑟_𝑘. Hence we have‖𝑟_𝑘+1‖²=⟨𝑟_𝑘+1,𝑟_𝑘+1⟩=⟨𝑏−𝐴𝑥_𝑘−𝛼𝐴𝑣₁,𝑏−𝐴𝑥_𝑘−𝛼𝐴𝑣₁⟩=‖𝑟_𝑘‖²−𝛼⟨𝑟_𝑘,𝐴𝑣₁⟩−𝛼⟨𝐴𝑣₁,𝑟_𝑘⟩+𝛼²⟨𝐴𝑣₁,𝐴𝑣₁⟩=‖𝑟_𝑘‖²−2𝛼𝑝+𝛼²𝑎=‖𝑟_𝑘‖²−𝑝²𝑎. (2.12)
From Lemma 2.2, easily, we can get the reduction of the residual error of MR as follows:‖𝑟_𝑘‖²−‖𝑟_𝑘+1‖²=𝑝²𝑎. (2.13)

3. An Interpretation of 1V-DSMR with Projection Technique

In this section, we will give the modification of the minimal residual iterative method, which is abbreviated to 1V-DSMR; we can present this method in our uniform notation as follows:𝑥_𝑘+1=𝑥_𝑘+𝛼𝑣₁+𝛽𝑣₂, (3.1)where 𝛼=𝑝/𝑎 and 𝛽=(𝑎𝑞−𝑐𝑝)/𝑎𝑑.

We will have a two-step investigation of 1V-DSMR.

The first step is to choose 𝒦₁=span{𝑣₁} and ℒ₁=𝐴𝒦₁=span{𝐴𝑣₁}, then it turns into the proceeding of MR, so we have 𝛼=𝑝/𝑎.

The next step is a similar way to choose 𝒦₂=span{𝑣₂} and ℒ₂=𝐴𝒦₂=span{𝐴𝑣₂}; denote 𝑥_𝑘+1=𝑥_𝑘+𝛼𝑣₁ and (1.2) turns to find𝑥_𝑘+1∈𝑥_𝑘+1+𝒦₂,suchthat𝑏−𝐴𝑥_𝑘+1⟂ℒ₂, (3.2)where 𝑥_𝑘+1=𝑥_𝑘+1+𝛽𝑣₂.

Equation (3.2) can be represented as in terms of inner product as⟨𝑏−𝐴𝑥_𝑘+1,𝐴𝑣₂⟩=0, (3.3)which is⟨𝑏−𝐴𝑥_𝑘+1−𝛽𝐴𝑣₂,𝐴𝑣₂⟩=⟨𝑏−𝐴𝑥_𝑘−𝛼𝐴𝑣₁−𝛽𝐴𝑣₂,𝐴𝑣₂⟩=⟨𝑏−𝐴𝑥_𝑘,𝐴𝑣₂⟩−𝛼⟨𝐴𝑣₁,𝐴𝑣₂⟩−𝛽⟨𝐴𝑣₂,𝐴𝑣₂⟩=𝑞−𝛼𝑐−𝛽𝑑=0. (3.4)This gives rise to 𝛽=(𝑎𝑞−𝑐𝑝)/𝑎𝑑, which is the same as in (3.1).

If we choose a special 𝑣₁=𝑟_𝑘, then (3.1) is a modification of the minimal residual iteration, which is named as 1V-DSMR; up to now, it is clear that 1V-DSMR is also a special case of projection methods.

As 1V-DSMR, we have the following relation of residual errors.

Theorem 3.1. Let {𝑥_𝑘+1} generated by (3.1) and 𝑟̂_𝑘+1=𝑏−𝐴𝑥_𝑘+1, then we have‖𝑟̂_𝑘+1‖²=‖𝑟_𝑘̂‖²−1𝑎²𝑑[(𝑎𝑑+𝑐²)𝑝²+𝑎²𝑞²−2𝑎𝑐𝑝𝑞], (3.5)where 𝑟_𝑘̂ is the same as in Lemma 2.2, and 𝑎,𝑐,𝑑,𝑝,𝑞 are defined in (2.2) and (2.3).

Proof. Using (3.1), we have𝑟̂_𝑘+1=𝑏−𝐴𝑥_𝑘+1=𝑏−𝐴𝑥_𝑘−𝛼𝐴𝑣₁−𝛽𝐴𝑣₂=𝑟_𝑘̂−𝛼𝐴𝑣₁−𝛽𝐴𝑣₂. (3.6)By deduction, we get‖𝑟̂_𝑘+1‖²=⟨𝑟_𝑘̂−𝛼𝐴𝑣₁−𝛽𝐴𝑣₂,𝑟_𝑘̂−𝛼𝐴𝑣₁−𝛽𝐴𝑣₂⟩=⟨𝑟_𝑘̂,𝑟_𝑘̂⟩−2𝛼⟨𝑟_𝑘̂,𝐴𝑣₁⟩−2𝛽⟨𝑟_𝑘̂,𝐴𝑣₂⟩+𝛼²⟨𝐴𝑣₁,𝐴𝑣₁⟩+2𝛼𝛽⟨𝐴𝑣₁,𝐴𝑣₂⟩+𝛽²⟨𝐴𝑣₂,𝐴𝑣₂⟩=‖𝑟_𝑘̂‖²−2𝛼𝑝−2𝛽𝑞+𝛼²𝑎+2𝛼𝛽𝑐+𝛽²𝑑. (3.7)If we substitute 𝛼=𝑝/𝑎 and 𝛽=(𝑎𝑞−𝑐𝑝)/𝑎𝑑 into (3.7), then we obtain‖𝑟̂_𝑘+1‖²=‖𝑟_𝑘̂‖²−1𝑎²𝑑[(𝑎𝑑+𝑐²)𝑝²+𝑎²𝑞²−2𝑎𝑐𝑝𝑞]. (3.8)From Theorem 3.1, we also get a reduction of residual error of 1V-DSMR as follows:‖𝑟_𝑘̂‖²−‖𝑟̂_𝑘+1‖²=1𝑎²𝑑[(𝑎𝑑+𝑐²)𝑝²+𝑎²𝑞²−2𝑎𝑐𝑝𝑞]. (3.9)

Next we will depict the comparison results with respect to residual error reduction between 1V-DSMR and MR.

Theorem 3.2. 1V-DSMR gives a better (at least the same) reduction of the residual error than MR.

Proof. From the equalities (2.13) and (3.9), we have(‖𝑟_𝑘̂‖²−‖𝑟̂_𝑘+1‖²)−(‖𝑟_𝑘‖²−‖𝑟_𝑘+1‖²1)=𝑎²𝑑[(𝑎𝑑+𝑐²)𝑝²+𝑎²𝑞²𝑝−2𝑎𝑐𝑝𝑞]−²𝑎=1𝑎²𝑑(𝑐𝑝−𝑎𝑞)²≥0, (3.10)which proves the assertion of Theorem 3.2.

Theorem 3.2 implies that the residual error of MR is bigger than that of 1V-DSMR at 𝑘+1 step if the residual vectors 𝑟_𝑘 and 𝑟_𝑘̂ at the 𝑘th iteration are equal to each other.

4. A Particular Method of 1V-DSMR

In this section, particular 𝑣₁ and 𝑣₂ will be chosen, and an algorithm to interpret 1V-DSMR is obtained.

Since 1V-DSMR is a modification of minimal residual iteration, take 𝑣₁=𝑟_𝑘̂. In general, 𝑣₂ may be chosen in different ways. Here, we choose a particular 𝑣₂=𝑥_𝑘−1, then from (3.1), each step of 1V-DSMR is as follows:𝑥_𝑘+1=𝑥_𝑘+𝛼𝑟_𝑘̂+𝛽𝑥_𝑘−1,for𝑘=0,1,…, (4.1)where 𝑟₋₁=0 and𝑝𝛼=𝑎=⟨𝐴𝑟_𝑘̂,𝑟_𝑘̂⟩⟨𝐴𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩,𝛽=𝑎𝑞−𝑐𝑝=𝑟̂𝑎𝑑⟨𝐴_𝑘,𝐴𝑟_𝑘̂⟩⟨𝑟_𝑘̂,𝐴𝑥_𝑘−1⟩−⟨𝐴𝑟_𝑘̂,𝐴𝑥_𝑘−1⟩⟨𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩⟨𝐴𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩⟨𝐴𝑥_𝑘−1,𝐴𝑥_𝑘−1⟩. (4.2)In this case, the first step of 1V-DSMR is the same as minimal residual iteration.

The above results provide the following algorithm of 1V-DSMR.

Algorithm 4.1 (A particular implementation of 1V-DSMR in a generalized way). (1) Choose an initial guess solution 𝑥₀ of (2.1), 𝑘∶=0.
(2) Calculate𝑟₀̂=𝑏−𝐴𝑥₀,𝛼₀=⟨𝐴𝑟₀̂,𝑟₀̂⟩⟨𝐴𝑟₀̂,𝐴𝑟₀̂⟩. (4.3)
(3) Calculate𝑥₁=𝑥₀+𝛼₀𝑟₀̂=𝑥₀+⟨𝐴𝑟₀̂,𝑟₀̂⟩⟨𝐴𝑟₀̂,𝐴𝑟₀̂⟩𝑟₀̂,𝑟₁̂=𝑏−𝐴𝑥₁=𝑟₀̂−⟨𝐴𝑟₀̂,𝑟₀̂⟩⟨𝐴𝑟₀̂,𝐴𝑟₀̂⟩𝐴𝑟₀̂. (4.4)
(4) If 𝑟_𝑘̂=0, then stop; else,
(5) calculate𝛼_𝑘=⟨𝐴𝑟_𝑘̂,𝑟_𝑘̂⟩⟨𝐴𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩,𝛽_𝑘=⟨𝐴𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩⟨𝑟_𝑘̂,𝐴𝑥_𝑘−1⟩−⟨𝐴𝑟_𝑘̂,𝐴𝑥_𝑘−1⟩⟨𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩⟨𝐴𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩⟨𝐴𝑥_𝑘−1,𝐴𝑥_𝑘−1⟩,𝑥_𝑘+1=𝑥_𝑘+𝛼_𝑘𝑟_𝑘̂+𝛽_𝑘𝑥_𝑘−1,𝑟̂_𝑘+1=𝑏−𝐴𝑥_𝑘+1,𝑘∶=𝑘+1. (4.5)
(6) Goto step (4).

About Algorithm 4.1, we have the following basic property.

Theorem 4.2. If 𝐴 is a positive matrix, then the sequence of the iterations of Algorithm 4.1 converges to the solution of the linear system 𝐴𝑥=𝑏.

Proof. From the equality (3.9), we get‖𝑟_𝑘̂‖²−‖𝑟̂_𝑘+1‖²≥𝑝²𝑎=⟨𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩²⟨𝐴𝑟_𝑘̂,𝐴𝑟_𝑘̂⟩≥𝜆²_min‖𝑟_𝑘̂‖⁴𝜆²_max‖𝑟_𝑘̂‖²=𝜆²_min𝜆²_max‖𝑟_𝑘̂‖². (4.6)
This means that the sequence ‖𝑟_𝑘̂‖² is a decreasing and bounded one. Thus, the sequence in question is convergent implying that the left-hand side tends to zero. Obviously, ‖𝑟_𝑘̂‖² tends to zero, and the proof is complete.

5. Numerical Examples

In this section, we use examples to further examine the effectiveness and show the advantages of 1V-DSMR over MR.

We compare the numerical behavior of Algorithm 4.1 with Algorithm 2.1. All the tests are performed by MATLAB 7.0. Because of the influence of the error of roundoff, we regard the matrix 𝐴 as zero matrix if ‖𝐴‖<10⁻¹⁰.

For convenience of comparison, consider the two-dimensional partial differential equations on the unit square region Ω=[0,1]×[0,1] of the form−(𝑝𝑢_𝑥)_𝑥−(𝑞𝑢_𝑦)_𝑦+𝑟𝑢_𝑥+(𝑟𝑢)_𝑥+𝑠𝑢_𝑦+(𝑠𝑢)_𝑦+𝑡𝑢=𝑓, (5.1)where 𝑝,𝑞,𝑟,𝑠, and 𝑡 are all given real valued function of 𝑥 and 𝑦, which are as follows: 𝑝=𝑒^−𝑥𝑦,𝑞=𝑒^𝑥𝑦,𝑟=20(𝑥+𝑦),𝑠=10(𝑥+𝑦), and 𝑡=1/(1+𝑥+𝑦).

Here, we use a five-point finite difference scheme to discretize the above problem with a uninform grid of mesh spacing Δ𝑥=Δ𝑦=1/𝑚 in 𝑥 and 𝑦 directions, respectively; we can obtain a matrix of order 𝑚×𝑚 as 𝑚 varies, which is called PDE matrix. Now, we take 𝑚=30, then we get a matrix, which is called PDE900, and denoted by 𝑃. It is easy to check that 𝑃 is real unsymmetrical and nonsingular.

If we take the coefficient matrix 𝐴 of linear system (2.1) as 𝑃, the right vector 𝑏=(1,1,…,1)^𝑇, and initial iterative vector 𝑥₀=𝑏, then we use Algorithms 2.1 and 4.1 to compute the linear system (2.1), respectively. The comparison results between MR and 1V-DSMR are shown in Figure 1.

Figure 1 

Comparison of convergence curve of residual norm (1.1) between Algorithm 2.1 and Algorithm 4.1.

From Figure 1, we can see that the convergence velocity of Algorithm 4.1 is always faster than that of Algorithm 2.1. In fact, when we use Algorithm 2.1 to compute the linear system (2.1), we only need to iterate 814 steps, and the residual norm is ‖𝑟₈₁₄‖≤9.8625×10⁻¹¹. While using Algorithm 4.1, we only need to iterate 647 steps, and the residual norm is ‖𝑟₆₄₇‖≤9.8465×10⁻¹¹.

Acknowledgment

This project was granted financial support from Shanghai Science and Technology Committee (no. 062112065), Shanghai Priority Academic Discipline Foundation, The University Young Teacher Sciences Foundation of Anhui Province (no. 2006jql220zd) and PhD Program Scholarship Fund of ECNU2007.