Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 548487, 8 pages

http://dx.doi.org/10.1155/2013/548487

Research Article

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

Yiqin Lin,¹ Liang Bao,² and Yanhua Cao³

¹Department of Mathematics and Computational Science, Hunan University of Science and Engineering, Yongzhou 425100, China

²Department of Mathematics, East China University of Science and Technology, Shanghai 200237, China

³Department of Mathematics, North China Electric Power University, Beijing 102206, China

Received 1 November 2012; Revised 18 March 2013; Accepted 19 March 2013

Academic Editor: Hung Nguyen-Xuan

Copyright © 2013 Yiqin Lin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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 ,

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

where 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

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

where and 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:

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

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

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

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

where , , has orthonormal columns, and 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

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) Compute QR-factorization .

    (2) Compute , , and set .

    (3) For Do:

             ;

             For Do:

                    ;

                    ;

             End For

             For Do:

                    ;

                    ;

             End For

             ;

             ;

             ;

           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

in which is an orthogonal projector onto . Obviously, the reduced minimization problem (12) is equivalent to

The normal equations of the minimization problem (13) is

We denote the solution of the minimization problem (13) by . Then, from (14) it follows that

The approximate solution of (3) is

Since the matrix has a larger condition number than the matrix , we apply the QR factorization of to obtain the solution of (13) instead. The QR factorization of can be implemented by a sequence of Givens rotations.

Define

Substituting into (17) and using the modified Arnoldi decomposition (10) yield

Substituting in (15) into (18), we obtain

Note that

Therefore, the function can be expressed as

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

Assume that and . Then is strictly decreasing and convex for with . Moreover, the equation

has a unique solution , such that , for any with

where 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 defined by

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

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 , of the modified Arnoldi process is chosen so that

We can improve the quality of the computed solution by choosing the practical iterative step number somewhat larger than .

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: , , , , , .

Output: , , .

    (1) Compute the modified Arnoldi decomposition , with ,

         where is the smallest number such that .

    (2) Compute the solution of the equation 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 . Let

that is, satisfies the system of linear equations

Note that the above linear system is the normal equations of the least-squares problem

For numerical stability, we compute the vector by solving the least-squares problem. Then, the is evaluated by computing

It is easy to show that the first derivative of can be written as

where

Hence, we may compute by solving a least-squares problem analogous to the above with the vector replaced by .

The algorithm for implementing the Newton iteration for solving the nonlinear equation is presented in Algorithm 3.

(1) Set and .

(2) Solve the least-squares problem

      to obtain .

(3) Compute .

(4) Solve the least-squares problem

      to obtain .

(5) Compute .

(6) Compute the new approximation

                                .

(7) If , stop; else , , and go to 2.

Algorithm 3: Newton's method for .

Let be the th iterate produced by the augmented RRGMRES in [16]. Then, satisfies

As , the reduced minimization problem (12) is the same as the above minimization problem, which shows that

Therefore, we obtain the following result.

Theorem 3. Let be the th iterate determined by augmented RRGMRES applied to (1) with initial iterate . Then

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 . In all the examples, .

All the numerical experiments are performed in Matlab on a PC with the usual double precision, where the floating point relative accuracy is .

Example 4. The first example considered is the Fredholm integral equation of the first kind, which takes the generic form

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

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 and the solution . The condition number of is . 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 . The augmentation subspace is one-dimensional and is spanned by the vector . 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.

Method

MV
RERR

AATRM
0
2

AATRM
1
3

AATRM
2
4

AATRM
3
5

ATRM
0
4

ATRM
1
5

ATRM
2
6

ATRM
3
7

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

Figure 1: Example 4. Exact solution ; approximation solutions with .

Example 5. This example comes again from the regularization package [26] and is the inversion of Laplace transform

where the right-hand side and the exact solution are given by

The system of linear equations (2) with the coefficient matrix 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 . Numerical results for the example are reported in Table 2 for . Table 2 shows that AATRM with has the smallest relative error, and AARTM works slightly better than ATRM for this problem.

Table 2: Computational results of Example 5.

Method

MV
RERR

AATRM
1
5

AATRM
2
6

AATRM
3
7

ATRM
1
6

ATRM
2
7

ATRM
3
8

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

By using the same Matlab program as Example 5, we generate a system with . The augmentation subspace is one-dimensional and is spanned by the vector . In Table 3, we report numerical results for .

Table 3: Computational results of Example 6.

Method

MV
RERR

AATRM
1
10

AATRM
2
11

AATRM
3
12

ATRM
1
10

ATRM
2
11

ATRM
3
12

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 .

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.

References

H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, vol. 375 of Mathematics and Its Applications, Kluwer Academic, Dordrecht, The Netherlands, 1996.

A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems, John Wiley & Sons, New York, NY, USA, 1977.

A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov, and A. G. Yagola, Numerical Methods for the Solution of Ill-Posed Problems, vol. 328 of Mathematics and Its Applications, Kluwer Academic, Dordrecht, The Netherlands, 1995.

P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM Monographs on Mathematical Modeling and Computation, SIAM, Philadelphia, Pa, USA, 1998.

G. H. Golub and C. F. van Loan, Matrix Computations, Johns Hopkins Studies in the Mathematical Sciences, Johns Hopkins University Press, Baltimore, Md, USA, 3rd edition, 1996.

Å. Björck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia, Pa, USA, 1996.

Å. Björck, “A bidiagonalization algorithm for solving large and sparse ill-posed systems of linear equations,” BIT: Numerical Mathematics, vol. 28, no. 3, pp. 659–670, 1988.

D. Calvetti, S. Morigi, L. Reichel, and F. Sgallari, “Tikhonov regularization and the -curve for large discrete ill-posed problems,” Journal of Computational and Applied Mathematics, vol. 123, no. 1-2, pp. 423–446, 2000.

D. Calvetti and L. Reichel, “Tikhonov regularization of large linear problems,” BIT: Numerical Mathematics, vol. 43, no. 2, pp. 263–283, 2003.

D. P. O'Leary and J. A. Simmons, “A bidiagonalization-regularization procedure for large scale discretizations of ill-posed problems,” SIAM Journal on Scientific and Statistical Computing, vol. 2, no. 4, pp. 474–489, 1981.

B. Lewis and L. Reichel, “Arnoldi-Tikhonov regularization methods,” Journal of Computational and Applied Mathematics, vol. 226, no. 1, pp. 92–102, 2009.

T. F. Chan and K. R. Jackson, “Nonlinearly preconditioned Krylov subspace methods for discrete Newton algorithms,” SIAM Journal on Scientific and Statistical Computing, vol. 5, no. 3, pp. 533–542, 1984.

L. Reichel, F. Sgallari, and Q. Ye, “Tikhonov regularization based on generalized Krylov subspace methods,” Applied Numerical Mathematics, vol. 62, no. 9, pp. 1215–1228, 2012.

D. Calvetti, B. Lewis, and L. Reichel, “On the choice of subspace for iterative methods for linear discrete ill-posed problems,” International Journal of Applied Mathematics and Computer Science, vol. 11, no. 5, pp. 1069–1092, 2001, Numerical analysis and systems theory (Perpignan, 2000)

L. Reichel and Q. Ye, “Breakdown-free GMRES for singular systems,” SIAM Journal on Matrix Analysis and Applications, vol. 26, no. 4, pp. 1001–1021, 2005.

J. Baglama and L. Reichel, “Augmented GMRES-type methods,” Numerical Linear Algebra with Applications, vol. 14, no. 4, pp. 337–350, 2007.

K. Morikuni, L. Reichel, and K. Hayami, “FGMRES for linear discrete problems,” Tech. Rep., NII, 2012.

D. Calvetti, B. Lewis, and L. Reichel, “On the regularizing properties of the GMRES method,” Numerische Mathematik, vol. 91, no. 4, pp. 605–625, 2002.

R. B. Morgan, “A restarted GMRES method augmented with eigenvectors,” SIAM Journal on Matrix Analysis and Applications, vol. 16, no. 4, pp. 1154–1171, 1995.

Y. Saad, “A flexible inner-outer preconditioned GMRES algorithm,” SIAM Journal on Scientific Computing, vol. 14, no. 2, pp. 461–469, 1993.

B. N. Parlett, The Symmetric Eigenvalue Problem, vol. 20 of Classics in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1998.

C. T. Kelley, Solving Nonlinear Equations with Newton's Method, vol. 1 of Fundamentals of Algorithms, SIAM, Philadelphia, Pa, USA, 2003.

J. A. Ezquerro and M. A. Hernández, “On a convex acceleration of Newton's method,” Journal of Optimization Theory and Applications, vol. 100, no. 2, pp. 311–326, 1999.

W. Gander, “On Halley's iteration method,” The American Mathematical Monthly, vol. 92, no. 2, pp. 131–134, 1985.

L. Reichel and A. Shyshkov, “A new zero-finder for Tikhonov regularization,” BIT: Numerical Mathematics, vol. 48, no. 3, pp. 627–643, 2008.

P. C. Hansen, “Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems,” Numerical Algorithms, vol. 6, no. 1-2, pp. 1–35, 1994.