An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems

Applied Mathematical Modelling xxx (2013) xxx–xxx

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems Shervan Erfani a, Ali Tavakoli a,⇑, Davod Khojasteh Salkuyeh b a b

Mathematics Department, Vali-e-Asr University of Rafsanjan, Iran Faculty of Mathematical Sciences, University of Guilan, P.O. Box 1914, Rasht, Iran

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 13 April 2012 Received in revised form 4 February 2013 Accepted 24 March 2013 Available online xxxx Keywords: Lanczos bidiagonalization Ill-posed problems Preconditioner Regularization Singular value decomposition

Rezghi and Hosseini [M. Rezghi, S.M. Hosseini, Lanczos based preconditioner for discrete ill-posed problems, Computing 88 (2010) 79–96] presented a Lanczos based preconditioner for discrete ill-posed problems. Their preconditioner is constructed by using few steps (e.g., k) of the Lanczos bidiagonalization and corresponding computed singular values and right Lanczos vectors. In this article, we propose an efficient method to set up such preconditioner. Some numerical examples are given to show the effectiveness of the method. Ó 2013 Published by Elsevier Inc.

1. Introduction We consider large scale discrete ill-posed problems, i.e., linear systems of equations of the from

Ax ¼ b;

A 2 Rnn ;

x; b 2 Rn :

ð1Þ

The coefficient matrix A is of full-rank and typically ill-conditioned, the right-hand side vector b is perturbed by an error such ^ þ e, which may occur by measurement or discretization errors. Here e and b ^ denote the Gaussian noise and unthat b ¼ b known error-free right-hand side vectors, respectively. These problems typically arise, for instance, in discretizing linear ill-posed problems, such as Fredholm integral equation of the first kind with a smooth kernel. Also in these problems the two following criteria are satisfied: (1) The singular values of A decay gradually to zero with no particular gap in the spectrum. (2) The ratio between the largest and the smallest nonzero singular values is large [1]. ^ with error-free unknown right-hand side. We would like to obtain an approximate solution ^ x of the linear system A^ x¼b Since the matrix A is severely ill-conditioned, the solution of system (1) typically is not a significant approximation of ^ x even if the noise e is small. x. So, we use the regularization methods to determine a solution that approximates the exact solution ^ The conjugate gradient method applied to normal equations (CGLS) is a well-known iterative method for solving ill-posed problems. Since, iterative methods, such as CGLS, have low rate of convergence, it is possible to speed up the convergence by a suitable preconditioner [2]. ⇑ Corresponding author. E-mail addresses: [email protected] (S. Erfani), [email protected] (A. Tavakoli), [email protected] (D.K. Salkuyeh). URL: http://tavakoli.vru.ac.ir (A. Tavakoli). 0307-904X/$ - see front matter Ó 2013 Published by Elsevier Inc. http://dx.doi.org/10.1016/j.apm.2013.03.064

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

2

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

All regularization methods use one or more regularization parameters specific to the regularization method that controls the amount of the stabilization imposed on the solution, and in most cases it is necessary to choose this parameter from given the problem and the given set of data. A large regularization parameter makes a new well-conditioned problem, but its solution may be far from the exact solution. A small regularization parameter generally yields a solution very close to the noise-contaminated exact solution of (1), and hence its distance from the noise-free solution also can be large. Thus, one can choose a regularization parameter to balance the error due to noise with the error due to regularization. A good choice of regularization parameter is clearly crucial to determine useful approximate solution for ill-posed problems [3]. In [2], Rezghi and Hosseini proposed a new preconditioner for discrete ill-posed problems which is produced by a few steps of Lanczos bidiagonalization process. The computation of their preconditioner depends on the large singular values of A. In this paper, we propose a new and efficient numerical method to estimate the large singular values of the matrix A and construct the new form of Lanczos based preconditioner, which can be obtained in a very small number of operations. The paper is organized as follows. In Section 2, we study Fredholm integral equation of the first kind and introduce the quadrature rule and Galerkin method for discretization of these equations. Also, we study the treatment of the singular values on the matrix resulted by discretization. In Section 3, we define Lanczos bidiagonalization for constructing Least squares via Lanczos bidiagonalization (LSQR) method and preconditioner to approximate the solution of the ill-posed problems. In Section 4, we introduce a new method for constructing the new form of Lanczos based preconditioner and finding regularization parameter based on large singular values of the singular spectrum. Finally, in Section 5 some numerical experiments are given. 2. Fredholm integral equation of the first kind 2.1. The singular value expansion A classical example in linear ill-posed problems is a Fredholm integral equation of the first kind with a square integrable kernel

Z

b

Kðs; tÞf ðtÞdt ¼ gðsÞ;

c 6 s 6 d;

ð2Þ

a

where the right-hand side g and kernel K are known functions and f is an unknown solution. Fredholm equation of the first kind can be classified into categories. The first one arises if the kernel function Kðs; tÞ is smooth. In this case, Fredholm equation is very often extremely ill-conditioned, i.e., small changes in the data cause huge changes in the unknown function f. That means the solution f is extremely sensitive to small changes in gðsÞ and so special numerical methods called regularization methods are required. In the second category, the kernel kðs; tÞ is a singular function [4]. The superior analytical tool for analysis of the first-kind Fredholm integral Eq. (2) with square integrable kernels is the singular value expansion (SVE) of the kernel. By means of the SVE, any square integrable kernel K can be written as the following infinite sum [1]: 1 X

li ui ðsÞv i ðtÞ:

Kðs; tÞ ¼

i¼1

The functions ui and v i are said to be the singular functions of K. They are orthonormal with respect to the L2 inner product. The numbers li ’s are the singular values of K, which are nonnegative and can always be ordered in non-increasing order so that

l1 P l2 P l3 P    P 0 a

and, if there exists a positive real number a such that the singular values satisfy li ¼ Oði Þ, then a is called the degree of illposedness, and problems is characterized as mildly or moderately ill-posed if a 6 1 or a > 1, respectively. On the other hand, if li ¼ Oðeai Þ, then the problem is termed severely ill-posed (see [1]). The singular values satisfy the relation

kKk2 ¼

Z

Z

b

a

c

1 dX

1 X

i¼1

j¼1

li ui ðsÞv i ðtÞ

lj uj ðsÞv j ðtÞdsdt ¼

Z a

b

Z c

1 X 1 dX

li lj ui ðsÞuj ðsÞv i ðtÞv j ðtÞdsdt;

ð3Þ

i¼1 j¼1

So, we have

kKk2 ¼

1 X

l2i < 1:

i¼1

P a 1 Also, we have 1 is convergent if a > 1 and divergent if a 6 1. This shows that l2i must decay faster than i . In general, i¼1 i p12 if K is a continuous partial derivatives of order p then li is approximately Oði Þ [5,6]. So, we conclude that if the kernel K is smoother, then the singular values li ’s decay to zero rapidly. This yields an increase in the number of small singular values. In order to solve Fredholm integral equation of the first kind we have to discretize it. There are essentially two main classes of methods, namely, quadrature and Galerkin methods (see for more details [7,1,8]). Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

3

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

2.2. The singular value decomposition Let A 2 Rnn . Then, the SVD of A is a decomposition of the from

A ¼ U RV T ¼

n X ui ri v Ti ; i¼1

where U ¼ ðu1 ; u2 ; . . . ; un Þ and V ¼ ðv 1 ; v 2 ; . . . ; v n Þ are unitary matrices, i.e., U T U ¼ V T V ¼ In , R ¼ diagðr1 ; r2 ; . . . ; rn Þ has non-negative diagonal elements appearing in non-increasing order such that

and

where

r1 P r2 P    P rn : The numbers ri ; i ¼ 1; . . . ; n are the singular values of A. In continuation, we describe that an increase in the dimensions of A will only increase the number of small singular values. To this end, let l1 ; l2 ; . . . ; lk be the large singular values of the kernel K. As we cited in the last subsection, if K is a smooth function, there exist a few large singular values li that yield k  n. Moreover, the singular values ri of A, are approximations to the singular values li of the Kernel K. This property has been studied by Wing and his coworkers in the papers [9,10], and by Hansen in [11]. Then, one can write:

r1  l1 P r2  l2 P . . . P rk  lk  0: On the other hand, since li ’s are fixed, then the large singular values ri would be almost unchanged whenever the dimension of the matrix A is increased. In connection with discrete ill-posed problems, two characteristic features of SVD are very often found (1) The singular value ri decay gradually to zero with no particular gap in the spectrum. An increase of the dimensions of A will increase the number of small singular values. (2) The left and right singular vectors ui and v i tend to have more sign changes in their elements as the index i increases (i.e., when ri decreases). Both features are consequences of the fact that the SVD of A is closely related to SVE of the underlying kernel K. In fact, in the singular values ri of A denote to many cases approximations of the singular values li of K [12]. Let us consider, for example, the test problem derive2(100,1) in [4]. The kernel K of the integral Eq. (2) is Greens function for the second derivative:

 Kðs; tÞ ¼

sðt  1Þ; s 6 t; tðs  1Þ; s P t

and both integration intervals are ½0; 1. The singular values and functions are given by (see for example [12])

8 2 ; > < li ¼ ðipÞ pffiffiffi ui ðsÞ ¼ 2 sinðipsÞ; > pffiffiffi : v i ðtÞ ¼ 2 sinðiptÞ:

ð4Þ

i ¼ 1; 2; . . . ; 2

Since the singular values are proportional to i , The problem is moderately ill-posed. Fig. 1 shows that the singular values ri of discretized matrix of the test problem deriv2(100,1) with dimension 100. It is clear from this figure that the singular values of matrix A and the corresponding li of kernel K are almost the same. Since A is a full rank matrix, then its inverse is given by

A1 ¼

n X

T v i r1 i ui

i¼1

and therefore the solution of Ax ¼ b is

x ¼ A1 b ¼

n X uT b i

i¼1

ri

vi:

ð5Þ

Moreover, since there are discretization errors or linear approximation errors, the discrete ill-posed problems always contain perturbations and error components created in tension all singular vectors of A. So, we cannot compute a stabilized solution. By (5), the errors in the vector b are amplified with coefficient r1 i . If ri is close to zero then the solution is affected by the error e in the vector b,

^ þ e; b¼b

^ P kek: kbk

In this case, we have

^ þ V R1 U T e ¼ ^x þ xe : x ¼ V R1 U T b ¼ V R1 U T b

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

4

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx 0

10

σ

i

μi −1

10

−2

10

−3

10

−4

10

−5

10

−6

10

0

10

20

30

40

Fig. 1. A comparison between

50 i

60

70

80

90

100

ri and li for derive test problem.

^ and xe ¼ V R1 U T e. If the singular values decay gradually to zero, then the term xe overcomes the exact x ¼ V R1 U T b where ^ solution. Therefore, Fourier coefficients juTi bj corresponding to small singular values with lower rate tend to zero. Hence, the terms corresponding to small singular values overshadow the solution. We should use the regularization methods to subtract or eliminate the solution corresponding to small singular values. 3. Lanczos Bidiagonalization For a given matrix A 2 Rnn , the Lanczos Bidiagonalization algorithm generates two orthogonal matrices U n and V n such that U Tn AV n ¼ Bn , where Bn is a real lower bidiagonal matrix. This algorithm is described as follows:

Algorithm 1. Lanczos Bidiagonalization 1. Let b1 :¼ kbk2 ; u1 :¼ b=b1 and 2. For i ¼ 1 to n do

v 0 :¼ 0

3. pi :¼ AT ui  bi v i1 4. ai :¼ kpi k2 5. v i :¼ pi =ai 6. qi :¼ Av i  ai ui 7. biþ1 :¼ kqi k2 8. uiþ1 :¼ qi =biþ1 9. Endfor The vectors ui and v i produced by Lanczos Bidiagonalization algorithm are called the left and right Lanczos vectors, respectively. Denoting U k ¼ ½u1 ; u2 ; . . . ; uk  and V k ¼ ½v 1 ; v 2 ; . . . ; v k , the following relations can be established: (1) (2) (3) (4) (5)

b ¼ b1 u1 ¼ b1 U kþ1 e1 , AV k ¼ U kþ1 Bk , AT U kþ1 ¼ V k BTk þ akþ1 v kþ1 eTkþ1 , V Tk V k ¼ Ik , U Tkþ1 U kþ1 ¼ Ikþ1 ,

where Ij denotes the identity matrix of order j; ei is the ith unit vector and Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

0

1

a1

Bb B 2 B B Bk ¼ B B B B @

5

C C C C C: . C C .. C . ak A bkþ1

a2

..

b3

Now, suppose we want to solve

minkb  Axk2 ; x2S

where S denotes the k-dimensional subspace spanned by the first k Lanczos vectors v i ; i ¼ 1; . . . ; n. The solution which we seek is of the from xðkÞ ¼ V k yðkÞ for some vectors yðkÞ 2 Rk . Let rðkÞ ¼ b  AxðkÞ be the corresponding residual. From the above relations, we get

rðkÞ ¼ b1 u1  AV k yðkÞ ¼ U kþ1 ðb1 e1  Bk yðkÞ Þ: Since the columns of U kþ1 are orthonormal, we have

kr ðkÞ k2 ¼ kb1 e1  Bk yðkÞ k2 : In the step k of the LSQR algorithm we are going to solve

minkb1 e1  Bk yðkÞ k2 : yðkÞ

Now, we write the exact SVD of the computed Bk as

Bk ¼ Hk Ck Q Tk ¼

k X hi ci qTi : i¼1

Then we obtain

xðkÞ ¼ b1 V k

k X h1i i¼1

ci

qi ;

ð6Þ

where k is typically small. The best value of k can be computed by one of the following four ways: (1) In [2] the parameter k is considered as the smallest integer for which

rk < sr1 ;

ð7Þ

where s is the square root of the machine precision. Here, ri ’s are the singular values of Bk [2]. (2) Generalized Cross-Validation (GCV) The parameter k is chosen by GCV method that minimizes the following function (see [3]):

wk ¼

kb1 e1  Bk yðkÞ k22 ðn  kÞ2

(3) L-Curve To determine the L-curve associated with LSQR, estimates of krðkÞ k and kxðkÞ k are needed for several values of k. In the method, the corner of the L-curve gives a good balance of the solution size and residual size, so the parameter k can be chosen by the corner of this curve (see [3]). (4) Discrete Picard condition A standard tool for analyzing the discrete ill-posed problems is the discrete Picard plot, which is a plot of the quantities   T ri ; uTi b and juri ibj that arise in (5). In order to derive a meaningful regularization solution in a discrete ill-posed problem,   it must satisfies the discrete picard condition, i.e the Fourier coefficient uT b on the average should decay to zero faster i

than the singular values

ri [4]. In this paper, we use the discrete Picard condition to find the parameter k.

juT bj Fig. 2 shows plots of the first 200 singular values, Fourier coefficients juTi bj and coefficients ri of the solution for peri 3 1 turbed problem shaw [4] with Gaussian noise in which kkek and kkek ^ ¼ 1:10 ^ ¼ 1:10 . As we observe, one can choose k ¼ 8 bk bk 3 1 and k ¼ 7 for shaw problem with kkek and kkek ^ ¼ 1:10 ^ ¼ 1:10 , respectively. We moreover see that the different perturbations bk bk

contain a little change in the value of parameter k. e RV e T be the SVD of A and Now we review the preconditioner presented by Rezghi and Hosseini in [2]. Let A ¼ U e e e e e ~ ~ ~ ~ V ¼ ½ V 1 ; V 2 , where V 1 ¼ ½v 1 ; . . . ; v k  and V 2 ¼ ½v kþ1 ; . . . ; v n . In [2], it has been shown that if Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

6

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx Picard plot(a)

20

Picard plot(b)

20

10

10 σi

σi

T i T |u b|/σ i i

T i T i

|u b| 15

10

|u b| 15

10

10

i

10

10

10

5

5

k=8

10

k=7

10

0

0

10

10

−5

−5

10

10

−10

−10

10

10

−15

−15

10

10

−20

10

|u b|/σ

−20

0

1

10

2

10

10

3

10

10

0

10

i

1

2

10

10

3

10

i juT bj

3 Fig. 2. The first 200 singular values ri , Fourier coefficients juTi bj, and coefficients ri for the shaw test problem with two perturbations kkek (a) and ^ ¼ 1:10 bk i 1 kek (b). ^ ¼ 1:10 kbk

e 1 R2 V e T þ ðI  V e1V e T Þ; P¼V 1 1 1

ð8Þ

then

e PAT A ¼ V



I1

0

0

R22



eT; V

ð9Þ

where R1 ¼ diagðr1 ; . . . ; rk Þ and R2 ¼ diagðrkþ1 ; . . . ; rn Þ. Eq. (9) shows that the preconditioner P improves the first k largest singular values of AT A, those are important in reconstruction of the exact solution for discrete ill-posed problems, without affecting the small ones. According the above observation, Rezghi and Hosseini proposed a preconditioner of the form

M ¼ V k ðBTk Bk Þ1 V Tk þ ðI  V k V Tk Þ 2 Rnn ;

ð10Þ

which gives an approximation of P. In fact, from

Bk ¼ Hk C

T kQ k

¼ Hk

k C 0

! Q Tk ;

 2 Q T and substituting this in Eq. (10) gives we have ðBTk Bk Þ1 ¼ Q k C k k

 2 Q T V T þ ðI  V k V T Þ ¼ ðV k Q k ÞC  2 ðV k Q k ÞT þ ðI  ðV k Q k ÞðV k Q k ÞT Þ: M ¼ V kQ kC k k k k k  k (see [2]), we see that M gives a reasonable approximation of P. Now, since V 1  V k Q k and R1  C Eq. (10) provides a new regularized preconditioner obtained by k steps of Lanczos bidiagonalization for discrete ill-posed problems introduced in the previous sections. The proposed preconditioner clusters approximately the large singular values around 1 and leaves the others unchanged.The idea here is based on the fact that the preconditioner is constructed by k  n steps of Lanczos bidiagonalization. 4. An effective method for computing optimal Lanczos based preconditioner We know, the matrix A can be produced by discretization of Fredholm integral equation of the first kind. Solving these equations by methods like Galerkin or quadrature rules usually produce the ill-posed system Ax ¼ b that yields a huge condition number of A. The solution of the system Ax ¼ b converges well, if the dimensions of A are sufficiently large (see for Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

7

instance [7]). However, in this case, the singular values of A generated by SVD are not exact hence make the selection of parameter k difficult. In the sequel, we introduce a criterion for choosing a suitable regularization parameter k. To do so, we discretize Fredholm  Then the matrix A resulted from this kind of discretiza ¼ b. integral equation on a coarse grid that produces the system Ax tion has less dimensions than those of A (e.g, half of A) and its singular values can be absolutely computed more easily than the original matrix. On the other hand, as we mentioned in Section 2, an increase in the dimensions of A will increase the number of small singular values. Thus, the number of large singular values (i.e., k) for A and A are almost the same. Moreover the number k (for the matrix A) can be derived via the methods listed in Section 3. All we need is to compute the right Lanczos vectors V k . To do this, we first note that such vectors are available for A. Let v ¼ ðv 1 ; v 2 ; . . . ; v m ÞT be a right Lanczos vector corresponding to a coarse grid. We use an interpolation (prolongation) operator I k as follows:

I k v ¼ w; where

w2j1 :¼ v j

and w2j :¼

1 j ðv þ v jþ1 Þ 2

ð11Þ

for j ¼ 1; 2; . . . ; m. In fact, the values of coarse grid points are mapped unchanged to the fine grid and the values of the fine grid points which are not on the coarse grid, are equal to the average of the point values in their neighborhood of coarse grid. In this way, we obtain a vector of the form w ¼ ðw1 ; w2 ; . . . ; w2m1 ÞT which is corresponding to a fine grid. It should be noted that one may apply the operator I k several times (e.g., r) to obtain the right Lanczos vector for the matrix A. The following algorithm describes the procedure clearly:

Algorithm 2. Computing right Lanczos vectors of A from

v

 2 Rm ; v 2 Rm and r 2 N. 1 ;    ; a m ; b 2 Rn ; b 1. Input A ¼ ½a1 ;    ; an ; A ¼ ½a 2. For s ¼ 1 to r do 3. m :¼ lengthðv Þ 4. w :¼ zerosð2m  1; 1Þ 5. j :¼ 1 6. For i ¼ 1 to m do 7. wj :¼ v i 8. If i – m, then 9. wjþ1 ¼ ðv i þ v iþ1 Þ=2 10. j¼jþ2 11. Endif 12. Endfor v :¼ w 13. 14. EndFor 15. q :¼ 1  – 0 and aT b – 0, Tq b 16. If a q 17. 18. 19.

Tq 

v q1 ¼ kAa bbk  , T

aTq b

v q1 ¼ kA bk, T

a

v ¼ v q1 , q1

20. Else 21. q :¼ q þ 1 22. Go to Line 16 23. EndIf r

24. Return aw 2 R2 ðm1Þþ1 Remark 1. We note that 2r ðm  1Þ þ 1 ¼ n and so

m ¼ ðn  1Þ2r þ 1:

ð12Þ

Hence, in discretization of ill-posed problems like Fredholm equation of the first kind, one can choose n such that ðn  1Þ2r þ 1 is a natural number for some small r (e.g r < 10). Moreover, to find the suitable value of m, we should find the optimal value of r in (12). To do this, we select r the smallest integer value such that the right Lanczos vectors of A are exactly (not approximately) computable. Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

8

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

5. Convergence In order to show the convergence of Algorithm 2, we take simply for i ¼ 1; 2; . . . ; r, the vectors ui and v i to be the exact and approximated right Lanczos vectors of the discretized system on the i-th coarse grid, respectively. In addition, let I iþ1 be the i interpolation operator such that it takes vectors v i in the ith coarse grid and produces vectors v iþ1 in the ði þ 1Þth fine-grid according to the rule I iþ1 i v i ¼ v iþ1 . Now, we assume that

kuiþ1  I iiþ1 ui k 6 i ;

i

which

ð13Þ

denotes the accuracy of the approximation. We have:

v i ¼ I ii1 I iiþ1 . . . I rr1 v r1 ;

i ¼ 1; 2; . . . ; r:

ð14Þ

On the other hand, because of the boundedness of the interpolation operator, there exists a constant ai such that for the vectors v i and wi in the ith grid,

kI iiþ1 v i  I iiþ1 wi k ¼ kI iþ1 i ðv i  wi Þk 6 ai kv i  wi k

ð15Þ

holds. Then, by (13)–(15), we have:

kur  v r k ¼ kur  I rr1 v r1 k 6 kur  I rr1 ur1 k þ kI rr1 ur1  I rr1 v r1 k 6 r1 þ ar1 kur1  v r1 k: Similarly, we have:

kuiþ1  v iþ1 k 6 i þ ai kui  v i k;

i ¼ 1; 2; . . . ; r  1:

Hence,

kur  v r k 6 r1 þ ar1 ðr2 þ ar2 ð. . . ð1 þ a1 ku1  v 1 kÞ . . .ÞÞ: For the coarsest level 1, we assume that ku1  v 1 k 6 1 . So, kur  v r k. Let

¼

max

i¼1;...;r1

i ; a ¼

 ¼ 1 þ a1 ð2 þ a2 ð. . . ðr1 þ ar1 r Þ . . .ÞÞ is an upper bound for

max ai :

i¼1;...;r1

In this case, we see that

kur  v r k 6 Now, if

1  ar : 1a

 ! 0, then kur  v r k ! 0 and convergence is guaranteed.

Remark 2. We note that 2 kur  v r k ¼ kur  I rr1 I r1 r2 . . . I 1 v 1 k;

ð16Þ

where v 1 ¼ v . Then, by (16) we can measure the difference between the exact right Lanczos vector ur and that of the approximated vector obtained using v in Algorithm 2 for different values of m. Discussion: Without loss of generality, we can assume that u1 is the vector values of a sufficiently smooth function f1 : ½0; 1 ! R such that f1 ðxj Þ ¼ uj1 , in which

xj  xj1 ¼ h;

j ¼ 1; 2; . . .

with xj and h as a knot value and a constant value, respectively. On the other hand, the general form of the interpolation matrix operator in Algorithm 2 is as follows:

0

1 0 0 0 0 B 1=2 1=2 0 0 0 B B B 0 1 0 0 0 B B 0 1=2 1=2 0 0 B B B 0 0 1 0 0 B B 0 1=2 1=2 0 B 0 B B 0 0 0 1 0 @ .. .. .. .. .. . . . . .

...

1

...C C C ...C C ...C C C: ...C C C ...C C ...C A .. .

Here, I21 : u1 ! u2 plays the role of piecewise linear interpolation operator on the function f2 : ½0; 1 ! R such that f2 ðxj Þ ¼ uj2 , in which

xj  xj1 ¼

h ; 2

j ¼ 1; 2; . . .

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

9

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

The error estimation for this approximation is well known from any textbook (e.g. [13]) on interpolation as

ku2  I 21 u1 k 6

 2 M h  ; 8 2

where M is a constant value. Hence, in general we can have

kuiþ1  I iþ1 i ui k 6

 2 M h  ; 8 2i

But, in fact the value of

i <

 2 h

M  8 2i

i

i ¼ 1; 2; . . . ; r  1:

can be chosen such that

ð17Þ

;

because as we mentioned before the values of coarse grid points are mapped unchanged to the fine grid. In other words, the  2 odd components of the vector ui are derived exactly. Hence, it is clear that i is much less than M8  2hi . 6. Complexity

v

v

v

Let V k ¼ ½  1 ;    ;  k  with  i ; i ¼ 1;    ; k as the right Lanczos vectors for r I k V k ¼ ½I rk  1 ;    ; I rk  k , where I rk is the prolongation operator and it is applied r-times in 2, I rk  ¼ w. Now, the new form of the preconditioner can be defined as follows:

v

v

the coarsest grid. We define Algorithm 2. In fact, by Algorithm

v

MMG ¼ a2 ðI rk V k ÞðBTk Bk Þ1 ðI rk V k ÞT þ ðI  a2 ðI rk V k ÞðI rk V k ÞT Þ: We call this new preconditioner by subscript MG, because this preconditioner is based on an idea of prolongation and restric2 tion operators in Multigrid Algorithms. The construction of regularized inverse preconditioners M and M MG need Oðkn Þ and 2 Oðkm þ rmÞ operations, respectively. We know, k and r are very small in comparison with m and n. On the other hand, by (12), the value of m is much smaller than n. Therefore, computational cost of the new form Lanczos based preconditioner M MG for large n is more less than that of the preconditioner M. This preconditioner can be used as a left or right preconditioner with CGLS method [14,2]. Algorithm 3 shows the scheme of a left preconditioned CGLS method. In fact, CGLS is a semiconvergent method: For some j, in the first j iterations, the method converges to the exact solution, and then suddenly starts to diverge and the noise begins to enter the solution. Methods for finding the optimal value of such j (e.g. L-curve) can be found in [1,15].

Algorithm 3. Left-Preconditioned CGLS M MG AT Ax ¼ M MG AT b

1. Let r ð0Þ ¼ b  Axð0Þ ; pð0Þ ¼ sð0Þ ¼ MMG ðAT rð0Þ Þ; 2. For j ¼ o until convergence do 3. tðjÞ ¼ MMG pðjÞ , 4.

qðjÞ ¼ AtðjÞ ,

5.

aðjÞ ¼ kqðjÞj k2 ,

6.

xðjþ1Þ ¼ xðjÞ þ aj t ðjÞ ,

7.

rðjþ1Þ ¼ rðjÞ  aj qðjÞ ,

8.

sðjþ1Þ ¼ MMG ðAT r ðjþ1Þ Þ,

9.

cjþ1 ¼ ksðjþ1Þ k22 ,

10.

bj ¼ cjþ1 , j

c0 ¼ ksð0Þ k22

c

2

c

11. pðjþ1Þ ¼ sðjþ1Þ þ bk pðjÞ , 12. Endfor

To illustrate the efficiency of the proposed method, we consider the shaw test problem for two cases of dimensions 250 and 2000. Fig. 3 shows the singular values obtained by SVD from the discretized matrix of shaw problem, with the dimensions of 2000 and 250 for the matrices A and A, respectively. As we observe, the value and number of large singular values of A and A are almost the same. Also, Fig. 4 shows a comparison between right Lanczos vectors of shaw problem for different dimensions of A with the modifier coefficient a ¼ vv 11  1=3. 11 Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

10

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx

7. Numerical results In this section, we present some numerical examples. All experiments were carried out using MATLAB and Hansen Regularization [4]. We investigate the performance of our proposed regularization parameter (k) for introduced preconditioner and LSQR method. In the following examples, taken from [4] the contaminated vector b and the Gaussian noise e are con^ þ e, where sidered as b ¼ b

^ x e ¼ 103 kbk: kxk and x is Gaussian random vector with mean 0 and standard deviation 1. In all the following examples, we consider the coarse grid with dimensions 250  250. On this coarse grid, the parameter k is obtained by (7). Moreover, by using the interpolation operators, the ingredients of the preconditioner M, i.e. V k ; rk and Bk are obtained for fine grid. Example 7.1. Consider the integral equation

Z

1

Kðs; tÞf ðtÞdt ¼ gðsÞ;

0 6 s 6 1;

0

where

Kðs; tÞ ¼ ( f ðtÞ ¼



sðt  1Þ; s < t; tðs  1Þ; s P t; t < 12 ;

t;

1  t; t P 12 ; ( ð4s3 3sÞ

gðsÞ ¼

24

;

ð4s3 þ12s2 9sþ1Þ ; 24

s < 12 ; s P 12 :

This is a mildly ill-posed problem; i.e., its singular values decay slowly to zero. We discretize this integral equation, using the function derive2(2000,3) (see [4]), to obtain the matrix A 2 R20002000 :

5

10

σi for n=2000 σi for n=250 0

10

−5

10

−10

10

−15

10

−20

10

−25

10

0

10

1

10

2

10 i

3

10

4

10

Fig. 3. A comparison between singular values of shaw test problem for different dimensions of A.

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

11

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx 0.04

0.04 2000

v1∈ ℜ

0.03

0

v2∈ℜ

−0.02

0.01

−0.04 0

500

1000 i

1500

2000

0.06

−0.06

3

250

v ∈ℜ

0.02

500

1000 i

1500

2000

v ∈ℜ2000 4

0.02

v ∈ℜ250

3

4

0

0

−0.02

−0.02 0

500

1000 i

1500

2000

0.04

−0.04

0

500

1000 i

1500

2000

0.15 v ∈ℜ2000 5

0.02

250

v3∈ℜ

0.05

−0.02

0 0

500

1000 i

1500

2000

v ∈ℜ2000 6

0.1

0

−0.04

0

0.04 v ∈ℜ2000

0.04

−0.04

v2∈ℜ

250

v1∈ℜ

0.02

0

2000

0.02

250

−0.05

250

v6∈ℜ

0

500

1000 i

1500

2000

Fig. 4. A comparison between right Lanczos vectors of shaw test problem for different dimensions of A.

We apply the introduced methods in Section 4 for selected suitable regularization parameter on this coarse grid. In this problem, we have only used k ¼ 4 steps of Lanczos bidiagonalization for constructing the preconditioner. Fig. 5 shows the solution of Lanczos preconditioned obtained for k ¼ 4 steps of the bi-diagonalization Lanczos and un-preconditioned systems after two iteration steps. Also, Fig. 6 shows that the relative error for preconditioned linear systems with selection k ¼ 4 is better than that of the original system. Example 7.2. The inverse heat equation used here is a Volterra integral equation of the first kind with ½0; 1 as integration interval. The kernel is Kðs; tÞ ¼ Kðs  tÞ with 3

kðtÞ ¼

t2 1 pffiffiffiffi expð Þ: 2 2kappa p 4kappa t

Here, the parameter kappa controls the ill-conditioning of the matrix A: kappa = 5 gives a well-conditioned matrix, kappa = 1 gives an ill-conditioned matrix. The default is kappa = 1. The integral equation is discretized by means of simple collocation and the midpoint rule with n points (see [16,17]). An exact solution x is constructed, and then the right-hand side b is produced as b ¼ Ax. We use the code heat (n,kappa) from [4] to discretize this integral to obtain the matrix A 2 R20002000 . In this problem, we used the regularization parameter k ¼ 8. Fig. 7 compares the solution obtained of preconditioned system at 4 iterations with k ¼ 8 steps of Lanczos bidiagonalization and without preconditioner. Fig. 8 shows the relative error with and without preconditioner. Example 7.3. We discretize a Fredholm integral equation of the first kind (2) with ½ p2 ; p2  as both integration interval. The kernel K and solution f are given by

Kðs; tÞ ¼ ðcos s þ cos tÞ2



2 sin u ; u

u ¼ pðsin s þ sin tÞ;

f ðtÞ ¼ a1 expðc1 ðt  t1 Þ2 Þ þ a2 expðc2 ðt  t2 Þ2 Þ: The parameters a1 ; a2 ; c1 ; c2 ; t 1 and t2 are constants that determine the shape of the solution f. In this implementation, we use a1 ¼ 2; a2 ¼ 1; c1 ¼ 6; c2 ¼ 2; t 1 ¼ 0:8 and t2 ¼ 0:5. The kernel and the solution are discretized by simple collocation with n Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

12

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx 0.012 Exact Lanczos based preconditioner with k=4 No−preconditioner 0.01

0.008

0.006

0.004

0.002

0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Fig. 5. A comparison between the exact solution and the solution obtained by preconditioned CGLS method at step 2 and the selection k ¼ 4 of bidiagonalization Lanczos for derive2(2000,3) problem.

0.5 No−preconditioner Lanczos based preconditioner with k=4 0.45

0.4

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0

0

2

4

6

8

10

12

14

Fig. 6. Relative errors versus iteration steps for AT Ax ¼ AT b of derive2(2000,3) problem.

points to produce A and x. Then the discrete right-hand side is produced by b ¼ Ax. We used the code shaw from [1] to obtain A 2 R20002000 . In Fig. 9, we plotted the solution obtained by Lanczos based preconditioner with k ¼ 7 at step 2 and un-preconditioned systems. Fig. 10 shows, the relative errors from iteration steps of CGLS method for preconditioned systems with k ¼ 7 and original system. Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

13

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx 1.2 Exact Lanczos based preconditioner with k=8 No−preconditioner 1

0.8

0.6

0.4

0.2

0

−0.2

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Fig. 7. A comparison between the exact solution and the solution obtained by preconditioned CGLS method with 4 iterations and the selection k ¼ 8 of Lanczos bidiagonalization process for AT Ax ¼ AT b of heat problem.

1 No−preconditioner Lanczos based preconditioner with k=8 0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

0

2

4

6

8

10

12 T

14

16

18

20

T

Fig. 8. Relative errors versus iteration steps, for A Ax ¼ A b of heat problem.

Example 7.4. Define the function

 /ðxÞ ¼

1 þ cos p3x ; jxj < 3; 0;

jxj P 3:

Then kernel K, the solution f, and the right-hand side g are given by

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

14

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx 2.5 Exact Lanczos based preconditioner with k=7 No−preconditioner 2

1.5

1

0.5

0

−0.5

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Fig. 9. A comparison between the exact solution and the solution obtained by preconditioned CGLS method with 3 iterations and the selection k ¼ 7 of bidiagonalization Lanczos for AT Ax ¼ AT b of shaw problem.

1 No−preconditioner Lanczos based preconditioner with k=6 0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

1

2

3

4

5

6

7

8

9

10

Fig. 10. Relative errors versus iteration steps, for AT Ax ¼ AT b of shaw problem.

Kðs; tÞ ¼ /ðs  tÞ; f ðtÞ ¼ /ðtÞ;

  1 px 9 pjsj þ gðsÞ ¼ ð6  jsjÞ 1 þ cos sin : 2 3 2p 3

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064

15

S. Erfani et al. / Applied Mathematical Modelling xxx (2013) xxx–xxx 0.16 Exact LSQR method 0.14

0.12

0.1

0.08

0.06

0.04

0.02

0

−0.02

0

200

400

600

800

1000

1200

1400

1600

1800

2000

Fig. 11. A comparison between the exact solution and the solution obtained by LSQR method at k ¼ 6 for min kb  Axk of phillips problem.

Both integration intervals are ½6; 6. We discretize this problem, using the function phillips (2000) [3], to obtain the matrix A 2 R20002000 . In the example, we used (6) to compute the approximated solution. Fig. 11, displays an exact solution and the solution obtained by LSQR method for k ¼ 6 steps of Lanczos bidiagonalization. 8. Conclusion We have studied a recently proposed preconditioner for ill-posed linear system of equations which is based on the Lanczos bidiagonalization method. Then, we have presented an efficient method to produced such preconditioner. Finally, some numerical results have been presented to show the effectiveness of the proposed method. Acknowledgments The authors thank the anonymous referees for their valuable comments and suggestions. References [1] P.C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, Philadelphia, 1998. [2] M. Rezghi, S.M. Hosseini, Lanczos based preconditioner for discrete ill-posed problems, Computing 88 (2010) 79–96. [3] M.E. Kilmer, D.P. Oleary, Choosing regularization parameter in iterative methods for ill-posed problems, SIAM J. Matrix Anal. Appl. 22 (2001) 1204– 1221. [4] P.C. Hansen, Regularization tools: a MATLAB package for analysiss and solution of discrete ill-posed problems, Numer. Alg. 6 (1994) 1–35. [5] F.R. de Hoog, Review of Fredholm equations of the first kind, The Application and Numerical Solution of Integral Equations, Sijthoff and Noordhoff, Leyden, Netherlands, 1980. 119–134. [6] F. Smithies, The eigenvalues and singular values of integral equations, Proc. London Math. Soc. 43 (1937) 255–279. [7] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, third ed., Springer, 2008. [8] A. Tavakoli, M. Ferronato, On existence-uniqueness of the solution in a nonlinear Biot’s model, Appl. Math. Inf. Sci. 7 (2013) 333–341. [9] R.C. Allen, W.R. Boland, V. Faber, G.M. Wing, Singular values and condition numbers of Galerkin matrices arising from linear integral equations of the first kind, J. Math. Anal. Appl. 109 (1985) 564–590. [10] G.M. Wing, Condition numbers of matrices arising from the numerical solutin of linear integral equations of the first kind, 9, J. Integral Equ. 9 (Suppl.) (1985) 191–204. [11] P.C. Hansen, Computation of the singular value expansion, Computing 40 (1998) 185–199. [12] A. Bjorck, Numerical Methods for Least Square Problems, SIAM, Philadelphia, 1996. [13] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, Texts in Applied Mathematics, third ed., vol. 12, Springer-Verlag, New York, 2002. [14] D.P. Oleary, J.A. Simmons, A bidiagonalization-regularization procedure for large scale discretizations of ill-posed problems, SIAM J. Sci. Stat. Comput. 2 (1981) 474–489. [15] P.C. Hansen, T.K. Jensen, G. Rodriguez, An adaptive pruning algorithm for the discrete L-curve criterion, J. Comput. Appl. Math. 198 (2007) 483–492. [16] A.S. Carasso, Determining surface temperatures from interior observations, SIAM J. Appl. Math. 42 (1982) 558–574. [17] L. Elden, The numerical solution of a non-characteristic Cauchy problem for a parabolic equation, in: P. Deuflhart, E. Hairer (Eds.), Numerical Treatment of Inverse Problems in Differential and Integral Equations, Birkhauser, Boston, 1983.

Please cite this article in press as: S. Erfani et al., An efficient method to set up a Lanczos based preconditioner for discrete ill-posed problems, Appl. Math. Modell. (2013), http://dx.doi.org/10.1016/j.apm.2013.03.064