王彦飞---中国科学院地质与地球物理研究所.pdf
王彦飞---中国科学院地质与地球物理研究所.pdf
On the Regularity of a Trust Region-CG Algorithm for Nonlinear Ill-posed Inverse Problems Yan-fei Wang Ya-xiang Yuan State Key Laboratory of Scientic and Engineering Computing, Institute of Computational Mathematics and Scientic/Engineering computing, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, P.O.BOX 2719, Beijing, 100080, China Email: wyf@lsec.cc.ac.cn yyx@lsec.cc.ac.cn Abstract. In this paper we consider the regularity of the trust region-cg algorithm, when it is applied to nonlinear ill-posed iverse problems. The trust region algorithm can be viewed as a regularization method, but it diers from the traditional regularization method, because no penalty term is need. Thus, the determing of the so-called regularization parameter in a standard regularization method is avoided. Theoretical analysis of the trust region-cg method is presented, convergency and regularity of the trust region algorithm are proved, and numerical tests are also given. Key Words. Nonlinear ill-posed problems, trust region-cg method, convergence, regularity AMS Subject Classications: 65J15, 65J20, 65K10 1 Introduction In scienti c and engineering computing, we are often encountered with nonlinear inverse problems. An inverse problem consists of a direct problem and some unknown function(s) or parameters. Inverse problems are usually ill-posed in the sense of J. Hadmard, i.e., at leat one term of the existence, uniqueness, stability of the solution is vilated. Particularly we are concerned with the stability, since in many applications the solution does not depend continuously on the unknown quantities and the problem is ill-posed. A typical ill-posed problems is to determine these unknowns given measured, or condaminated data. We can outline the nonlinear ill-posed problems into an abstract operator equations F (x) = y (1) where F : D(F ) X ;! Y is a nonlinear mapping, X and Y are both seperable Hilbert spaces. We assume that F is continuous and compact for xed x 2 D(F ). Problem (1) is typically ill-posed in the sense that a solution x+ does not depend continuously on the obervation data y. Since in practice only approximate data with some error level Partially supported by Chinese NSF grant 19731010 and the Knowledge Innovation Program of CAS. The rst author also partially supported by Hebei Province NSF grant 101001. 1 2 Y.Wang, Y.Yuan , i.e., ky ; yk (2) are available, problem (1) has to be regularized (see e.g. 3, 7, 23]). Through out this paper we assume that a solution x+ of (1) exists, i.e. F (x+ ) = y: (3) Regularization methods are such kind of methods which replace the ill-posed problem with a stabilized problem whose solution depends on a parameter, named as the regularization parameter. The regularized problem is well-posed in the sense of J. Hadamard. For a complete theoretical analysis of such kind of method, please see some well-written books 23, 3, 13, 7, 14, 17]. Certainly the most well-known and most widely used regularization method for nonlinear illposed problems is the method of Tikhonov regularization. In which one solves the unconstrained minimization problem min J x y] := kF (x) ; y k2 + (x): x2X (4) > 0 is the regularization parameter. (x) serves as the stabilizer, i.e., stablizes the minimiza- tion process and provides a priori information about the solution. Replacing F (x) by rst order Taylor's expansion, i.e., (4) turns into min J y] := ky ; F (xk ) ; F 0 (xk ) k2 + ( ): 2X (5) If an approximate solution k of (5) is computed, we can let xk+1 = xk + k . Assume that xk is some approximation of the solution x+ , then F (x+ ) ; F (xk ) = F 0 (xk )(x+ ; xk ) + r(x+ xk ) (6) where r(x+ xk ) is the Taylor remainder. Denoting + = x+ ; xk and solving for it leads to F 0 (xk ) + = y ; F (xk ) ; r(x+ xk ): (7) The above relation can be rewritten as F 0 (xk ) + = y ; F (xk ) + y ; y ; r(x+ xk ): (8) Thus, the linearized problem F 0 (xk ) = y ; F (xk ) (9) is an approximation to the original problem up to up to an error err = y ; y + r(x+ xk ) with kerrk + kr(x+ xk )k: Clearly (5) is equivalent to applying Tikhonov regularization to problem (9). (10) Regularity of Trust Region-CG 3 Assumption 1.1 Assume that after k iterations, + = x+ ; xk satises ky ; F (xk ) ; F 0 (xk ) + k !ky ; F (xk )k 0 < ! < 1: Apart from the above analysis, we need the following assumption: Assumption 1.2 For a certain ball B D(F ) around the exact solution x+ of (1), and some 1 > d > 0 let kF (x) ; F (^x) ; F 0 (^x)(x ; x^)k dkF (x) ; F (^x)k (11) for all x x^ 2 B. This assumption is helpful for analyzing the properties of the trust region algorithm which we presented in this paper. Recently, optimization methods are becoming popular for solving nonlinear ill-posed inverse problems, for example, Gauss-Newton method (1, 12]), Broyden's method (11]), and Levenberg-Marquardt method (9]), which have been well developed in nonlinear programming. Trust region method has been used in parameter identi cation problem and image restoration problem (see 25, 26]) and seems promising. This paper will consider trust region method for nonlinear ill-posed inverse problems. 2 A Trust Region-CG Algorithm Considering the unconstrained optimization problem min J x y ] := kF (x) ; y k2: x2X (12) We denote by g(x) the gradient of the functional J , Hess(x) the approximate Hessian of J , i.e., g(x) = F 0 (x)T (F (x) ; y ) Hess(x) = F 0 (x)T F 0 (x): At the k-th iteration, a trust region subproblem (TRS) for (12) is T + 1 (Hessk ) := k ( ) min g k n x2R 2 s: t: k k (13) (14) where gk = g(xk ) Hessk = Hess(xk ) and > 0 is the trust region bound. (13){(14) is solved exactly or inexactly to obtain a trial step k . The ratio Aredk rk = Pred (15) Aredk = J xk y ] ; J xk + k y ] (16) k is used to decide whether the trial step k is acceptable and to adjust the trust region bound. 4 Y.Wang, Y.Yuan is called the actual reduction in the objective model, and Predk = k (0) ; k (k ) (17) is the predicted reduction. We outline the general trust region algorithm for unconstrained optimization as follows. Algorithm 2.1 (Trust region algorithm for nonlinear ill-posed problem) STEP 1 Given the initial guess value x1 2 Rn , 1 > 0, 0 < 3 < 4 < 1 < 1 , 0 0 2 < 1, 2 > 0, k := 1 STEP 2 If some stopping rule is satised then STOP Else, solve (13)-(14) giving k STEP 3 Compute rk xk+1 = Choose k+1 that satises k+1 = xk if rk 0 xk + k otherwise (18) 3 kk k 4 k ] if rk < 2 k 1 k ] otherwise (19) STEP 4 Evaluate gk and Hessk k:=k+1 GOTO STEP 2. The constant i (i = 0 4) can be chosen by users. Typical values are 0 = 0 1 = 2 2 = 3 = 0:25 4 = 0:5. For other choices of those constants, please see 4], 5], 15], 19], etc.. The parameter 0 is usually zero (see 4], 20]) or a small positive constant (see 2] and 21]). The advantage of using zero 0 is that a trial step is accepted whenever the objective function is reduced. Hence it would not throw away a \good point", which is a desirable property especially when the function evaluations are very expensive. In STEP 2, the stopping rule is based on some kind of so-called discrepancy principle, i.e., once the inequality kF (xk ) ; y k !~ with !~ > 1 is satis ed, no further iteration is needed. The following lemma is well known (for example, see 16] and 6]): Lemma 2.2 A vector 2 Rn is a solution of (13)-(14) if and only if there exists 0 such that (Hessk + I ) = ;gk (20) and that Hessk + I is positive semi-denite, k k and ( ; k k) = 0: (21) Regularity of Trust Region-CG 5 It is shown by Powell 20] that trust region algorithms for (12) is convergent if the trust region step satis es Pred( ) ckgk minf kgk=kHesskg (22) and some other conditions on Hess are satis ed. It is easy to see that 1 kgk minf kgk=kHesskg: (23) (0) ; kkmin ( s ) 2 2spanfgg Therefore it is quite common that in practice the trial step at each iteration of a trust region method is computed by solving the TRS (13)-(14) inexactly. One way to compute an inexact solution of (13)-(14) was the truncated conjugate gradient method proposed by Toint 24] and Steihaug 22] and analyzed by Yuan 30]. The conjugate gradient method for (13) generates a sequence as follows: l+1 = l + l dl (24) dl+1 = ;gl+1 + l dl (25) where gl = r k (l ) = Hessk l + gk with gk = g(xk ) = F 0 (xk )T (F (xk ) ; y ) Hessk = Hess(xk ) = F 0 (xk )T F 0 (xk ) and l = ;gl T dl =dTl Hessk dl l = kgl+1 k2=kgl k2 (26) with the initial values 1 = 0 d1 = ;g1 = ;gk . Toint 24] and Steihaug 22] were the rst to use the conjugate gradient method to solve the general trust region subproblem (13)-(14). Even without assuming the positive de niteness of Hess, we can continue the conjugate gradient method provided that dTl Hess dl is positive. If the iterate l + l dl computed is in the trust region ball, it can be accepted, and the conjugate gradient iterates can be continued to the next iteration. Whenever dTl Hess dl is not positive or l + l dl is outside the trust region, we can take the longest step along dl within the trust region and terminate the calculations. Algorithm 2.3 (Truncated conjugate gradient method for TRS) STEP 1 Given 1 = 0, 0 < < 1, (tolerance) > 0 and compute g1 = r (1 ), set l := 1 d1 = ;g1 = ;gk STEP 2 If kAk l ; u~k k ku~k k, stop, output = l Compute dTl Hess dl: if dTl Hessk dl 0 then goto step 4 Calculate l by (26). STEP 3 If kl + l dl k k then goto step 4 Set l+1 by (24) and gl+1 = gl + l Hessk dl Compute l by (26) and set dl+1 by (25) l:=l+1, goto step 2. 6 Y.Wang, Y.Yuan STEP 4 Compute l 0 satisfying kl + l dl k = Set = l + l dl , and stop. Note that l can be computed by choosing the positive root of the quadratic equation in : kdl k22 + 2(l dl ) + kl k2 ; 2k = 0: (27) Let be the inexact solution of (13)-(14) obtained by the above truncated CG method and ^ be the exact solution of (13)-(14). Recently Yuan 30] shows that (0) ; ( ) 1 (28) (0) ; (^) 2 which can be written as the following theorem: Theorem 2.4 For any > 0 g 2 Rn and any positive denite matrix Hess 2 Rnn, let s^ be the global solution of the trust region subproblem (13)-(14), and let be the solution obtained by the truncated CG method, then ( ) 21 (^s): (29) This theorem tells us that the reduction in the approximate model is at least half of the maximum reduction if we use the truncated conjugate gradient method for solving the subproblem (13)-(14). Applying Algorithm 2.3 to compute the trial step k in Step 2 of Algorithm 2.1, we obtain a trust region-cg algorithm for nonlinear ill-posed inverse problems. The algorithm consists of two stage iterations: the inner loop and the outer loop. The inner loop is the truncated conjugate gradient method, the outer loop is the trust region method. To avoid too many inner loop iterations in one out loop iteration, we terminate the inner loop iteration if itermax cg steps have been taken, where itermax is a given positive number. We also terminate the inner look iteration if a progress in function reduction in the cg step is smaller than . However, in our numerical tests, this termination rule was not activated. 3 Properties of Algorithm 2.3 In this section we give some properties of the truncated conjugate gradient method. The main result is the monotonicity of the iterates. First, we present an equivalent form of the conjugate gradient method. Denote Ak = F 0 (xk ) uk = y ; F (xk ) ; r(x+ xk ) u~k = y ; F (xk ) we have Hessk = Ak Ak gk = ;Ak u~k gl = Ak (Ak l ; u~k ) Regularity of Trust Region-CG 7 dl+1 = ;gl+1 + l dl = ;Ak (Ak l+1 ; u~k ) + l dl = ;Ak Ak l ; l Ak Ak dl + Ak u~k + l dl : If we let dl = Ak zl provided that such zl exists, then dl+1 = Ak (~uk ; Ak l ; l Ak dl + l zl ): Further if we denote rl = u~k ; Ak l , then clearly rl+1 = u~k ; Ak l+1 = rl ; l Ak dl and dl+1 = Ak (rl ; l Ak dl + l zl ) hold. We can generate the next search direction by dl+1 = Ak zl+1 with zl+1 = rl ; l Ak dl + l zl = rl+1 + l zl : Hence, the conjugate gradient iterates can be genereated in the following way: l+1 dl zl+1 rl+1 l + l dl Ak zl rl+1 + l zl rl ; l Ak dl T d 2 g l l = ; dT Hessl d = kkAAkdrl kk2 k l l l 2 k rl+1 k l = kA kAk rl k2 l+1 = 1 + l l = = = = (30) (31) (32) (33) (34) (35) (36) The initial values are 1 = 0 d1 = ;gk z1 = r1 r1 = u~k = y ; F (xk ) 1 = 1. Here another scalar l is added, which will be used for the analysis of the truncated conjugate gradient method. One tool for the analysis of the truncated conjugate gradient method is the so-called residual polynomials (see 10, 3]). Let l be the set of all polynomials of degree l or less, and set 0l := fp 2 l : p(0) = 1g: Then there is an 1-1 relation between elements 2 Kl (Ak u~k Ak Ak ) and p 2 0l via the representation u~k ; Ak = p(Ak Ak )~uk of the corresponding residual, where Kl (Ak u~k Ak Ak ) is the l-th Krylov subspace Kl (Ak u~k Ak Ak ) = spanfAk u~k (Ak Ak )Ak u~k (Ak Ak )l;1 Ak u~k g: (37) 8 Y.Wang, Y.Yuan For simplicity, pl 2 0l denotes the residual polynomial associated with l , the l-th CG iterate. The bilinear form < p q >:= (p(Ak Ak )~uk q(Ak Ak )~uk ) de nes the inner product for p q 2 l . If q 2 l;1 is an arbitrary polynomial of degree l ; 1 then the polynomial p given by p( ) = pl ( ) + t q( ) belongs to 0l for every t 2 R. Noticing that pl solves the minimization problem < p p >;! min for p 2 0l : Hence, < pl q >= 21 dtd < p p > jt=0 = 0 for all q 2 l;1 : (38) < pl 1 >=< pl pl > (39) If we de ne q = ql;1 by pl = 1 ; ql;1 , then clearly we have that, which will be used for later analysis. From (30)-(36) and the above de nitions, we have that zl = sl (Ak Ak )~uk sl ( ) := pl ( ) ; pl+1 ( ) 2 l : l (40) It was pointed out by 8] that in general sl does not belong to 0l . Instead, since the vectors zl are updated by zl+1 = rl+1 + l zl with rl+1 = u~k ; Ak l+1 , it follows from (40) that sl+1 ( ) = pl+1 ( ) + l sl ( ) (41) and hence, sl (0) and l of Algorithm 2.3 share the same recurrence relation (this in fact has been observed by 8]), i.e., sl (0) = l : (42) With the above analysis, we can now present the monotonicity of the iterates for perturbed right-hand side. Theorem 3.1 Let 2 l 2 N . If Assumption 1.1 holds and if ku~ ; A k2 + ku~ ; A k2 > ! kzl kku~k k 0 < ! < 1 l = 1 2 l k k l k k l+1 l (43) then k + ; l k is strictly monotonically decreasing for l = 1 2 l , and l X + 2 + 2 k k ; k ; l +1 k > ( ; 2)!ku~k k l kzlk: l=1 (44) Regularity of Trust Region-CG 9 Proof. By induction, we obtain k + ; l+1 k2 = k + ; l ; l dl k2 = k + ; l k2 ; 2( + ; l l Ak zl ) + (l Ak zl l Ak zl ) = k + ; l k2 ; l (2Ak + ; 2Ak l ; l Ak dl zl ) = k + ; l k2 ; l (~uk ; Ak l zl ) ; l (~uk ; Ak l+1 zl ) +2l(~uk ; Ak + zl ): From the de nitions of pl and zl , we have k + ; l k2 ; k + ; l+1 k2 = l < pl sl >+l< pl+1 sl >;2l< u~k ; Ak + zl > : By (42), sl ( ) = l + q( ) for some polynomial q 2 l;1 , and hence from (38) and (39) we nd that k + ; l k2 ; k + ; l+1 k2 = l l (< pl 1 >+< pl+1 1 >) ; 2l < u~k ; Ak + zl > = l l (< pl pl >+< pl+1 pl+1 >) ; 2l < u~k ; Ak + zl > : Since < pl pl >= ku~k ; Ak l k2 , hence it follows from the above relation and (43) that k + ; l k2 ; k + ; l+1 k2 > !l kzlkku~k k ; 2!lku~k kkzlk for all l = 1 l: (45) Thus, the monotonicity of k + ; l k follows from the above inequality and the assumption > 2. Relation (44) follows by taking the sum of (45) for l = 1 l and observing the fact that 1 = 0. Q.E.D Remark 3.2 We have noted that in general sl will not belong to 0l , however, s^l := sl=l 2 0l . Hence from the minimization property of the truncated CG we obtain ku~k ; Ak l+1 k ku~k ; Ak l k =< pl pl > 21 < s^l s^l > 21 = 1 < sl sl > 12 = 1 kzlk: l l (46) This together with (43) yield that 2ku~k ; Ak l k2 ku~k ; Ak l k2 + ku~k ; Ak l +1 k2 > ! kzl kku~k k > ! ku~k kku~k ; Ak l k l which indicates that ku~k ; Ak l k ! k2u~k k > !ku~k k: Since ku~k ; Ak + k !ku~k k according to Assumption 1.1, this shows that l can not exceed lk , the smallest index of the inner iteration. 10 Y.Wang, Y.Yuan 4 Convergence of Trust Region-CG for Exact Data Before presenting the proposition in the following paragraph, we rst give an assumption: Assumption 4.1 Assume that in each inner iteration, l satises ku~k ; Ak l k !1 ku~k k with !12 = ! (47) until convergence, where 0 < ! < ;1, > 2. Assumption 4.1 is closely related with the termination rule of Algorithm 2.3. Where, in STEP 2 serves as the number !1 here. Once the opposite inequality of Assumption 4.1 is satis ed, we stop the inner iteration. Proposition 4.2 Suppose that Assumption 4.1 holds. The inequality (47) indicates that (43) is true if l > 0. Furthermore, there are only nitely many l for which (47) holds. Proof. From the Remark 3.2 we know 1 kz k =< s^ s^ > 12 : l l l l It is proved by 10] that < s^l s^l > is strictly monotonically decreasing with l, consequently 1 kz k < 1 kz k = kr k = ku~ k: l l 1 1 1 k The above inequality together and (47) indicates that (43) is true for l > 0. From Remark 3.2, we see that (47) holds for only nitely many indices l. Q.E.D Now assume that y = y, we rst prove the monotonicity of the trust region-cg algorithm, i.e., xk + lk is a better approximation of x+ than xk . We also assumed that F (x) = y has a solution x+ 2 B D(F ). Proposition 4.3 Assume that Assumptions 1.2 and 4.1 holds, then the iteration error kx+ ; xk k is monotonically decreasing. Proof. According to Assumption 1.2, (11) holds for x = x+ x^ = xk , i.e. ky ; F (xk ) ; F 0 (xk )(x+ ; xk )k dky ; F (xk )k: Note that y = y, the above expression indicates that Assumption 1.1 is ful lled with 0 < d < 1. Due to Assumption 4.1, ku~k ; Ak l k2 !ku~k k2 is satis ed. Hence from Proposition 4.2, the requirement of Theorem 3.1 is ful lled. From our notations, we have that xk+1 = xk + lk . Thus, from + = x+ ; xk , + ; lk = x+ ; xk+1 and Theorem 3.1, we see that kx+ ; xk+1 k < kx+ ; xk k. Q.E.D Regularity of Trust Region-CG 11 Remark 4.4 Propostion 4.3 also implies two inequalities: and kuk kkwk k < d( 1; 2) (kx+ ; xk k2 ; kx+ ; xk+1 k2 ) (48) 2 kuk k2 < dk(Ak;k2) (kx+ ; xk k2 ; kx+ ; xk+1 k2 ): (49) (48) is straightforward as kuk ; Ak + k dkuk k with 0 < d < 1. (49) follows from relations kuk kkwk k > 1 kuk kkz1k = 1 kuk k2 and 2 2 1 = kkAAkdr1 kk2 = kAkAAk u kuk k2 kAk k;2: k 1 k k k Theorem 4.5 Given the exact data y = y and suppose that Assumptions 1.2 and 4.1 hold. Then the iterates fxk g generated by Algorithms 2.1 and 2.3 converge to a solution of (1) as k ! 1. Proof. First we prove that fxk g forms a Cauchy sequence. Let us denote the iteration errors by ek = x+ ; xk . Given k j 2 N with k > j , let 2 fj kg be chosen in such a way that ky ; F (x )k ky ; F (xi )k i = j k: Consider now ke ; ej k2 = kej k2 ; ke k2 + 2(e e ; ej ): Note that ej ; e = x ; xj x ; xj = P ;1 z . Hence we obtain where Ai = F 0 (xi ) wi = lli=1 l l j(e e ; ej )j = j X ;1 i=j X ;1 X ;1 i=1 i=1 (Ai wi e )j Ai wi kwi kkAi e k: We can estimate that kAi e k = kAi ei ; Ai (e ; ei )k ky ; F (xi ) ; F 0 (xi )ei k + kF (x ) ; F (xi ) ; F 0 (xi )(e ; ei )k +ky ; F (x )k dky ; F (xi )k + dkF (x ) ; F (xi )k + ky ; F (x )k 2dky ; F (xi )k + (1 + d)ky ; F (x )k (1 + 3d)ky ; F (xi )k: (50) 12 Y.Wang, Y.Yuan Relation (48) and the above expression give that j(e e ; ej )j (1 + 3d) X ;1 i=j kwi kky ; F (xi )k d1(+;3d2) (kx+ ; xj k2 ; kx+ ; x k2) which, together with (50), yields ke ; ej k2 C (kx+ ; xj k2 ; kx+ ; x k2 ) d) with C = 2(1+3 d( ;2) + 1 independent of j k . Similarly one can obtain kej ; e k2 C (kx+ ; x k2 ; kx+ ; xj k2) hence kxk ; xj k2 = kej ; ek k2 2(kek ; e k2 + ke ; ej k2 ) 2C (kx+ ; xj k2 ; kx+ ; xk k2 ): (51) Therefore, fxk g form a Cauchy sequence because the monotonicity of fkx+ ; xk kg. P 2 Denote the limit of xk by x. From (49) we know 1 k=1 kuk k converges, and therefore F (xk ) ! y as k ! 1. This indicates that x is a solution of (1). Q.E.D 5 Regularity of the Algorithm for Inexact Data Now we consider the case where inexact date y instead of y. It is assumed that (2) is satis ed. Our stopping rule is based on the discrepancy principle, i.e., we terminate the calculations at the smallest iteration index kD such that the discrepancy inequality ky ; F (xk )k !~ with !~ > 1 (52) holds. We denote xk the corresponding iterates and consider the regularity of the trust region-cg algorithm. Theorem 5.1 Assume that Assumptions 1.2 and 4.1 hold. Let x be a solution of (1) with F satises (11) for some 1 > d > 0 in a ball B D(F ) around x. Let !~ in (52) be chosen that d !~ > 1+ 1;d . Then kx ; xk k is monotonically decreasing. Moreover, Algorithm 2.1 terminates after kD < 1 iterations. Proof. We prove that kx ; xk+1 k kx ; xk k (53) Regularity of Trust Region-CG 13 with x a solution of (1). Using Assumption 1.2, we estimate that ky ; F (xk ; F 0 (xk )(x ; xk )k + kF (x) ; F (xk ; F 0 (xk )(x ; xk )k + dky ; F (xk )k (1 + d) + dky ; F (xk )k: According to the discrepancy principle, ky ; F (xk )k > !~ as k < kD , hence < !~1 ky ; F (xk )k and ky ; F (xk ; F 0 (xk )(x ; xk )k 1 + d!~ + !~ ky ; F (xk )k: By assumption, 0 < (1 + d + !~ )=!~ < 1, hence Assumption 1.1 is ful lled. Consequently Propostion 4.3 applies and the monotonicity assertion (53) follows as in the proof of Proposition 4.3. Next we show that there are only nite number of iterations. In fact as the same as in the proof of (49), we have ky ; F (xk )k2 < d(L; 2) (kx ; xk k2 ; kx ; xk+1 k2 ) (54) with L = supfkF 0 (xk )k2 g for all k < kD . Assume that (53) holds for x = x+ . Now taking the sum of (54) for k = 1 2 kD ; 1 we obtain kX D ;1 2 2 (kD ; 1)~! ky ; F (xk )k2 d(L; 2) kx+ ; x1 k2 < 1: k=1 This indicates that kD is a nite number. Q.E.D Theorem 5.2 Assume that F (xk ) ! F (xk ) as ! 0. If k kD for all suciently small, then xk ! xk for k kD as ! 0. Proof. Given su ciently small number , we want to prove kxk ; xk k as ! 0 for k kD . We proceed by induction. Assume that xk ! xk as ! 0, and that k + 1 kD . Note that xk+1 = xk + lk xk+1 = xk + lk lk = F (xk ) lX k ;1 i=1 i zi lk = F (xk ) lX k ;1 i=1 i zi 14 Y.Wang, Y.Yuan we can estimate that kxk+1 ; xk+1 k kxk ; xk k + klk ; lk k lX k ;1 kxk ; xk k + kF (xk ) (i zi ; i zi )k i zi k i=1 kxk ; xk k + kF (xk )k(lk ; 1) max fki zi ; i zi kg i lX k ;1 +kF (xk ) ; F (xk )kk i zi k i=1 +k(F (xk ) ; F (xk ) ) i=1 lX k ;1 (55) By the induction assumption xk ! xk , we have that F (xk ) ! F (xk ): Therefore, it follows that i ! i zi ! zi . Consequently, it from (55) that xk+1 ! xk+1 . Q.E.D Theorem 5.3 Assume that F satises (11) in some ball B D(F ) and let y x as before. Then the iterates xk generated by Algorithms 2.1 and 2.3 converge to a solution of (1) as k ! 1 and ! 0. Proof. For simplicity, We use k() instead of kD in the following analysis. From Theorem 4.5 we know that iterates xk converge to a solution of (1). Combining this fact with Theorem 5.2, we nd that the iterates xk converge to a solution of (1) for k k() as ! 0. Now assume that k() ! 1 as ! 0, and denote x+ the limit of the iterates xk , x+ is a solution of (1). It su eces to consider subsequences fk(n )gn which are monotinically increasing to in nity as n ! 1 and n ! 0. Without loss of generality, let us consider k(m ) > k(n ) for m > n. By the monotonicity of xk , i.e., Theorem 5.1, we have kxkm( m ) ; x+ k kxkm( n ) ; x+ k kxk((m()n) ; xk (n )k + kxk( n ) ; x+ k: Given a su ciently small number > 0 and for some su ciently large number n, kxk( n ) ; x+ k =2 by Theorem 4.5. On the other hand, for su ciently large number m and xed n, kxk((m()n) ; xk (n )k =2 by Theorem 5.1. This proves that kxkm( m ) ; x+ k for all m su ciently large, and thereafter xkm( m ) ! x+ as m ! 1. Hence, we see that xk ! x+ as k ! 1 and ! 0. Q.E.D. Regularity of Trust Region-CG 15 6 Numerical Test In this section, we give an example to test our algorithm. The example is the inverse Gravimetry problem (see 23]). We write it as F (x) : X ;! Y Zb F (x)(t) = k(t s x(s))ds = y(t) t 2 c d]: a (56) with k(t s x(s)) = ln (t;s()t2;+(s)x+(sH);H )2 . Clearly the kernel k is de ned on the set = fc d] a b] Rg and k(t s x(s)) 2 C 1 ( ). The rst derivative F 0 (x) : X ;! Y is de ned by 2 2 Z b @k F (x)u](t) = @x (t s x(s))u(s)ds t 2 c d] a where the kernel @k @x (t s x(s)) can be evaluated by @k (t s x(s)) = 2(H ; x(s)) @x (t ; s)2 + (x(s) ; H )2 : 0 (57) F 0 (x) is compact, since the kernel is square integrable. Now, we will set up the problem of approximate determination of normal pseudosolution to the equation (56). For simplicity, two equidistant grids on intervals a b] and c d] are applied: !n (s) = fsj : sj = a + hs (j ; 1) j = 1 2 ng hs = nb ;; a1 !m (t) = fti : ti = c + ht (i ; 1) i = 1 2 ng ht = md ;; c1 : In this way, the spaces of all grid functions de ned on !n (s) and !m (t), respectively, are treated as Xn and Ym . The integral operator F gives rise to an operator Fmn : Xn ;! Ym by Fmn (x)]i = Zb a k(ti s x(s))ds 1 i m: Similarly the derivative operator F 0 (x) yields an m n matrix: Z b @k Fmn (x)]ij = @x (ti s x(s))j (s)ds 1 i m 1 j n: a 0 Where j (s) we used is the standard linear basic functions 8 s;sj;1 < h if s 2 sj;1 sj ] j (s) = : sj+1h ;s if s 2 sj sj+1 ] 0 else: In which, sj = jh h = n1 j = 1 2 n. The integral can be computed numerically. 16 Y.Wang, Y.Yuan plot of the true and the computed solution plot of the true and the computed solution 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 1: Solution of the inverse problem: nonlinear Fredholm equation We take a b] = c d] = 0 1] H = 0:1 and di"erent m n to give a discretization. Our true function is xtrue (s) = 1:3s(1 ; s) + 0:2, and it is discretized by evaluating it at the points si to give the components xi of x. The right-hand side y is generated by integral (56). The numerical results are shown in gures 1 and 2. In all of these gures, the true solution is denoted by solid line, the approximate solution is denoted by dotted line. First we choose: n = 30, m = 30, !~ = 4:8 = 0:8 0 = 0:1 with perturbation error level = 0:01. The results are shown in the left of Figure 1 It needs 18 inner iterations and 16 outer iterations to generate convergence. Then we choose n = 40, m = 40, !~ = 6:7 = 0:8 0 = 0:1 with small perturbation = 0:005, The results are shown in the right of Figure 1. It needs 20 inner loops and 18 outer loops to generate convergence. Finally we choose n = m = 50, 0 = 0:1 with large perturbation = 0:05 and dominant parameter !~ = 6:8 = 0:9 to give a computation. The results are shown in Figure 2. It needs 7 inner iterations and 10 outer iterations to generate convergence. Remark 6.1 To safeguard Pred is not too small, we choose = 1:0 10;32. If Pred , then we regard it as zero and stop the inner iteration. However, this is not activated in our numerical test. Remark 6.2 In practical applications, the right-hand side is the observation data y which contains noise or error instead of the exact data ytrue . To give a reasonable simulation of the observation data, we add Gaussian white noise rand to the right-hand side ytrue , i.e., ynoise = ytrue + rand (58) where rand is a vector with its components some random numbers in 0 1]. Note that should not be too small or too large. If is too small, according to (58), the noise will not be enough important to give interesting results. However, if is too large, the Regularity of Trust Region-CG 17 plot of the true and the computed solution 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 2: Solution of the inverse problem: nonlinear Fredholm equation observation is a poor approximation to the original problem, it will also not enough to give a reasonable results. 7 Conclusion We have establised the convergence and regularity of the trust region-cg method for nonlinear illposed inverse problems. It deserved pointing out that Plato in 18] had establised the regularity property of the conjugate gradient method. Later on, Hanke in 8] had established the regularity of Newon-CG method. All of the methods are stable for solving ill-posed inverse problems. References 1] Bakushinsky A and Goncharsky A 1994, Ill-Posed Problems: Theory and Applications, Kluwer Academic Publishers. 2] Du" I S, Nocedal J 1987, Reid J K et al., The use of linear programming for the solution of sparse sets of nonlinear equations. SIAM J Sci Stat Comput, 8, 99-108. 3] Engl H W, Hanke M and Neubauer 1996, Regularization of Inverse Problems, Dordrecht: Kluwer. 4] Fletcher R 1987, Practical Methods of Optimization, 2nd ed., Chichester: John Wiley and Sons. 5] Fletcher R 1982, A model algorithm for composite NDO problem, Math Prog Study, 17, 67-76. 18 Y.Wang, Y.Yuan 6] Gay D M 1981, Computing optimal local constrained step, SIAM J Sci Stat Comp, 2, 186-197. 7] Groetsch C W 1984, The Theory of Tikhonov Regularization for Fredholm Equations of The First Kind (Boston,MA:Pitman) 8] Hanke M, Regularizing Properties of a Truncated Newton-CG Algorithm for Nonlinear Inverse Problems, submitted. 9] Hanke M 1997, A Regularizing Levenberg-Marquardt Scheme, with Applications to Inverse Groundwater Filteration Problems, Inverse Problems 13, 79-95. 10] Hanke M 1995, Conjugate Gradient Type Methods for Ill-posed Problems, Longman Scienti c and Technical, Harlow, Essex. 11] Kaltenbacher B 1998, On Broyden's Method for Nonlinear Ill-posed Problems, Numer Funct Anal Opt, 19. 12] Kaltenbacher B 1998, A Posteriori Parameter Choice Strategies for Some Newton Type Methods for the Regularization of Nonlinear Ill-posed Problems, Numer.Math. 79. 13] Kirsch A 1996, An Introduction to the Mathematical Theory of Inverse Problems, New York: Springer-Verlag. 14] Louis A K 1992, Inverse und Schlecht Gestellte Problems, Teubner, Stuttgart. 15] Mor#e J J 1983, Recent developments in algorithms and software for trust region methods. In: Bachem A, Gr$otschel, Korte B, eds. Mathematical Programming: The State of the Art. Berlin: Springer-Verlag, 258-287. 16] Mor#e J J and Sorensen D C 1983, Computing a Trust Region Step, SIAM J.Sci.Stat.Comput. 4, 553-572. 17] Morozov V A 1984 Methods for Solving Incorrectly Posed Problems (New York: Springer) 18] Plato R 1990, Optimal Algorithms for Linear Ill-posed Problems Yield Regularization Methods, Numer.Funct.Anal.and Optim. 11, 111-118. 19] Powell M J D 1984, Nonconvex minimization calculations and the conjugate gradient method. In: Gri ths, ed. Lecture Notes in Mathematics 1066: Numerical Analysis. Berlin: Springer-Verlag, 122-141. 20] Powell M J D 1975, Convergence properties of a class of minimization algorithms. In: Mangasarian O L, Meyer R R, Robinson S M, eds. Nonlinear Programming. New York: Academic Press, 2, 1-27. Regularity of Trust Region-CG 19 21] Sorensen D C 1982, Newton's Method with a Model Trust Region Modi cation, SIAM J.Numer.Anal. 19, 409-426. 22] Steihaug T 1983, The conjugate gradient method and trust regions in large scale optimization. SIAM J Numer Anal 20, 626-637. 23] Tikhonov A N and Arsenin V Y 1977, Solutions of Ill-Posed Problems, New York: Wiley. 24] Toint Ph L 1981, Towards an e cient sparsity exploiting Newton method for minimization. In: Du" I ed. Sparse Matrices and Their Uses. Berlin: Academic Press, 57-88. 25] Yan-fei Wang and Ya-xiang Yuan 2001, A Trust Region Method for Solving Distributed Parameter Identi cation Problems, Report, AMSS, Chinese Academy of Sciences, 2001. 26] Yan-fei Wang, Ya-xiang Yuan and Hong-chao Zhang 2001, A Trust Region-CG Algorithm for Deblurring Problem in Atmospheric Image Reconstruction, to appear in Science in China. 27] Y. Yuan 1994, Nonlinear Programming: Trust Region Algorithms, in: S.T. Xiao and F. Wu, eds., Proceedings of Chinese SIAM annual meeting (Tsinghua University, Beijing) pp. 83 {97. 28] Y. Yuan 1998, Matrix Computation Problems in Trust Region Algorithms for Optimization, in: Q.C. Zeng, T.Q. Li, Z.S. Xue and Q.S. Cheng, eds. Proceedings of the 5th CSIAM annual meeting, (Tsinghua Unviersity Press, Beijing) pp. 52{64. 29] Y. Yuan 1999, Problems on Convergence of Unconstrained Optimization Algorithms, in Y. Yuan ed., Numerical Linear Algebra and Optimization, (Science Press, Beijing, New York), pp. 95-107. 30] Y. Yuan 2000, On the truncated conjugate gradient method. Math Prog, 87, 561-571.