Augmented Arnoldi-Tikhonov Regularization Methods for Solving Large-Scale Linear Ill-Posed Systems

Abstract

We propose an augmented Arnoldi-Tikhonov regularization method for the solution of large-scale linear ill-posed systems. This method augments the Krylov subspace by a user-supplied low-dimensional subspace, which contains a rough approximation of the desired solution. The augmentation is implemented by a modified Arnoldi process. Some useful results are also presented. Numerical experiments illustrate that the augmented method outperforms the corresponding method without augmentation on some real-world examples.

1. Introduction

We consider the iterative solution of a large system of linear equations 𝐴𝑥=𝑏, (1) where 𝐴∈ℝ^𝑛×𝑛 is nonsymmetric and nonsingular, and 𝑏∈ℝ^𝑛. We further assume that the coefficient matrix 𝐴 is of ill-determined rank; that is, all its singular values decay gradually to zero, with no gap anywhere in the spectrum. Such systems are often referred to as linear discrete ill-posed problems and arise from the discretization of ill-posed problems such as Fredholm integral equations of the first kind with a smooth kernel. The right-hand side 𝑏 of (1) is assumed to be contaminated by an error 𝑒∈ℝ^𝑛, which may stem from discretization or measurement inaccuracies. Thus, 𝑏=𝑏̂+𝑒, where 𝑏̂ is the unknown error-free right-hand side vector.

We would like to compute the solution 𝑥 of the linear system of equations with the error-free right-hand side 𝑏̂, 𝐴𝑥=𝑏̂. (2) However, since the right-hand side in (2) is not available, we seek to determine an approximation of 𝑥 by solving the available system (1) or a modification. Due to the ill conditioning of 𝐴, the system (1) has to be regularized in order to compute a useful approximation of 𝑥. Perhaps the best known regularization method is Tikhonov regularization [1–3], which in its simplest form is based on the minimization problem min_{𝑥∈ℝ𝑛}‖𝐴𝑥−𝑏‖²+1𝜇‖𝑥‖², (3) where 𝜇>0 is a regularization parameter. Here and throughout this paper ‖⋅‖ denotes the Euclidean vector norm or the associated induced matrix norm.

After regulating the system (1), we need to compute the solution 𝑥_𝜇 of the minimization problem (3). Such a vector 𝑥_𝜇 is also the solution of 𝐴^𝑇1𝐴+𝜇𝐼𝑥=𝐴^𝑇𝑏. (4) Here and in the following, 𝐼 denotes the identity matrix, whose dimension is conformed with the dimension used in the context. If 𝜇 is far away from zero, then, due to the ill conditioning of 𝐴, 𝑥_𝜇 is badly computed while, if 𝜇 is close to zero, 𝑥_𝜇 is well computed, but the error 𝑥_𝜇−𝑥 is quite large. Thus, the choice of a good value for 𝜇 is fairly important. Several methods have been proposed to obtain an effective value for 𝜇. For example, if the norm of the error 𝑒 or a fairly accurate estimate is known, the regularization parameter is quite easy to determine by application of the discrepancy principle. The discrepancy principle proposes that the regularization parameter 𝜇 can be chosen so that the discrepancy 𝑏−𝐴𝑥_𝜇 satisfies 𝑏−𝐴𝑥_𝜇=𝜂𝜀, (5) where 𝜀=‖𝑒‖ and 𝜂>1 is a constant; see, for example, [4] for further details on the discrepancy principle.

The singular value decomposition [5] of 𝐴 can be used to determine the solution 𝑥_𝜇 of the minimization problem (3). For an overview of numerical methods for computing the SVD, we refer to [6]. We remark that the computational effort required to compute the SVD is quite high even for moderately sized matrices.

Many numerical methods using Krylov subspaces have been proposed for the solution of large-scale Tikhonov regularization problems (3). The main idea of such algorithms has been to first project the large problems onto some Krylov subspace to produce problems with small size and then solve the small-sized problems by the SVD. For instance, several well-established methods based on the Lanczos bidiagonalization process have been proposed for the solution of the minimization problem (3); see [7–10] and references therein. These methods use the Lanczos bidiagonalization process to construct a basis of the Krylov subspace: 𝒦_𝑚𝐴^𝑇𝐴,𝐴^𝑇𝑏=span𝐴^𝑇𝑏,𝐴^𝑇𝐴𝐴^𝑇𝑏,…,𝐴^𝑇𝐴^𝑚−1𝐴^𝑇𝑏. (6) We remark that each Lanczos bidiagonalization step needs two matrix-vector product evaluations, one with 𝐴 and the other with 𝐴^𝑇. Other methods using the Krylov subspace 𝒦_𝑚(𝐴,𝐴𝑏)=span𝐴𝑏,𝐴²𝑏,…,𝐴^𝑚𝑏 (7) as the projection subspace have been also designed. For example, Lewis and Reichel [11] proposed to exploit the Arnoldi process to produce a basis of the Krylov subspace 𝒦_𝑚(𝐴,𝐴𝑏) to obtain an approximation of the solution of the Tikhonov regularization problem (3). Since each Arnoldi decomposition step requires only one matrix-vector evaluation with 𝐴, the approach based on the Arnoldi process may require fewer matrix-vector product evaluations than that based on the Lanczos bidiagonalization process. Moreover, the methods based on the Arnoldi process do not require the adjoint matrix 𝐴^𝑇 and, hence, are more appropriate to problems for which the adjoint is difficult to evaluate. For such problems we refer to [12]. A similar Tikhonov regularization method based on generalized Krylov subspace is proposed in [13].

Some numerical methods without using the Tikhonov regularization technique have been already proposed to solve the large-scale linear discrete ill-posed problem (1). These methods include the range-restricted GMRES (RRGMRES) method [14, 15], the augmented range-restricted GMRES (ARRGMRES) method [16], and the flexible GMRES (FGMRES) method [17]. The RRGMRES method determines the 𝑚th approximation 𝑥_𝑚 of (1) by solving the minimization problem min_{𝑥∈𝒦𝑚(𝐴,𝐴𝑏)}‖𝐴𝑥−𝑏‖,𝑚=1,2,…,𝑥₀=0. (8) The regularization is implemented by choosing a suitable dimension number 𝑚; see, for example, [18]. The ARRGMRES method augments the Krylov subspace 𝒦_𝑚(𝐴,𝐴𝑏) by a low-dimensional user-supplied subspace. The low-dimensional subspace is determined by vectors that are able to represent the known features of the desired solution. The augmented method can yield approximate solutions of higher accuracy than the RRGMRES method if the Krylov subspace 𝒦_𝑚(𝐴,𝐴𝑏) does not allow representation of the known features.

In this paper, we propose a new iterative method, named augmented Arnoldi-Tikhonov regularization method, for solving large-scale linear ill-posed systems (1). The new method is deduced by combining the Tikhonov regularization technique and the augmentation technique. The augmentation is implemented by a modified Arnoldi process.

The following summarizes the structure of this paper. Section 2 describes the augmented Arnoldi-Tikhonov regularization method and some useful results. Some real-world examples are presented in Section 3. Section 4 contains the conclusions.

2. Augmented Arnoldi-Tikhonov Regularization Method

We attempt to improve the Arnoldi-Tikhonov regularization method proposed in [11] by augmenting the Krylov subspace 𝒦_𝑚(𝐴,𝐴𝑏) by a 𝑘-dimensional subspace 𝒲, which contains a rough approximation of the desired solution of (2). Then, the subspace of projection we will exploit in the following is of the form 𝒦_𝑚(𝐴,𝐴𝑏)+𝒲=span𝐴𝑏,𝐴²𝑏,…,𝐴^𝑚𝑏+𝒲. (9)

Let 𝑊 be the 𝑛 by 𝑘 matrix whose columns form an orthonormal basis of the space 𝒲. For the purpose of augmentation by 𝒲, we apply the modified Arnoldi process [16] to construct the modified Arnoldi decomposition 𝐴𝑊,𝑉_𝑚=𝑈,𝑉_𝑚+1𝐻, (10) where 𝑈∈ℝ^𝑛×𝑘, 𝑉_𝑚+1∈ℝ^{𝑛×(𝑚+1)}, [𝑈,𝑉_𝑚+1] has orthonormal columns, and 𝐻∈ℝ^{(𝑚+𝑘+1)×(𝑚+𝑘)} is an upper Hessenberg matrix. We point out that the leading principal 𝑘 by 𝑘 submatrix of 𝐻 is the upper triangular matrix in the QR factorization [5] of 𝐴𝑊; that is, 𝐻 is of the formℎ𝐻=¹¹ℎ₁₂⋯ℎ_1𝑘ℎ_1,𝑘+1ℎ_1,𝑘+2⋯ℎ_1,𝑘+𝑚ℎ⋱⋱⋮⋮⋮⋯⋮_{𝑘−1,𝑘−1}ℎ_{𝑘−1,𝑘}ℎ_{𝑘−1,𝑘+1}ℎ_{𝑘−1,𝑘+2}⋯ℎ_{𝑘−1,𝑘+𝑚}ℎ_𝑘𝑘ℎ_𝑘,𝑘+1ℎ_𝑘,𝑘+2⋯ℎ_{𝑘,𝑘+𝑚}ℎ_{𝑘+1,𝑘+1}ℎ_{𝑘+1,𝑘+2}⋯ℎ_{𝑘+1,𝑘+𝑚}ℎ_{𝑘+2,𝑘+1}ℎ_{𝑘+2,𝑘+2}⋯ℎ_{𝑘+2,𝑘+𝑚}ℎ⋱⋱⋮_{𝑘+𝑚,𝑘+𝑚−1}ℎ_{𝑘+𝑚,𝑘+𝑚}ℎ_{𝑘+𝑚+1,𝑘+𝑚}. (11)The modified Arnoldi relation (10) is shared by other methods such as GMRES-E [19] and FGMRES [20], which are also augmented type methods.

The modified Arnoldi process is outlined in Algorithm 1. We remark that the modified Arnoldi process in this paper is slightly different from the one used by Morgan [19] to augment the Krylov subspace with some approximate eigenvectors. In his method, the augmenting vectors are put after the Krylov vectors while in Algorithm 1, the augmenting subspace containing a rough approximation solution is included in the projection subspace from the beginning. In general, this can give better results at the start of the regularization method proposed in this paper.

Input: 𝐴∈ℝ𝑛×𝑛, 𝑏∈ℝ𝑛, 𝑊∈ℝ𝑛×𝑘, 𝑚. Output: 𝑈∈ℝ𝑛×𝑘, 𝑉𝑚+1∈ℝ𝑛×(𝑚+1), 𝐻∈ℝ(𝑚+𝑘+1)×(𝑚+𝑘).  (1) Compute QR-factorization 𝐴𝑊=𝑈𝐻.  (2) Compute 𝑣1=𝐴𝑏−𝑈(𝑈𝑇𝐴𝑏), 𝑣1=𝑣1/𝑣1, and set 𝑉1=[𝑣1].  (3) For 𝑗=1,2,…,𝑚 Do:     𝑣𝑗+1∶=𝐴𝑣𝑗;     For 𝑖=1,2,…,𝑘 Do:      ℎ𝑖,𝑗+𝑘∶=𝑈(∶,𝑖)𝑇𝑣𝑗+1;      𝑣𝑗+1∶=𝑣𝑗+1−ℎ𝑖,𝑗+𝑘𝑈(∶,𝑖);     End For     For 𝑖=1,2,…,𝑗 Do:      ℎ𝑖+𝑘,𝑗+𝑘∶=𝑣𝑇𝑖𝑣𝑗+1;      𝑣𝑗+1∶=𝑣𝑗+1−ℎ𝑖+𝑘,𝑗+𝑘𝑣𝑖;     End For     ℎ𝑗+𝑘+1,𝑗+𝑘∶=𝑣𝑗+1/𝑣𝑗+1;     𝑣𝑗+1∶=𝑣𝑗+1/ℎ𝑗+𝑘+1,𝑗+𝑘;     𝑉𝑗+1∶=[𝑉𝑗,𝑣𝑗+1];      End For

Algorithm 1 

Modified Arnoldi process.

We remark that a loss of orthogonality can occur when the algorithm progresses; see [21]. A remedy is the so-called reorthogonalization where the current vector has to be orthogonalized against previously created vectors. One can choose between a selective reorthogonalization or a full reorthogonalization against all vectors in the current augmented subspace. In this paper we only use the full reorthogonalization. The full reorthogonalization can be done as a classical or modifed Gram-Schmidt orthogonalization; see [21] for details.

We now seek to determine an approximate solution 𝑥_𝜇,𝑚 of (3) in the augmented Krylov subspace 𝒦_𝑚(𝐴,𝐴𝑏)+𝒲. After computing the QR factorization [𝑊,𝑉_𝑚]=𝑄𝑅 with 𝑄∈ℝ^{𝑛×(𝑚+𝑘)} having orthonormal columns and 𝑅∈ℝ^{(𝑚+𝑘)×(𝑚+𝑘)} being an upper triangular matrix, we substitute 𝑥=[𝑊,𝑉_𝑚]𝑦, 𝑦∈ℝ^𝑚+𝑘, into (3). It yields the reduced minimization problem min_{𝑦∈ℝ𝑚+𝑘}𝐴𝑊,𝑉_𝑚𝑦−𝑏²+1𝜇𝑊,𝑉_𝑚𝑦²=min_{𝑦∈ℝ𝑚+𝑘}𝑈,𝑉_𝑚+1𝐻𝑦−𝑏²+1𝜇𝑄𝑅𝑦²=min_{𝑦∈ℝ𝑚+𝑘}𝐻𝑦−[𝑈,𝑉_𝑚+1]^𝑇𝑏²+‖(𝐼−𝑃)𝑏‖²+1𝜇𝑅𝑦², (12) in which 𝑃=[𝑈,𝑉_𝑚+1][𝑈,𝑉_𝑚+1]^𝑇 is an orthogonal projector onto span([𝑈,𝑉_𝑚+1]). Obviously, the reduced minimization problem (12) is equivalent to min_{𝑦∈ℝ𝑚+𝑘}𝐻𝑦−[𝑈,𝑉_𝑚+1]^𝑇𝑏²+1𝜇𝑅𝑦²=min_{𝑦∈ℝ𝑚+𝑘}𝐻1√𝜇𝑅𝑦−[𝑈,𝑉_𝑚+1]^𝑇𝑏0². (13) The normal equations of the minimization problem (13) is 𝐻^𝑇1𝐻+𝜇𝑅^𝑇𝑅𝑦=𝐻^𝑇𝑈,𝑉_𝑚+1^𝑇𝑏. (14) We denote the solution of the minimization problem (13) by 𝑦_𝜇,𝑚. Then, from (14) it follows that 𝑦_𝜇,𝑚=𝐻^𝑇1𝐻+𝜇𝑅^𝑇𝑅⁻¹𝐻^𝑇𝑈,𝑉_𝑚+1^𝑇𝑏. (15) The approximate solution of (3) is 𝑥_𝜇,𝑚=𝑊,𝑉_𝑚𝑦_𝜇,𝑚. (16)

Since the matrix 𝐻^𝑇𝐻+(1/𝜇)𝑅^𝑇𝑅 has a larger condition number than the matrix [𝐻^𝑇,(1/√𝜇)𝑅^𝑇]^𝑇, we apply the QR factorization of [𝐻^𝑇,(1/√𝜇)𝑅^𝑇]^𝑇 to obtain the solution 𝑦_𝜇,𝑚 of (13) instead. The QR factorization of [𝐻^𝑇,(1/√𝜇)𝑅^𝑇]^𝑇 can be implemented by a sequence of Givens rotations.

Define 𝜙_𝑚𝜇=𝑏−𝐴𝑥_𝜇,𝑚². (17) Substituting 𝑥_𝜇,𝑚=[𝑊,𝑉_𝑚]𝑦_𝜇,𝑚 into (17) and using the modified Arnoldi decomposition (10) yield 𝜙_𝑚𝜇=𝑏−𝐴𝑊,𝑉_𝑚𝑦_𝜇,𝑚²=𝑏−𝑈,𝑉_𝑚+1𝐻𝑦_𝜇,𝑚²=[𝑈,𝑉_𝑚+1]^𝑇𝑏−𝐻𝑦_𝜇,𝑚²+‖(𝐼−𝑃)𝑏‖². (18) Substituting 𝑦_𝜇,𝑚 in (15) into (18), we obtain 𝜙_𝑚𝜇=𝐼−𝐻𝐻^𝑇1𝐻+𝜇𝑅^𝑇𝑅⁻¹𝐻^𝑇𝑈,𝑉_𝑚+1^𝑇𝑏²+‖(𝐼−𝑃)𝑏‖². (19) Note that 𝐼−𝐻𝐻^𝑇1𝐻+𝜇𝑅^𝑇𝑅⁻¹𝐻^𝑇=𝜇𝐻𝑅⁻¹𝐻𝑅^−1𝑇+𝐼⁻¹. (20) Therefore, the function 𝜙_𝑚(𝜇) can be expressed as 𝜙_𝑚𝜇=𝜇𝐻𝑅⁻¹𝐻𝑅^−1𝑇+𝐼⁻¹𝑈,𝑉_𝑚+1^𝑇𝑏²+‖(𝐼−𝑃)𝑏‖². (21)

Concerning the properties of 𝜙_𝑚(𝜇), we have the following results, which are similar to those of Theorem 2.1 in [11] for the Arnoldi-Tikhonov regularization method.

Theorem 1. The function 𝜙_𝑚(𝜇) has the representation 𝜙_𝑚𝜇=𝑏^𝑇𝑈,𝑉_𝑚+1𝜇𝐻𝑅⁻¹𝐻𝑅^−1𝑇+𝐼⁻²×𝑈,𝑉_𝑚+1^𝑇𝑏+‖(𝐼−𝑃)𝑏‖². (22) Assume that 𝐴𝑏≠0 and (𝐻𝑅⁻¹)^𝑇[𝑈,𝑉_𝑚+1]^𝑇𝑏≠0. Then 𝜙_𝑚 is strictly decreasing and convex for 𝜇≥0 with 𝜙_𝑚(0)=‖𝑏‖². Moreover, the equation 𝜙_𝑚𝜇=𝜏 (23) has a unique solution 𝜇_𝜏,𝑚, such that 0<𝜇_𝜏,𝑚<∞, for any 𝜏 with 𝑃_{𝒩((𝐻𝑅−1)}^𝑇₎𝑈,𝑉_𝑚+1^𝑇𝑏²+‖(𝐼−𝑃)𝑏‖²<𝜏<‖𝑏‖², (24) where 𝑃_{𝒩((𝐻𝑅−1)}^𝑇₎ denotes the orthogonal projector onto 𝒩((𝐻𝑅⁻¹)^𝑇).

Proof. The proof follows the same argument of the proof of Theorem 2.1 in [11] and therefore is omitted.

We easily obtain the following theorem, of which the proof is almost the same as that of Corollary 2.2 in [11].

Theorem 2. Assume that the modified Arnoldi process breaks down at step 𝑗. Then the sequence {𝑠_𝑖}^𝑗_𝑖=0 defined by 𝑠₀=‖𝑏‖²,𝑠_𝑖=𝑃_{𝒩((𝐻𝑅−1)}^𝑇₎𝑈,𝑉_𝑚+1^𝑇𝑏²+‖(𝐼−𝑃)𝑏‖²𝑠,0<𝑖<𝑗,_𝑗=0, (25) is decreasing.

We apply the discrepancy principle to the discrepancy 𝑏−𝐴𝑥_𝜇,𝑚 to determine an appropriate regularization parameter 𝜇 so that it satisfies (5). To make the equation 𝜙_𝑚𝜇=𝜂²𝜀² (26) have a solution, it follows from Theorem 1 that the input parameter 𝑚 of the modified Arnoldi should be chosen so that 𝑠_𝑚<𝜂²𝜀². To simplify the computations, we ignore the first term in 𝑠_𝑚. Then, the smallest iterative step number, denoted by 𝑚_min, of the modified Arnoldi process is chosen so that ‖(𝐼−𝑃)𝑏‖²<𝜂²𝜀². (27) We can improve the quality of the computed solution by choosing the practical iterative step number 𝑚 somewhat larger than 𝑚_min.

After choosing the number 𝑚 of the modified Arnoldi iterative steps, the regularization parameter 𝜇_𝑚 is determined by solving the nonlinear equation 𝜙_𝑚(𝜇)=𝜂²𝜀². Many numerical methods have been proposed for the solution of a nonlinear equation, including Newton's method [22], super-Newton's [23] method, and Halley's method [24]. For the specific nonlinear equation 𝜙_𝑚(𝜇)=𝜂²𝜀², Reichel and Shyshkov proposed a new zero-finder method in their new paper [25]. In this paper, we still make use of Newton's method to obtain the regularization parameter 𝜇_𝑚.

In Algorithm 2, we outline the augmented Arnoldi-Tikhonov regularization method, which is used to solve large-scale linear ill-posed systems (1).

Input: 𝐴∈ℝ𝑛×𝑛, 𝑏∈ℝ𝑛, 𝑊∈ℝ𝑛×𝑘, 𝜀, 𝜂, 𝑚0. Output: 𝑥𝜇,𝑚, 𝜇𝑚, 𝑚=𝑚min+𝑚0.  (1) Compute the modified Arnoldi decomposition 𝐴[𝑊,𝑉𝑚]=[𝑈,𝑉𝑚+1]𝐻, with 𝑚=𝑚min+𝑚0,    where 𝑚min is the smallest number such that ‖(𝐼−𝑃)𝑏‖2<𝜂2𝜀2.  (2) Compute the solution 𝜇𝑚 of the equation 𝜙𝑚(𝜇)=𝜂2𝜀2 by Newton's method.  (3) Compute the solution 𝑦𝜇,𝑚 of the least-squares problem (13) and obtain the approximate    solution 𝑥𝜇,𝑚=[𝑊,𝑉𝑚]𝑦𝜇,𝑚.

Algorithm 2 

Augmented Arnoldi-Tikhonov regularization method.

Newton's method requires to evaluate 𝜙_𝑚(𝜇) and its first derivative with respect to 𝜇 for computing approximations 𝜇_𝑚^(𝑗) of 𝜇_𝑚 for 𝑗=0,1,2,…. Let 𝑧_𝜇,𝑚=𝜇𝐻𝑅⁻¹𝐻𝑅^−1𝑇+𝐼⁻¹𝑈,𝑉_𝑚+1^𝑇𝑏; (28) that is, 𝑧_𝜇,𝑚 satisfies the system of linear equations 𝜇𝐻𝑅⁻¹𝐻𝑅^−1𝑇+𝐼𝑧=𝑈,𝑉_𝑚+1^𝑇𝑏. (29) Note that the above linear system is the normal equations of the least-squares problem min_{𝑧∈ℝ𝑚+𝑘+1}𝜇^1/2𝐻𝑅^−1𝑇𝐼0𝑧−U,𝑉_𝑚+1^𝑇𝑏. (30) For numerical stability, we compute the vector 𝑧_𝜇,𝑚 by solving the least-squares problem. Then, the 𝜙_𝑚(𝜇) is evaluated by computing 𝜙_𝑚𝜇=𝑧^𝑇_𝜇,𝑚𝑧_𝜇,𝑚+‖(𝐼−𝑃)𝑏‖². (31) It is easy to show that the first derivative of 𝜙^′_𝑚(𝜇) can be written as 𝜙^′_𝑚𝜇=−2𝑧^𝑇_𝜇,𝑚𝑤_𝜇,𝑚, (32) where 𝑤_𝜇,𝑚=𝜇𝐻𝑅⁻¹𝐻𝑅^−1𝑇+𝐼⁻¹𝐻𝑅⁻¹𝐻𝑅^−1𝑇𝑧_𝜇,𝑚. (33) Hence, we may compute 𝑤_𝜇,𝑚 by solving a least-squares problem analogous to the above with the vector [𝑈,𝑉_𝑚+1]^𝑇𝑏 replaced by 𝐻𝑅⁻¹(𝐻𝑅⁻¹)^𝑇𝑧_𝜇,𝑚.

The algorithm for implementing the Newton iteration for solving the nonlinear equation 𝜙_𝑚(𝜇)=𝜂²𝜀² is presented in Algorithm 3.

(1) Set 𝜇0=0 and 𝑖=0. (2) Solve the least-squares problem       min𝑧∈ℝ𝑚+𝑘+1𝜇𝑖1/2(𝐻𝑅−1)𝑇𝐼0𝑧−(𝑈,𝑉𝑚+1)𝑇𝑏    to obtain 𝑧𝜇𝑖,𝑚. (3) Compute 𝜙𝑚(𝜇𝑖)=𝑧𝑇𝜇𝑖,𝑚𝑧𝜇𝑖,𝑚+‖(𝐼−𝑃)𝑏‖2. (4) Solve the least-squares problem      min𝑧∈ℝ𝑚+𝑘+1𝜇𝑖1/2(𝐻𝑅−1)𝑇𝐼0𝑧−𝐻𝑅−1(𝐻𝑅−1)𝑇𝑧𝜇𝑖,𝑚    to obtain 𝑤𝜇𝑖,𝑚. (5) Compute 𝜙′𝑚(𝜇𝑖)=−2𝑧𝑇𝜇𝑖,𝑚𝑤𝜇𝑖,𝑚. (6) Compute the new approximation         𝜇𝑖+1=𝜇𝑖−𝜙𝑚(𝜇𝑖)−𝜂2𝜀2𝜙′𝑚(𝜇𝑖). (7) If |𝜇𝑖+1−𝜇𝑖|<10−6, stop; else 𝜇𝑖∶=𝜇𝑖+1, 𝑖∶=𝑖+1, and go to 2.

Algorithm 3 

Newton's method for 𝜙_𝑚(𝜇)=𝜂²𝜀².

Let 𝑥_𝑚 be the 𝑚th iterate produced by the augmented RRGMRES in [16]. Then, 𝑥_𝑚 satisfies 𝐴𝑥_𝑚−𝑏²=min_{𝑥∈𝒦𝑚()𝐴,𝐴𝑏+𝒲}‖𝐴𝑥−𝑏‖²=min_{𝑦∈ℝ𝑚+𝑘}𝐴𝑊,𝑉_𝑚𝑦−𝑏²=min_{𝑦∈ℝ𝑚+𝑘}𝑈,𝑉_𝑚+1𝐻𝑦−𝑏²=min_{𝑦∈ℝ𝑚+𝑘}𝐻𝑦−[𝑈,𝑉_𝑚+1]^𝑇𝑏²+‖(𝐼−𝑃)𝑏‖². (34) As 𝜇→∞, the reduced minimization problem (12) is the same as the above minimization problem, which shows that 𝑥_𝑚=lim_𝜇→∞𝑥_𝜇,𝑚. (35) Therefore, we obtain the following result.

Theorem 3. Let 𝑥_𝑚 be the 𝑚th iterate determined by augmented RRGMRES applied to (1) with initial iterate 𝑥₀=0. Then 𝐴𝑥_𝑚−𝑏²=𝜙_𝑚(∞). (36)

3. Numerical Experiments

In this section, we present some numerical examples to illustrate the performance of the augmented Arnoldi-Tikhonov regularization method for the solution of large-scale linear ill-posed systems. We compare the augmented Arnoldi-Tikhonov regularization method implemented by Algorithm 2 to the Arnoldi-Tikhonov regularization method proposed in [11]. The Arnoldi-Tikhonov regularization method is denoted by ATRM while the augmented Arnoldi-Tikhonov regularization method is denoted by AATRM. In all the following tables, we denote by MV the number of matrix-vector products and by RERR the relative error ‖𝑥_𝜇,𝑚−𝑥‖/‖𝑥‖, where 𝑥 is the exact solution of the linear error-free system of (2). Note that the number of matrix-vector products in Algorithm 2 is 𝑘+𝑚_min+𝑚₀. In all the examples, 𝜂=1.01.

All the numerical experiments are performed in Matlab on a PC with the usual double precision, where the floating point relative accuracy is 2.22⋅10⁻¹⁶.

Example 4. The first example considered is the Fredholm integral equation of the first kind, which takes the generic form ¹₀𝑘(𝑠,𝑡)𝑥(𝑡)𝑑𝑡=𝑏(𝑠),0≤𝑠≤1. (37) Here, both the kernel 𝑘(𝑠,𝑡) and the right-hand side 𝑏(𝑠) are known functions, while 𝑥(𝑡) is the unknown function. For test, the kernel 𝑘(𝑠,𝑡) and the right-hand side 𝑏(𝑠) are chosen as 𝑘()=𝑠()𝑡()𝑏(𝑠)=𝑠𝑠,𝑡𝑡−1,𝑠<𝑡,𝑠−1,𝑠≥𝑡,³−𝑠6. (38) With this choice, the exact solution of (37) is 𝑡. We use the Matlab program deriv2 from the regularization package [26] to discretize the integral equation (37) and to generate a system of linear equations (2) with the coefficient matrix 𝐴∈ℝ^200×200 and the solution 𝑥∈ℝ²⁰⁰. The condition number of 𝐴 is 4.9⋅10⁴. The right-hand side 𝑏 is given by 𝑏=𝐴𝑥+𝑒, where the elements of the error vector 𝑒 are generated from normal distribution with mean zero and the norm of 𝑒 is 10⁻²⋅‖𝐴𝑥‖. The augmentation subspace is one-dimensional and is spanned by the vector 𝑤=[1,2,…,200]^𝑇. Numerical results for the example are reported in Table 1 for several choice of the number 𝑚₀ of additional iterations.

Table 1 

Computational results of Example 4.

From Table 1, we can see that for 0≤𝑚₀≤3, AATRM has smaller relative errors than ATRM, and the smallest relative error is given by AATRM with 𝑚₀=1. The exact solution 𝑥 and the approximate solutions generated by AATRM and ARTM with 𝑚₀=1 are depicted in Figure 1.

Figure 1 

Example 4. Exact solution 𝑥; approximation solutions with 𝑚0=1.

Example 5. This example comes again from the regularization package [26] and is the inversion of Laplace transform ^∞₀𝑒^−𝑠𝑡𝑥(𝑠)𝑑𝑠=𝑏(𝑡),𝑡≥1, (39) where the right-hand side 𝑏(𝑡) and the exact solution 𝑥(𝑡) are given by 𝑏(𝑡)=1(𝑡)𝑡+1/2,𝑥=𝑒^−𝑡/2. (40) The system of linear equations (2) with the coefficient matrix 𝐴∈ℝ^250×250 and the solution 𝑥∈ℝ²⁵⁰ is obtained by using the Matlab program ilaplace from the regularization package [26]. In the same way as Example 4, we construct the right-hand side 𝑏 of the system of linear equations (1). The augmentation subspace is also one-dimensional and is spanned by the vector 𝑤=[1,1/2,…,1/250]^𝑇. Numerical results for the example are reported in Table 2 for 𝑚₀=1,2,3. Table 2 shows that AATRM with 𝑚₀=2 has the smallest relative error, and AARTM works slightly better than ATRM for this problem.

Table 2 

Computational results of Example 5.

Example 6. This example considered here is the same as Example 5 except that the right-hand side 𝑏(𝑡) and the exact solution 𝑥(𝑡) are given by 𝑏(𝑡)=2(𝑡+1/2)³,𝑥(𝑡)=𝑡²𝑒^−𝑡/2. (41) By using the same Matlab program as Example 5, we generate a system with 𝐴∈ℝ^1000×1000. The augmentation subspace is one-dimensional and is spanned by the vector 𝑤=[1,1,…,1]^𝑇. In Table 3, we report numerical results for 𝑚₀=1,2,3.

Table 3 

Computational results of Example 6.

We observe from Table 3 that for this example AATRM has almost the same relative errors as ATRM, and the approximate solution can be slightly improved by the augmentation space spanned by 𝑤=[1,1,…,1]^𝑇.

4. Conclusions

In this paper we propose an iterative method for solving large-scale linear ill-posed systems. The method is based on the Tikhonov regularization technique and the augmented Arnoldi technique. The augmentation subspace is a user-supplied low-dimensional subspace, which should contain a rough approximation of the desired solution. Numerical experiments show that the augmented method is effective for some practical problems.

Acknowledgments

Yiqin Lin is supported by the National Natural Science Foundation of China under Grant 10801048, the Natural Science Foundation of Hunan Province under Grant 11JJ4009, the Scientific Research Foundation of Education Bureau of Hunan Province for Outstanding Young Scholars in University under Grant 10B038, the Science and Technology Planning Project of Hunan Province under Grant 2010JT4042, and the Chinese Postdoctoral Science Foundation under Grant 2012M511386. Liang Bao is supported by the National Natural Science Foundation of China under Grants 10926150 and 11101149 and the Fundamental Research Funds for the Central Universities.