An Optimally Generalized Steepest-Descent Algorithm for Solving Ill-Posed Linear Systems

Abstract

It is known that the steepest-descent method converges normally at the first few iterations, and then it slows down. We modify the original steplength and descent direction by an optimization argument with the new steplength as being a merit function to be maximized. An optimal iterative algorithm with 𝑚-vector descent direction in a Krylov subspace is constructed, of which the 𝑚 optimal weighting parameters are solved in closed-form to accelerate the convergence speed in solving ill-posed linear problems. The optimally generalized steepest-descent algorithm (OGSDA) is proven to be convergent with very fast convergence speed, accurate and robust against noisy disturbance, which is confirmed by numerical tests of some well-known ill-posed linear problems and linear inverse problems.

1. Introduction

The steepest-descent method (SDM), which can be traced back to Cauchy (1847), is the simplest gradient method for solving positive definite linear equations system. The SDM is effective for well-posed and low-dimensional linear problems; however, for large scale linear system and ill-posed linear system it converges very slowly.

Several modifications to the SDM have been addressed. These modifications have led to a new interest in the SDM that the gradient vector itself is not a bad choice but rather that the original steplength used in the SDM leads to the slow convergence behavior. Barzilai and Borwein [1] presented a new choice of steplength through two-point stepsize. Although their method did not guarantee the monotonic descent of the residual norms, Barzilai and Borwein [1] were able to produce a substantial improvement of the convergence speed for a certain test of a positive linear system. The results of Barzilai and Borwein [1] have encouraged many researches on the SDM, for example, Raydan [2, 3], Friedlander et al. [4], Raydan and Svaiter [5], Dai et al. [6], Dai and Liao [7], Dai and Yuan [8], Fletcher [9], and Yuan [10].

The iterative method for solving the system of algebraic equations can be derived from the discretization of certain ordinary differential equations (ODEs) system [11]. In particular, some descent methods can be interpreted as the discretizations of gradient flows [12]. Indeed, the continuous algorithms have been investigated in many literature works for a long time, for example, Gavurin [13], Alber [14], Hirsch and Smale [15], and Chu [16]. The Lyapunov methods used in the analysis of iterative methods have been made by Ortega and Rheinboldt [17] and Bhaya and Kaszkurewicz [11, 18, 19]. Liu [20] has developed a state feedback controlling method together with a Lie-group differential algebraic equations method to solve ill-posed linear systems.

In this paper we solve 𝐁𝐱=𝐛₁, (1) where 𝐱∈ℝ^𝑛 is an unknown vector determined from a given coefficient matrix 𝐁∈ℝ^𝑛×𝑛, which might be unsymmetric, and the input 𝐛₁∈ℝ^𝑛 which might be polluted by random noise. A measure of the ill-posedness of (1) is the condition number of 𝐁: Cond(𝐁)=‖𝐁‖𝐁⁻¹, (2) where ‖𝐁‖ is the Frobenius norm of 𝐁.

Instead of (1), we can solve a normal linear system: 𝐂𝐱=𝐛, (3) where 𝐛=𝐁^T𝐛₁,𝐂=𝐁^T𝐁>𝟎. (4) Throughout this paper the superscript “T” signifies the transpose.

We consider iterative method for solving (3) and define for any vector 𝐱_𝑘 the steepest-descent vector: 𝐫_𝑘∶=𝐂𝐱_𝑘−𝐛. (5) Ascher et al. [21] have viewed the gradient-descent method: 𝐱_𝑘+1=𝐱_𝑘−𝜂_𝑘𝐫_𝑘, (6) as a forward Euler scheme of a time-dependent ODE: 𝐱̇=𝐛−𝐂𝐱. (7) The absolute stability bound 𝜂_𝑘≤2max_𝜎(𝐂)𝜆 (8) must be obeyed if a uniform stepsize is employed, where 𝜎(𝐂) is the set of all the eigenvalues of 𝐂.

Specifically, (6) presents a steepest-descent method (SDM), if the steplength is taken to be 𝜂_𝑘=𝐫_𝑘²𝐫^T_𝑘𝐂𝐫_𝑘. (9) When ‖𝐫_𝑘‖ is rather small the calculated 𝐫_𝑘 may deviate from the real steepest-descent direction to a great extent due to a round-off error of computing machine, which usually leads to the numerical instability of SDM. An improvement of SDM is the conjugate gradient method (CGM), which enhances the search direction of the minimum by imposing an orthogonality to the residual vector at each iterative step.

The relaxed SDM (RSDM) for solving (3) is summarized as follows [22, 23].(i)Give 𝛾 and an initial value of 𝐱₀.(ii)For 𝑘=0,1,…, we repeat the following iterations: 𝐫_𝑘=𝐂𝐱_𝑘−𝐛,𝜂_𝑘=𝐫_𝑘²𝐫^T_𝑘𝐂𝐫_𝑘,𝐱_𝑘+1=𝐱_𝑘−1−𝛾𝜂_𝑘𝐫_𝑘. (10)If 𝐱_𝑘+1 converges according to a given stopping criterion ‖𝐫_𝑘+1‖<𝜀, then stop; otherwise, go to step (ii).

The steepest-descent method is the basis of several gradient-based methods [22, 23], and it is one of the most prominent iterative methods for solving positive definite linear equations system. Although the SDM works very well for most linear systems, the SDM does lose some of its luster for some ill-posed problems like inverse problems, image processing, and box-constrained optimization. Both the efficiency and robustness of the SDM iterative techniques can be improved by using a suitable preconditioning technique.

However, we do not go to the details of the preconditioning techniques of SDM, and instead we will propose an alternative approach by modifying the search direction and steplength used in the SDM by an argument from the optimality of a suitably defined merit function. In this paper we explore a generalization of the SDM by introducing the concept of optimal descent vector to solve linear equations, which is an optimal combination of the steepest-descent vector −𝐫 and the 𝑚-dimensional search vector in a Krylov subspace. Here we modify the direction 𝐫 as well as the steplength 𝜂 from a theoretical optimization being performed on the derived steplength in an 𝑚-dimensional Krylov subspace.

The remaining parts of this paper are arranged as follows. A numerical iterative scheme with −𝐮 as a descent direction is proposed, and the convergence to the minimal point 𝐱^∗=𝐂⁻¹𝐛 of a quadratic functional is proven in Section 2. In Section 3 we propose a generalization of the steepest-descent method to solve positive linear equations systems, by using a new steplength as a merit function to be optimized in an 𝑚-dimensional Krylov subspace. Then, a genuine algorithm is constructed in Section 4, resulting in an optimally generalized steepest-descent algorithm (OGSDA) in terms of 𝑚 optimal weighting parameters which is being derived explicitly to maximize the introduced steplength. The numerical examples are given in Section 5 to display some advantages of the OGSDA, which is compared with CGM, BBM, GMRES, and other algorithms. Finally, the conclusions are drawn in Section 6.

2. The Convergence of New Algorithm

Solving (3) by the steepest-descent method (SDM) is equivalent to solving the following minimum problem: min_{𝐱∈ℝ𝑛}𝜑(𝐱)=min_{𝐱∈ℝ𝑛}12𝐱^T𝐂𝐱−𝐛^T𝐱. (11) In the SDM, we search the next 𝐱(𝑡+Δ𝑡) from the current state 𝐱(𝑡) by minimizing the functional 𝜑 along the direction −𝐫(𝑡); that is, min_𝜂𝜑𝐱(𝑡)−𝜂𝐫(𝑡). (12)

Through some calculations we can obtain ‖‖𝜂=𝐫(𝑡)²𝐫^T(𝑡)(𝑡).𝐂𝐫 (13) Thus we have the following iterative algorithm of SDM: 𝐱()(𝑡)−‖‖𝑡+Δ𝑡=𝐱𝐫(𝑡)²𝐫^T(𝑡)(𝑡)𝐫(𝑡).𝐂𝐫 (14) A straightforward calculation reveals that 𝜑(𝐱(𝑡+Δ𝑡))<𝜑(𝐱(𝑡)). (15) It means that the iterative sequences of SDM converge to the minimal point.

Then, we consider a more general scheme with the descent direction −𝐮(𝑡); that is, 𝐱()(𝑡)−𝐫(𝑡)(𝑡)𝑡+Δ𝑡=𝐱⋅𝐮𝐮^T(𝑡)(𝑡)𝐮(𝑡).𝐂𝐮 (16) Similarly we have the following result.

Lemma 1. Equation (16) is a result of the following minimization along the descent direction −𝐮: min_𝜂𝜑𝐱(𝑡)−𝜂𝐮(𝑡). (17) The functional 𝜑 is decreasing stepwise and the iterative scheme (16) is convergent to the minimal point 𝐱^∗=𝐂⁻¹𝐛.

Proof. Inserting 𝐱(𝑡+Δ𝑡)=𝐱(𝑡)−𝜂𝐮(𝑡) (18) into the functional 𝜑 we have 𝜑()=1𝑡+Δ𝑡2𝐱^T𝐂𝐱−2𝜂𝐮^T𝐂𝐱+𝜂²𝐮^T𝐂𝐮−𝐛^T𝐱+𝜂𝐮^T𝐛, (19) where we simplify the notation of 𝜑(𝐱(𝑡+Δ𝑡)) by 𝜑(𝑡+Δ𝑡) and omit the terms (𝑡) in 𝐱(𝑡) and 𝐮(𝑡). Taking the derivative of (19) with respect to 𝜂 and equating it to zero we can derive 𝐮𝜂=^T𝐫𝐮^T.𝐂𝐮 (20) Because 𝜂 is a steplength, we demand it to be positive; that is, 𝐫⋅𝐮>0. By feeding the above 𝜂 into (18) we can derive the iterative scheme (16).
Inserting (20) into (19) and subtracting the resultant by 𝜑(𝑡)=12𝐱^T𝐂𝐱−𝐛^T𝐱, (21) we can obtain 𝜑()(𝑡)=1𝑡+Δ𝑡−𝜑2𝜂²𝐮^T𝐂𝐮−𝜂𝐮^T𝐮𝐫=−^T𝐫²2𝐮^T.𝐂𝐮 (22) It implies 𝜑()(𝑡)−𝐮𝑡+Δ𝑡=𝜑^T𝐫²2𝐮^T(𝑡),𝐂𝐮<𝜑 (23) and the functional 𝜑 is strictly decreasing. Because the minimal point 𝐱^∗=𝐂⁻¹𝐛 of 𝜑 is just the exact solution of (3), such that when the iterative sequences are generated from scheme (16), the functional 𝜑 decreases step-wisely to its minimal value and 𝐱(𝑡) tends to the minimal point gradually. Thus we can conclude that the new scheme is convergent absolutely.

3. A Generalized SDM and Its Optimization

In Lemma 1 we have proven that the new scheme (16) is convergent, but it cannot tell us what descent direction −𝐮 is the best. Although the SDM uses the steepest-descent direction −𝐫, it is well known that the steplength 𝜂 given in (13) quickly tends to a small value, causing the slowdown of the SDM scheme.

In order to overcome the slowness of SDM, Barzilai and Borwein [1] have presented a new choice of the steplength through two-point stepsize. The algorithm of the Barzilai-Borwein method (BBM) is 𝐱_𝑘+1=𝐱_𝑘−Δ𝐫_𝑘−1^TΔ𝐱_𝑘−1Δ𝐫_𝑘−1²𝐫_𝑘, (24) where Δ𝐫_𝑘−1=𝐫_𝑘−𝐫_𝑘−1 and Δ𝐱_𝑘−1=𝐱_𝑘−𝐱_𝑘−1, with initial guesses 𝐫₀=𝟎 and 𝐱₀=𝟎. Nowadays the improvements of the original BBM have been developed in many literature works [24–30] to treat different ill-posed and inverse problems.

In the numerical solution of linear equations system the Krylov subspace method is one of the most important classes of projective type numerical methods [31–35]. The iterative algorithms that are applied to solve large scale linear systems are mostly the preconditioned Krylov subspace methods [36].

The GMRES used to solve (1) can be summarized as follows [37, 38].(i)Select 𝑚 and give an initial 𝐱₀.(ii)For 𝑘=0,1,…, we repeat the following computations: 𝐫_𝑘=𝐛₁−𝐁𝐱_𝑘,𝐯Arnoldiproceduretosetup_𝑘𝑗,𝑗=1,…,𝑚,from𝐯^𝑘₁=𝐫_𝑘𝐫_𝑘,𝐕_𝑘=𝐯^𝑘₁,…,𝐯^𝑘_𝑚,Solve𝐇^T_𝑘𝐇_𝑘𝐲_𝑘=𝐫_𝑘𝐇^T_𝑘𝐞₁,obtaining𝐲_𝑘,𝐳_𝑘=𝐕_𝑘𝐲_𝑘,𝐱_𝑘+1=𝐱_𝑘+𝐳_𝑘. (25)If 𝐱_𝑘+1 converges according to a given stopping criterion ‖𝐫_𝑘+1‖<𝜀, then stop; otherwise, go to step (ii). In the above, 𝐕_𝑘 is an 𝑛  ×  𝑚 matrix, while 𝐇_𝑘 is an augmented Heissenberg upper triangular matrix with dimension (𝑚  +  1)  ×  𝑚, and 𝐞₁ is the first column of 𝐈_𝑚+1.

3.1. A Generalized SDM

We suppose that the original steepest-descent direction −𝐫 in the SDM is replaced by a new descent vector −𝐮: 𝐮=𝐫+^𝑚_𝑘=1𝛼_k𝐯_𝑘, (26) which is an optimal combination of 𝐫 and the 𝑚-vector 𝐯_𝑘,𝑘=1,…,𝑚. The coefficients 𝛼_𝑘 will be optimized in Section 3.2.

Now we describe how to set up the 𝑚-vector 𝐯_𝑘,𝑘=1,…,𝑚, by the Krylov subspace method. Suppose that we have an 𝑚-dimensional Krylov subspace generated by the coefficient matrix 𝐂 from the steepest-descent vector 𝐫: 𝒦_𝑚∶=span𝐂𝐫,…,𝐂^𝑚𝐫. (27) Then the Arnoldi process is used to normalize and orthogonalize the Krylov vectors 𝐂^𝑗𝐫,𝑗=1,…,𝑚, such that the resultant vectors 𝐯_𝑖, 𝑖=1,…,𝑚, satisfy 𝐯_𝑖⋅𝐯_𝑗=𝛿_𝑖𝑗,𝑖,𝑗=1,…,𝑚, where 𝛿_𝑖𝑗 is the Kronecker delta tensor.

The other way to obtain 𝐯_𝑖,𝑖=1,…,𝑚, is simply by using the unit vectors in the Euclidean space, where 𝐯_𝑖 is the 𝑖th column vector of the identity matrix 𝐈_𝑛. The resultant 𝑚-dimensional space is named the unit subspace.

3.2. Optimization of Algorithm

Let 𝐉 be an 𝑛×𝑚 matrix, which consists of 𝐉∶=𝐯₁,…,𝐯_𝑚. (28) That is, the 𝑗th column of 𝐉 is the vector 𝐯_𝑗. Because 𝐯₁,…,𝐯_𝑚 are linearly independent vectors and 𝑚<𝑛, the expansion matrix 𝐉 has a full rank with rank(𝐉)=𝑚. Then, (26) can be written as 𝐮=𝐫+𝐉𝜶, (29) where 𝜶∶=(𝛼₁,…,𝛼_𝑚)^T.

Now, in view of (16) the new steplength becomes 𝐫𝜂=^T𝐮𝐮^T.𝐂𝐮 (30) When 𝐮=𝐫, the above 𝜂 recovers to that defined in (13). In order to speed up the convergence we can search a way to maximize the steplength, such that we have the following maximization problem: 𝐫max𝜂=^T𝐮𝐮^T.𝐂𝐮 (31) Here, we use the steplength 𝜂 as a merit function to be optimized. Hence, we can prove the following result.

Theorem 2. For rank(𝐉)=𝑚 and 𝐮∈𝐫+𝒦_𝑚, the descent direction −𝐮 which maximizes (31) is given by 𝐮=𝐫−𝐄𝐂𝐫+𝜆𝐄𝐫,𝐀=𝐉^T𝐂𝐉,𝐄=𝐉𝐀⁻¹𝐉^T,𝑏₀=𝐫^T()𝑏𝐂𝐄𝐂−𝐂𝐫<0,₁=𝐫^T𝐈_𝑛𝑏−𝐂𝐄𝐫>0,₂=𝐫^T𝐄𝐫>0,𝜆=𝑏²₁−𝑏₀𝑏₂−𝑏₁𝑏₂. (32)(33)(34)(35)(36) Moreover, the steplength 𝜂 defined in (30) is positive; that is, 1𝜂==𝑏2𝜆²2𝑏²₁−𝑏₀𝑏₂−2𝑏₁>0. (37)

Proof. We divide this proof into three parts.(A)A Minimization Problem. Inserting (29) into (31) we can obtain the following minimization problem: min_{𝛼1,…,𝛼𝑚}𝐮⋅(𝐂𝐮),𝐫⋅𝐮 (38)where 𝐫⋅𝐮=‖𝐫‖²+𝐫^T()𝐉𝜶,𝐮⋅𝐂𝐮=𝐫^T𝐂𝐫+2𝐫^T𝐂𝐉𝜶+𝜶^T𝐀𝜶, (39)(40) in which 𝐀∶=𝐉^T𝐂𝐉 (41) is an 𝑚×𝑚 positive definite matrix. Because the matrix 𝐂 has a full rank with rank(𝐂)=𝑛 and positive definite and 𝐉 has a full rank with rank(𝐉)=𝑚, from the above equation it follows that 𝐀 has a full rank with rank(𝐀)=𝑚. Then, we confirm that 𝐀 is an 𝑚×𝑚 positive definite matrix by (41).
From (39) and (40) we can derive 𝐲₁∶=∇_𝜶𝐫⋅𝐮=𝐉^T𝐲𝐫,₂∶=∇_𝜶𝐮⋅(𝐂𝐮)=2𝐉^T𝐂𝐫+2𝐀𝜶, (42)(43) where ∇_𝜶 denotes the gradient with respect to 𝜶.
By using ∇_𝜶𝐮⋅(𝐂𝐮)𝐫⋅𝐮=𝟎⟹𝐫⋅𝐮∇_𝜶𝐮⋅(𝐂𝐮)−𝐮⋅(𝐂𝐮)∇_𝜶𝐫⋅𝐮=𝟎, (44) we can derive the following governing equation of vector form to solve the optimal parameter 𝜶: 𝐫⋅𝐮𝐲₂−𝐮⋅(𝐂𝐮)𝐲₁=𝟎. (45)
Now we prove that 𝛼 solved from the above equation is a minimal point of the functional in (38). We can compute the Jacobian matrix of (45), which is denoted by 𝐉_𝜶: 𝐉_𝜶=2𝐫⋅𝐮𝐀+𝐲₂𝐲^T₁−𝐲₁𝐲^T₂. (46) It is an 𝑚×𝑚 positive definite matrix, because for any 𝐳≠𝟎∈ℝ^𝑚 we have 𝐳^T𝐉_𝜶𝐳=2𝐫⋅𝐮𝐳^T𝐀𝐳>0, (47) where we have used the positive definite property of 𝐀, and 𝐫⋅𝐮>0. Because 𝜂 defined by (20) is a steplength, we demand it to be positive; that is, 𝐫⋅𝐮>0 (see the proof given below for 𝐫⋅𝐮>0).(B)A Quadratic Equation. Equations (40) and (43) can be written as 𝐮⋅(𝐂𝐮)=𝐫^T𝐂𝐫+2𝐲^T₀𝐉𝜶+𝜶^T𝐲𝐀𝜶,₂=2𝐉^T𝐲₀+2𝐀𝜶, (48)(49)where 𝐲₀∶=𝐂𝐫. (50)
From (45) we can observe that 𝐲₂ is proportional to 𝐲₁, which is supposed to be 𝐲₂=𝐮⋅(𝐂𝐮)𝐲𝐫⋅𝐮¹=2𝜆𝐲₁, (51) where 2𝜆 is a multiplier to be determined. From (51), two equations follow 𝐲₂=2𝜆𝐲₁,𝐮⋅(𝐂𝐮)=2𝜆𝐫⋅𝐮. (52)(53) Then, by (42), (49), and (52) we can solve 𝜶 by 𝜶=𝜆𝐀⁻¹𝐉^T𝐫−𝐀⁻¹𝐉^T𝐲₀. (54)
Inserting (54) for 𝜶 into (39) and (48) we have that 𝐫⋅𝐮=‖𝐫‖²−𝐲^T₀𝐄𝐫+𝜆𝐫^T()𝐄𝐫,𝐮⋅𝐂𝐮=𝐫^T𝐂𝐫−𝐲^T₀𝐄𝐲₀+𝜆²𝐫^T𝐄𝐫, (55) where 𝐄∶=𝐉𝐀⁻¹𝐉^T (56) is an 𝑛×𝑛 positive semi-definite matrix.
Now, from (53) and (55) we can derive a quadratic equation to solve 𝜆: 𝑎₂𝜆²+𝑎₁𝜆+𝑎₀=0, (57) where 𝑎₀∶=𝐲^T₀𝐄𝐲₀−𝐫^T𝑎𝐂𝐫,₁∶=2‖𝐫‖²−2𝐲^T₀𝐄𝐫,𝑎₂∶=𝐫^T𝐄𝐫. (58)(C)A Closed-Form Solution of Optimal Parameter 𝜶. Now we prove 𝑎₀<0. By using (50) we have 𝑎₀=𝐫^T[𝐂𝐄𝐂−𝐂]𝐫<0, (59)if we can prove 𝐂𝐄𝐂−𝐂<𝟎. (60) From (56) and (41) it follows that 𝐄𝐂𝐉=𝐉𝐉^T𝐂𝐉⁻¹𝐉^T𝐉𝐂𝐉=𝐉,_T𝐂𝐄=𝐉^T𝐂𝐉𝐉^T𝐂𝐉⁻¹𝐉^T=𝐉^T,𝐉^T[𝐂𝐄𝐂]𝐉=𝐉^T𝐂𝐉𝐉^T𝐂𝐉⁻¹𝐉^T𝐂𝐉=𝐉^T𝐂𝐉. (61)(62)(63) They indicate that 𝐄 is a Penrose pseudoinverse of 𝐂 in the Krylov subspace; hence, (60) follows.
According to (58), it is obvious that 𝑎₂>0; hence, by (59) we have 𝑎²₁−4𝑎₀𝑎₂>0. (64) As a consequence, (57) has a closed-form positive real solution: 𝜆=𝑎²₁−4𝑎₀𝑎₂−𝑎₁2𝑎₂. (65) Inserting it into (54) we can obtain the closed-form solution of 𝜶: 𝜶=𝑎²₁−4𝑎₀𝑎₂−𝑎₁2𝑎₂𝐀⁻¹𝐉^T𝐫−𝐀⁻¹𝐉^T𝐲₀. (66)
Upon comparing (53) and (30) we can derive (37) with the steplength 𝜂 being positive. Moreover, as a direct result of (30) and (37) we have 𝐫⋅𝐮>0.

Remark 3. Because of 𝑎₂>0 and 𝑎₀<0 in (57), it has a positive solution as shown in (65) and in addition a negative solution as follows: −𝜆=𝑎²₁−4𝑎₀𝑎₂−𝑎₁2𝑎₂, (67) where 𝑎₁=2𝐫^T𝐈_𝑛−𝐂𝐄𝐫>0, (68) in view of (58), (50), and (62). Hence, by (30) and (37) we have 𝐫⋅𝐮=‖𝐫‖²−𝐲^T₀𝐄𝐫+𝜆𝐫^T𝐄𝐫<0, (69) which contradicts (20), where 𝜂 is a steplength being positive. Below we will give a numerical example to demonstrate this fact.

Corollary 4. If 𝑚=𝑛, then one has 𝜶=𝐉⁻¹𝐂⁻¹𝐫−𝐉⁻¹𝐫. (70)

Proof. If 𝑚=𝑛, then 𝐉 defined by (28) is an 𝑛×𝑛 nonsingular matrix. Simultaneously, 𝐄 defined by (56) is equal to 𝐂⁻¹, 𝑎₀=0, 𝑎₁=0, 𝑎₂>0, and meanwhile (66) reduces to (70), where we have taken 𝜆=1 according to the above condition.

Remark 5. If 𝑚=𝑛, by (29) and (70) we have 𝐮=𝐫+𝐉𝜶=𝐫+𝐉𝐉⁻¹𝐂⁻¹𝐫−𝐉⁻¹𝐫=𝐂⁻¹𝐫, (71) and inserting it into (31) we can prove 𝜂=1. In any way, under this condition we have the fastest convergent behavior, because the solution 𝐱=𝐂⁻¹𝐛 is already obtained. Usually, we cannot take 𝑚=𝑛, which leads to a quite difficult task to find the inverse of 𝐀=𝐉^T𝐂𝐉 when 𝑚=𝑛 is a large number.

4. An Optimal Generalized Steepest-Descent Algorithm

Thus, we arrive at the following optimally generalized steepest-descent algorithm (OGSDA).(i)Select 𝑚 and 0≤𝛾<1, and give an initial value of 𝐱₀.(ii)For 𝑘=0,1,…, we repeat the following computations: 𝐫_𝑘=𝐂𝐱_𝑘𝐯−𝐛,Arnoldiproceduretosetup_𝑘𝑗,𝑗=1,…,𝑚,from𝐯^𝑘₁=𝐂𝐫_𝑘𝐂𝐫_𝑘,𝐉_𝑘=𝐯^𝑘₁,…,𝐯^𝑘_𝑚,𝐀_𝑘=𝐉^T_𝑘𝐂𝐉_𝑘,𝐄_𝑘=𝐉_𝑘𝐀_𝑘⁻¹𝐉^T_𝑘,Compute(36),obtaining𝜆_𝑘,𝐮^𝑘=𝐫_𝑘−𝐄_𝑘𝐂𝐫_𝑘+𝜆_𝑘𝐄_𝑘𝐫_𝑘,𝜂_𝑘=12𝜆_𝑘,𝐱_𝑘+1=𝐱_𝑘−1−𝛾𝜂_𝑘𝐮^𝑘. (72)If 𝐱_𝑘+1 converges according to a given stopping criterion ‖𝐫_𝑘+1‖<𝜀, then stop; otherwise, go to step (ii).

Remark 6. Depending on the different choice of 𝐯_𝑗,𝑗=1,…,𝑚, the optimal algorithm is further classified into two types. When 𝐯_𝑗,𝑗=1,…,𝑚, are constructed from the Krylov subspace and followed by the Arnoldi process to orthonormalize the Krylov vectors, the resultant optimal algorithm is named the OGSDA with Krylov subspace method. On the other hand, if 𝐯_𝑗 is simply the 𝑗th column of the unit matrix 𝐈_𝑛, then the resultant optimal algorithm is named the OGSDA with unit subspace method. Both methods have their own advantage. In the Krylov subspace method, the expansion matrix 𝐉 is not fixed, and it can adjust its configuration at each iterative step to accelerate the convergence, which however leads to a consumption of computational time in the calculation of the inverse of 𝐀=𝐉^T𝐂𝐉 to obtain 𝐀⁻¹ and 𝐄. On the other hand, in the unit subspace method, the expansion matrix 𝐉 is fixed, and we only need to calculate 𝐀=𝐉^T𝐂𝐉, 𝐀⁻¹ and 𝐄 one time, which is much time saving than that by using the Krylov subspace method; however, it cannot adjust the configuration of the descent direction at each iterative step, which might lead to an ill-condition in the course of computation.
We can evaluate the cost of computing the approximate solution 𝐱_𝑘 by OGSDA with 𝑘 steps. We can estimate the cost of the computation of the scalar 𝜆, which at each iterative step requires an inversion of 𝑚×𝑚 matrix, the matrix-vector multiplications of 𝑛×𝑛-matrix and 𝑛-vector three times, and the inner products of two 𝑛-vectors five times. This portion is very time saving with totally 3𝑘𝑛+𝑘𝑚², because of 𝑚≪𝑛. The computation of the Arnoldi vectors needs totally 𝑘(𝑘+1)𝑛+𝑘𝑛² multiplications. In the computation of 𝐮 it requires 2𝑘𝑛 multiplications. The 𝑘 steps OGSDA therefore requires totally 𝑘(𝑘+6)𝑛+𝑘𝑛²+𝑘𝑚² multiplications. Dividing by the total number of steps 𝑘 we can obtain that each step requires (𝑘+6)𝑛+𝑛²+𝑚² multiplications on the average. The computational cost of OGSDA is slightly larger than that of GMRES and BBM per step. However, because the convergence speed of the OGSDA is very fast, the total computational time is less than that of GMRES and BBM.

5. Numerical Examples

In order to evaluate the performance of the OGSDA, we test some well-known ill-posed linear problems of the Hilbert problem, the backward heat conduction problem, the heat source identification problem, and the inverse Cauchy problem.

5.1. Example 1

First by testing the performance of OGSDA on the solution of linear equations system, we consider the following convex quadratic programming problem with equality constraint: 1min𝑓=2𝐱^T𝐏𝐱+𝐪^T𝐱,𝐐𝐱=𝐛₀, (73)(74) where 𝐏 is an 𝑛₁×𝑛₁ matrix, 𝐐 is an 𝑚₁×𝑛₁ matrix, and 𝐛₀ is an 𝑚₁-vector, which means that (74) provides 𝑚₁ linear constraints. According to the Lagrange theory we need to solve (1) with the following 𝐛 and 𝐁: 𝐛₁𝐛∶=−𝐪₀,𝐁=𝐏𝐐^T𝐐𝟎_𝑚1×𝑚1. (75)

For definite we take 𝑛₁=3 and 𝑚₁=2 with min𝑓=𝑥²₁+2𝑥²₂+𝑥²₃−2𝑥₁𝑥₂+𝑥₃,𝑥₁+𝑥₂+𝑥₃=4,2𝑥₁−𝑥₂+𝑥₃=2. (76) Under the following parameters 𝑚=2, 𝛾=0.2, and 𝜀=10⁻⁵ and starting from the initial values of 𝑥₁=𝑥₂=𝑥₃=1 and two Lagrange multipliers 𝜆₁=−2 and 𝜆₂=2 we apply the OGSDA to solve the resultant five-dimensional linear equations system, which is convergent with 38 steps as shown in Figure 1(a) for the residual. In Figures 1(b) and 1(c) we show, respectively, the optimal parameters 𝛼₁ and 𝛼₂ and the values of 𝜂_𝑘 and 𝑎^𝑘₀. It can be seen that 𝑎^𝑘₀∈[−22.16,−5.72×10⁻¹³], as proven by (59) 𝑎^𝑘₀<0. Now we raise 𝑚 to 𝑚=𝑛=5 and also the convergence criterion to 𝜀=10⁻⁵ and 𝛾=0 is fixed. The numerical results obtained by the OGSDA with Krylov subspace is convergent with only three steps. The final solution is a minimal value min𝑓=3.977273 at the point (𝑥₁,𝑥₂,𝑥₃)=(1.909090909092,1.954545454479,0.13636363636).

(a)

(b)

(c)

(d)

Figure 1 

Example 1 showing (a) residual, ((b) and (c)) optimal parameters, and (d) steplength and 𝑎0.

In order to claim that the solution of 𝜆 given by (67) is incorrect, we use the corresponding algorithm to solve the above problem under the parameters 𝑚=2 and 𝛾=0.2. First the algorithm does not converge within 500 steps under a weak convergence criterion with 𝜀=10⁻² as shown in Figure 2(a). Then the algorithm gives negative steplength of 𝜂 as shown in Figure 2(b).

(a)

(b)

Figure 2 

For Example 1, using (67) the algorithm does not converge and gives a negative steplength.

5.2. Example 2

Finding an 𝑛-order polynomial function 𝑝(𝑥)=𝑎₀+𝑎₁𝑥+⋯+𝑎_𝑛𝑥^𝑛 to best match a continuous function 𝑓(𝑥) in the interval of 𝑥∈[0,1]min_{deg(𝑝)≤𝑛}¹₀𝑓(𝑥)−𝑝(𝑥)²𝑑𝑥, (77) leads to a problem governed by (3), where 𝐂 is the (𝑛+1)×(𝑛+1) Hilbert matrix defined by 𝐶_𝑖𝑗=1,𝑖+𝑗−1 (78)𝐱 is composed of the 𝑛+1 coefficients 𝑎₀,𝑎₁,…,𝑎_𝑛 that appeared in 𝑝(𝑥), and ∫𝐛=¹⁰𝑓(𝑥)∫𝑑𝑥₁₀(𝑥)⋮∫𝑥𝑓𝑑𝑥₁₀𝑥^𝑛𝑓(𝑥)𝑑𝑥 (79) is uniquely determined by the function 𝑓(𝑥).

The Hilbert matrix is a notorious example of highly ill-conditioned matrices. Equation (3) with the matrix 𝐂 having a large condition number usually displays that an arbitrarily small perturbation of data on the right-hand side may lead to an arbitrarily large perturbation to the solution on the left-hand side.

In this example we consider a highly ill-conditioned linear system (3) with 𝐂 given by (78). The ill-posedness of (3) increases fast with 𝑛. We consider an exact solution with 𝑥_𝑗=1 and 𝑏_𝑖 is given by 𝑏_𝑖=^𝑛_𝑗=11(𝑖),𝑖+𝑗−1+𝜎𝑅 (80) with 𝑅(𝑖) being random numbers between [−1,1].

Liu [22] has applied the relaxed steepest-descent method (RSDM) with 𝛾=0.06, starting from 𝑥₁=⋯=𝑥₉=0.5, to solve the Hilbert problem with 𝑛=9 with a stopping criterion 𝜀=10⁻⁸ for the relative residual. Through 50000 iterations the numerical solution converges to the exact solution very accurately as shown in Table 1 with the maximum error being 1.44×10⁻³. When we apply the OGSDA with the unit subspace method under 𝑚=𝑛=9 and 𝛾=10⁻⁵ to solve this problem it converges only three steps and with the maximum error being 4.68×10⁻⁶. It is very time saving because we only need to calculate 𝐀⁻¹ one time, by which using the matrix conjugate gradient method (MCGM) developed by Liu et al. [39] we only spend 103 steps under the convergence criterion 𝜀₁=10⁻⁵ to find 𝐀⁻¹. For the purpose of comparison, the values obtained by other numerical methods are also listed in Table 1.

Table 1 

Comparing the numerical results for Example  2 with different methods.

As shown in Corollary 4 when 𝑚 is taken to be 𝑚=𝑛 we have the fastest convergent solution. Usually, we can take a smaller 𝑚 for the easy solution of 𝐀⁻¹. For example, if we take 𝑚=5 in the OGSDA with the Krylov subspace method, we can obatin the solution with the maximum error being 4.45×10⁻⁴ through only four steps. Under the same conditions, the CGM is convergent with six steps with the maximum error being 7.857×10⁻³.

Next we consider a more difficult case with 𝑛=300. In the computation a random noise in the level of 𝜎=10⁻⁶ is imposed on the right-hand side, and the convergence criteria used in CGM and OGSDA are both 𝜀=10⁻² for the relative residual. The CGM converges very fast with 4 iterations as shown in Figure 3(a); however, the maximum error of CGM is large up to 0.25. For this noised case the OGSDA is applicable and convergent with 4 iterations, where 𝛾=0.15 and 𝑚=10 are used. The maximum error of OGSDA is 0.0113. In Figure 3 we compare (a) relative residuals, (b) steplength of OGSDA, and (c) numerical errors of OGSDA, CGM, GMRES, and BBM. It is obvious that for this highly ill-posed problem the OGSDA is convergent fast, and accurate. The solutions obtained by the CGM, GMRES, and BBM are not good, although the GMRES with 𝑚=9 converges with only two steps, CGM with four steps, and the BBM with ten steps.

(a)

(b)

(c)

Figure 3 

Example 2 of the Hilbert problem with 𝑛=300, showing (a) relative residuals, (b) steplength of OGSDA, and (c) comparing the numerical errors of CGM, GMRES, OGSDA, and BBM.

5.3. Example 3

When we apply a central difference scheme to the following two-point boundary value problem: −𝑢^′′(𝑥)(𝑥)𝑢(0)(1)=𝑓,0<𝑥<1,=𝑎,𝑢=𝑏, (81) we can derive a linear equations system 𝑢𝐂𝐮=2−1−12−1⋅⋅⋅⋅⋅⋅⋅⋅⋅−12¹𝑢₂⋮𝑢_𝑛=(Δ𝑥)²𝑓()()Δ𝑥+𝑎Δ𝑥²𝑓()⋮()2Δ𝑥Δ𝑥²𝑓))()((𝑛−1Δ𝑥Δ𝑥²𝑓(),𝑛Δ𝑥+𝑏 (82) where Δ𝑥=1/(𝑛+1) is the spatial length and 𝑢_𝑖=𝑢(𝑖Δ𝑥),𝑖=1,…,𝑛, are unknown values of 𝑢(𝑥) at the grid points 𝑥_𝑖=𝑖Δ𝑥. 𝑢₀=𝑎 and 𝑢_𝑛+1=𝑏 are the given boundary conditions. The above matrix 𝐂 is known as a central difference matrix.

The eigenvalues of 𝐂 are found to be 4sin²𝑘𝜋2(𝑛+1),𝑘=1,2,…,𝑛, (83) which together with the symmetry of 𝐂 indicates that 𝐂 is positive definite, and Cond(𝐂)=sin²(𝑛𝜋/2(𝑛+1))sin²(𝜋/2(𝑛+1)) (84) is a large number when the grid number 𝑛 is very large.

In this numerical test we fix 𝑛=200 and thus the condition number of 𝐂 is 16210.7. Let us consider the boundary value problem in (81) with 𝑓(𝑥)=sin𝜋𝑥. The exact solution is 𝑢(𝑥)()1=𝑎+𝑏−𝑎𝑥+𝜋²sin𝜋𝑥, (85) where we fix 𝑎=1 and 𝑏=2.

A relative random noise with intensity 𝜎=0.01 is added by 𝑓(𝑖Δ𝑥)=𝑓(𝑖Δ𝑥)[1+𝜎𝑅(𝑖)] into the right-hand side data of (82). Under a moderate convergence criterion with 𝜀=10⁻⁷ for the relative residual, we compare (a) relative residuals and (c) numerical errors of OGSDA, CGM, GMRES, and BBM in Figure 4. The number of iterations of CGM is 1325, and the maximum error is 3.14×10⁻⁵. Then we apply the OGSDA with 𝛾=0.25 and 𝑚=30 to this problem, whose number of iterations is 66, and the maximum error is 1.9×10⁻⁵. The number of iterations of GMRES with 𝑚=10 is 811, and the maximum error is 4.21×10⁻⁵, while the number of iterations of BBM is 1708, and the maximum error is 1.06×10⁻⁵.

(a)

(b)

(c)

Figure 4 

Example 3, comparing (a) relative residuals, (b) convergence rates, and (c) numerical errors of CGM, OGSDA, GMRES, and BBM.

It can be seen that the OGSDA is more convergent, faster, and more accurate than CGM and GMRES and much faster than the BBM; however, the accuracy of BBM is slightly better than that obtained by the OGSDA.

In order to demonstrate the efficiency of the present OGSDA with that of the CGM we can compare the convergence rate of these two algorithms. According to Liu [40, 41] the convergence rate is given by ConvergenceRate_OGSDA1:=√𝑠_𝑘,𝑠_𝑘=1−1−𝛾²𝐂𝐱_𝑘−𝐛⋅𝐂𝐮^𝑘2𝐂𝐱_𝑘−𝐛²𝐂𝐮^𝑘2,ConvergenceRate_CGM1:=√𝑠_𝑘,𝑠_𝑘𝜂=1−^𝑘𝐁𝐱_𝑘−𝐛₁², (86) where 𝜂_𝑘=‖𝐫_𝑘‖²/𝐩^T_𝑘𝐂𝐩_𝑘, in which 𝐩_𝑘 is the supplemental vector used in the CGM. In Figure 4(b) we compare the convergence rates of CGM and OGSDA, of which we can see that the convergence rate of OGSDA is much larger than that of CGM.

5.4. Example 4

When the backward heat conduction problem (BHCP) is considered in a spatial interval of 0<𝑥<ℓ by subjecting to the boundary conditions at two ends of a slab: 𝑢_𝑡(𝑥,𝑡)=𝜅𝑢_𝑥𝑥()𝑢()𝑥,𝑡,0<𝑡<𝑇,0<𝑥<ℓ,0,𝑡=𝑢₀(𝑡),𝑢(ℓ,𝑡)=𝑢_ℓ(𝑡), (87)(88) we solve 𝑢 under a final time condition: 𝑢(𝑥,𝑇)=𝑢^𝑇(𝑥). (89)

The fundamental solution of (87) is given by 𝐾()=𝐻(𝑡)𝑥,𝑡2√𝜅𝜋𝑡exp−𝑥²,4𝜅𝑡 (90) where 𝐻(𝑡) is the Heaviside function.

The method of fundamental solutions (MFS) has a serious drawback that the resulting linear equations system is always highly ill-conditioned, when the number of source points is increased, or when the distances of source points are increased.

In the MFS the solution of 𝑢 at the field point 𝐳=(𝑥,𝑡) can be expressed as a linear combination of the fundamental solutions 𝑈(𝐳,𝐬_𝑗): 𝑢(𝐳)=^𝑛_𝑗=1𝑐_𝑗𝑈𝐳,𝐬_𝑗,𝐬_𝑗=𝜂_𝑗,𝜏_𝑗∈Ω^𝑐, (91) where 𝑛 is the number of source points, 𝑐_𝑗 are unknown coefficients, and 𝐬_𝑗 are source points being located in the complement Ω^𝑐 of Ω=[0,ℓ]×[0,𝑇]. For the heat conduction equation we have the basis functions: 𝑈𝐳,𝐬_𝑗=𝐾𝑥−𝜂_𝑗,𝑡−𝜏_𝑗. (92)

It is known that the location of source points in the MFS has a great influence on the accuracy and stability. In a practical application of MFS to solve the BHCP, the source points are uniformly located on two vertical straight lines parallel to the 𝑡-axis, not over the final time, which was adopted by Hon and Li [42] and Liu [43], showing a large improvement than the line location of source points below the initial time. After imposing the boundary conditions and the final time condition to (91) we can obtain a linear equations system: 𝐁𝐱=𝐛₁, (93) where 𝐵_𝑖𝑗=𝑈𝐳_𝑖,𝐬_𝑗,𝐱=𝑐₁,…,𝑐_𝑛^T,𝐛₁=𝑢_ℓ𝑡_𝑖,𝑖=1,…,𝑚₁;𝑢^𝑇𝑥_𝑗,𝑗=1,…,𝑚₂;𝑢₀(𝑡_𝑘),𝑘=𝑚₁,…,1^T, (94) and 𝑛=2𝑚₁+𝑚₂.

Since the BHCP is highly ill-posed, the ill-condition of the coefficient matrix 𝐁 in (93) is serious. To overcome the ill-posedness of (93) we can use the OGSDA to solve this problem. Here we compare the numerical solution with an exact solution: 𝑢(𝑥,𝑡)=cos(𝜋𝑥)exp−𝜋²𝑡. (95) For the case with 𝑇=1 the value of final time data is in the order of 10⁻⁴, which is small by comparing with the value of the initial temperature 𝑓(𝑥)=𝑢₀(𝑥)=cos(𝜋𝑥) to be retrieved, which is 𝑂(1). First we impose a relative random noise with an intensity 𝜎=10% being imposed on the final time data. Under the following parameters 𝑚₁=11, 𝑚₂=6, 𝑚=10, 𝛾=0.25, and 𝜀=𝜀₁=10⁻⁵, we solve this problem by the OGSDA with the Krylov subspace method. With 28 iterations the OGSDA can obtain a very accurate solution with the maximum error being 9.877×10⁻⁴, which is better than 1.247×10⁻³ calculated by Liu [35] using the optimal multivector iterative algorithm (OMVIA). In Figures 5(a) and 5(b) we plot the relative residual and steplength, of which we can see that the minimal steplength is 1. The numerical result is compared with the exact initial value in Figure 5(c), whose numerical error as shown in Figure 5(d) indicates that the present OGSDA is very robust against noise. Under the same parameters we also apply the BBM and GMRES with 𝑚=10 to solve the BHCP, of which both algorithms cannot converge to satisfy the convergence criterion with 𝜀=10⁻⁵, and we let BBM run 100 steps and the GMRES run 30 steps with the residuals being shown in Figure 5(a). The maximum error obtained by the BBM is 0.151, while that for the GMRES is 0.809 as shown in Figure 5(d). For this problem we can quickly provide a very accurate numerical result by using the OGSDA, which is much better than the other three algorithms.

(a)

(b)

(c)

(d)

Figure 5 

Example 4 showing (a) relative residual and (b) steplength of OGSDA, (c) comparing numerical and exact initial value, and (d) comparing numerical errors of OMVIA, OGSDA, BBM, and GMRES.

5.5. Example 5

Let us consider the inverse Cauchy problem for the Laplace equation: Δ𝑢=𝑢_𝑟𝑟+1𝑟𝑢_𝑟+1𝑟²𝑢_𝜃𝜃𝑢(𝜃)𝑢=0,𝜌,𝜃=ℎ,0≤𝜃≤𝜋,_𝑛𝜌,𝜃=𝑔(𝜃),0≤𝜃≤𝜋, (96)(97)(98) where ℎ(𝜃) and 𝑔(𝜃) are given function. The inverse Cauchy problem is specified as follows.

To seek an unknown boundary function 𝑓(𝜃) on the part Γ₂∶={(𝑟,𝜃)∣𝑟=𝜌(𝜃),𝜋<𝜃<2𝜋} of the boundary under (96)–(98) with the overspecified data being given on Γ₁∶={(𝑟,𝜃)∣𝑟=𝜌(𝜃),0≤𝜃≤𝜋}.

It is well known that the method of fundamental solutions (MFS) can be used to solve the Laplace equation when a fundamental solution is known. In the MFS the solution of 𝑢 at the field point 𝐱=(𝑟cos𝜃,𝑟sin𝜃) can be expressed as a linear combination of fundamental solutions 𝑈(𝐱,𝐬_𝑗): 𝑢(𝐱)=^𝑛_𝑗=1𝑐_𝑗𝑈𝐱,𝐬_𝑗,𝐬_𝑗∈Ω^𝑐. (99) For the Laplace equation (96) we have the fundamental solutions: 𝑈𝐱,𝐬_𝑗=ln𝑟_𝑗,𝑟_𝑗=𝐱−𝐬_𝑗. (100)

In the practical application of MFS, by imposing the boundary conditions (97) and (98) at 𝑁 points on (99) we can obtain a linear equations system: 𝐁𝐜=𝐛₁, (101) where 𝐱_𝑖=𝑥¹_𝑖,𝑥²_𝑖=𝜌𝜃_𝑖cos𝜃_𝑖,𝜌𝜃_𝑖sin𝜃_𝑖,𝐬_𝑗=𝑠¹_𝑗,𝑠²_𝑗=𝑅𝜃_𝑗cos𝜃_𝑗,𝑅𝜃_𝑗sin𝜃_𝑗,𝐵_𝑖𝑗=ln𝐱_𝑖−𝐬_𝑗𝐵,if𝑖isodd,_𝑖𝑗=𝜂𝜃_𝑖𝐱_𝑖−𝐬_𝑗²𝜌𝜃_𝑖−𝑠¹_𝑗cos𝜃_𝑖−𝑠²_𝑗sin𝜃_𝑖−𝜌^′𝜃_𝑖𝜌𝜃_𝑖×𝑠¹_𝑗sin𝜃_𝑖−𝑠²_𝑗cos𝜃_𝑖𝑐,if𝑖iseven,𝐜=₁,…,𝑐_𝑛^T,𝐛₁=ℎ𝜃₁,𝑔𝜃₁,…,ℎ𝜃_𝑁,𝑔𝜃_𝑁^T, (102) in which 𝑛=2𝑁, and 𝜂(𝜃)=𝜌(𝜃)𝜌²(𝜃)+𝜌^′(𝜃)². (103) The above 𝑅(𝜃)=𝜌(𝜃)+𝐷 with an offset 𝐷 can be used to locate the source points along a contour with a radius 𝑅(𝜃). When the linear equations system (101) is established, we can apply the OGSDA to solve it.

For the purpose of comparison we consider the following exact solution: 𝑢𝑥,𝑦=cos𝑥cosh𝑦+sin𝑥sinh𝑦, (104) defined in a domain with a complex amoeba-like irregular shape as a boundary: 𝜌(𝜃)=exp(sin𝜃)sin²(2𝜃)+exp(cos𝜃)cos²(2𝜃). (105) After imposing the boundary conditions (97) and (98) at 𝑁 points on (99) we can obtain a linear equations system. The noise being imposed on the measured data ℎ and 𝑔 is 𝜎=0.01.

We solve this problem by the OGSDA with 𝑚=10 and 𝛾=0.2. Through 15 iterations the residual is smaller than the convergence criterion 𝜀=10⁻⁵ as shown in Figure 6(a), and also the steplength is shown there. We compare the recovered data computed by the OGSDA with the exact one in Figure 6(b). The numerical error as shown in Figure 6(c) is smaller than 0.062. It can be seen that the OGSDA with Krylov subspace can accurately recover the unknown boundary condition.

(a)

(b)

(c)

Figure 6 

Example 5 showing (a) residual and steplength of OGSDA, (b) comparing numerical and exact solutions, and (c) numerical errors of OGSDA, GMRES, and BBM.

We also apply the GMRES with 𝑚=5 and the BBM to this inverse Cauchy problem, whose residuals are shown in Figure 7. They fast tend to a plateau and then the residuals are no more reduced. The numerical solutions are very bad with the maximum error of BBM being 3.79 and the maximum error of GMRES being 2.81. As compared with those obtained by OGSDA, the accuracy and convergence speed of OGSDA are much better than those obtained by the GMRES and the BBM.

Figure 7 

Example 5 showing residuals of GMRES and BBM.

5.6. Example 6

In this section we apply the OGSDA to identify an unknown space-dependent heat source function 𝐻(𝑥) for a one-dimensional heat conduction equation: 𝑢_𝑡(𝑥,𝑡)=𝑢_𝑥𝑥(𝑥,𝑡)+𝐻(𝑥),0<𝑥<ℓ,0<𝑡<𝑡_𝑓,𝑢(0,𝑡)=𝑢₀(𝑡),𝑢(ℓ,𝑡)=𝑢_ℓ(𝑡),𝑢(𝑥,0)=𝑓(𝑥). (106)(107)(108) In order to identify 𝐻(𝑥) we impose a Neumann type boundary condition: 𝑢_𝑥(0,𝑡)=𝑞(𝑡). (109)

We propose a numerical differential method by letting 𝑣=𝑢_𝑡. Taking the differentials of (106), (107), and (109) with respect to 𝑡 and letting 𝑣=𝑢_𝑡 we can derive 𝑣_𝑡(𝑥,𝑡)=𝑣_𝑥𝑥(𝑥,𝑡),0<𝑥<ℓ,0<𝑡<𝑡_𝑓,𝑣(0,𝑡)=𝑢̇₀(𝑡),𝑣(ℓ,𝑡)=𝑢̇_ℓ(𝑡),𝑣_𝑥(0,𝑡)=𝑞̇(𝑡). (110)(111)(112)(113) This is an inverse heat conduction problem (IHCP) for 𝑣(𝑥,𝑡) without using the initial condition.

Therefore as being a numerical method, we can first solve the above IHCP for 𝑣(𝑥,𝑡) by using the MFS in Section 5.4 to obtain a linear equations system and then the method introduced in Section 4 to solve the resultant linear equations system. Thus, we can construct 𝑢(𝑥,𝑡) by 𝑢(𝑥,𝑡)=^𝑡₀𝑣(𝑥,𝜉)𝑑𝜉+𝑓(𝑥), (114) which automatically satisfies the initial condition in (108).

From (114) it follows that 𝑢_𝑥𝑥(𝑥,𝑡)=^𝑡₀𝑣_𝑥𝑥(𝑥,𝜉)𝑑𝜉+𝑓^′′(𝑥), (115) which together with 𝑢_𝑡=𝑣 being inserted into (106) leads to 𝑣(𝑥,𝑡)=^𝑡₀𝑣_𝑥𝑥(𝑥,𝜉)𝑑𝜉+𝑓^′′(𝑥)+𝐻(𝑥). (116) Inserting (110) for 𝑣_𝑥𝑥=𝑣_𝑡 into the above equation and integrating it we can derive the following equation to recover 𝐻(𝑥): 𝐻(𝑥)=𝑣(𝑥,0)−𝑓^′′(𝑥). (117) This approach exhibits threefold ill-posedness: one is the use of the IHCP to solve 𝑣(𝑥,𝑡) which is used to provide the data of 𝑣(𝑥,0) used in (117), one is all the boundary conditions being obtained from the first-order differentials of measured data as shown in (111)–(113), and another is the second-order differential of the data 𝑓(𝑥) in (117).

We consider 𝑢(𝑥,𝑡)=𝑥²(),𝐻(𝑥)+2𝑥𝑡+sin2𝜋𝑥=2𝑥−2+4𝜋²sin(2𝜋𝑥). (118) In (117) we disregard the ill-posedness of 𝑓^′′(𝑥) and suppose that the data 𝑓^′′(𝑥) are given exactly. We solve this problem by the OGSDA with unit subspace method, where we fix 𝑚=10 and 𝛾=0.05. The maximum number of iterations is set to be 100. A random noise with intensity 𝜎=0.05 is added on the data 𝑞(̇𝑡). We have obtained the numerical solution with the steplength being shown in Figure 8(a), which is very small. We compare the heat source computed by the OGSDA with the exact one in Figure 8(b). The numerical error is smaller than 0.568 as shown in Figure 8(c).

(a)

(b)

(c)

Figure 8 

Example 6 showing (a) steplength of OGSDA, (b) comparing numerical and exact solutions, and (c) numerical error of OGSDA.

6. Conclusions

In the present paper, we have derived an optimal algorithm including an 𝑚-vector optimal search direction in an 𝑚-dimensional Krylov subspace or a unit subspace as an optimally generalized steepest-descent algorithm (OGSDA) to solve highly ill-posed linear systems. This algorithm has good computational efficiency and accuracy in solving the ill-posed linear equations for the linear Hilbert problem with 𝑛=300 and three famous linear inverse problems of backward heat conduction problem, inverse Cauchy problem, and inverse heat source problem. In particular, the OGSDA has a better filtering effect against noise, such that for all the numerical examples of inverse problems the numerical results recovered are quite smooth. The robustness of the OGSDA was confirmed by imposing noisy disturbance on the given data even with an intensity being large up to 10%. In all the numerical examples the convergence steps are among 4 to 100 steps, which is very time saving, no one case running over one second in the PC computer. For highly ill-posed problems the performance of OGSDA without needing of large value 𝑚 of the Krylov subspace dimension is better than other algorithms investigated in this paper.

Conflict of Interests

This paper is a purely academic research, and the author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The project NSC-102-2221-E-002-125-MY3 and the 2011 Outstanding Research Award from the National Science Council of Taiwan and the 2011 Taiwan Research Front Award from Thomson Reuters are highly appreciated. The author thanks anonymous referee for the comments which improved the quality of this paper.