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

Authors: Chein-Shan Liu
Journal: Journal of Applied Mathematics
 (2013)

---

## 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 withm-vector descent direction in a Krylov subspace is constructed, of which the m 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.


---

## Body

## 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)Bx=b1,
where x∈ℝn is an unknown vector determined from a given coefficient matrix B∈ℝn×n, which might be unsymmetric, and the input b1∈ℝn which might be polluted by random noise. A measure of the ill-posedness of (1) is the condition number of B:
(2)Cond(B)=∥B∥∥B-1∥,
where ∥B∥ is the Frobenius norm of B.Instead of (1), we can solve a normal linear system:
(3)Cx=b,
where
(4)b=BTb1,C=BTB>0.
Throughout this paper the superscript “T” signifies the transpose.We consider iterative method for solving (3) and define for any vector xk the steepest-descent vector:
(5)rk≔Cxk-b.
Ascher et al. [21] have viewed the gradient-descent method:
(6)xk+1=xk-ηkrk,
as a forward Euler scheme of a time-dependent ODE:
(7)x˙=b-Cx.
The absolute stability bound
(8)ηk≤2maxσ(C)λ
must be obeyed if a uniform stepsize is employed, where σ(C) is the set of all the eigenvalues of C.Specifically, (6) presents a steepest-descent method (SDM), if the steplength is taken to be
(9)ηk=∥rk∥2rkTCrk.
When ∥rk∥ is rather small the calculated rk 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 x0.(ii)
Fork=0,1,…, we repeat the following iterations:
(10)rk=Cxk-b,ηk=∥rk∥2rkTCrk,xk+1=xk-(1-γ)ηkrk.If xk+1 converges according to a given stopping criterion ∥rk+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-r and the m-dimensional search vector in a Krylov subspace. Here we modify the direction r as well as the steplength η from a theoretical optimization being performed on the derived steplength in an m-dimensional Krylov subspace.The remaining parts of this paper are arranged as follows. A numerical iterative scheme with-u as a descent direction is proposed, and the convergence to the minimal point x*=C-1b 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 m-dimensional Krylov subspace. Then, a genuine algorithm is constructed in Section 4, resulting in an optimally generalized steepest-descent algorithm (OGSDA) in terms of m 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:
(11)minx∈ℝnφ(x)=minx∈ℝn[12xTCx-bTx].
In the SDM, we search the next x(t+Δt) from the current state x(t) by minimizing the functional φ along the direction -r(t); that is,
(12)minηφ(x(t)-ηr(t)).Through some calculations we can obtain(13)η=∥r(t)∥2rT(t)Cr(t).
Thus we have the following iterative algorithm of SDM:
(14)x(t+Δt)=x(t)-∥r(t)∥2rT(t)Cr(t)r(t).
A straightforward calculation reveals that
(15)φ(x(t+Δt))<φ(x(t)).
It means that the iterative sequences of SDM converge to the minimal point.Then, we consider a more general scheme with the descent direction-u(t); that is,
(16)x(t+Δt)=x(t)-r(t)·u(t)uT(t)Cu(t)u(t).
Similarly we have the following result.Lemma 1.
Equation (16) is a result of the following minimization along the descent direction -u:
(17)minηφ(x(t)-ηu(t)).
The functional φ is decreasing stepwise and the iterative scheme (16) is convergent to the minimal point x*=C-1b.Proof.
Inserting(18)x(t+Δt)=x(t)-ηu(t)
into the functional φ we have
(19)φ(t+Δt)=12[xTCx-2ηuTCx+η2uTCu]-bTx+ηuTb,
where we simplify the notation of φ(x(t+Δt)) by φ(t+Δt) and omit the terms (t) in x(t) and u(t). Taking the derivative of (19) with respect to η and equating it to zero we can derive
(20)η=uTruTCu.
Because η is a steplength, we demand it to be positive; that is, r·u>0. By feeding the above η into (18) we can derive the iterative scheme (16).
Inserting (20) into (19) and subtracting the resultant by
(21)φ(t)=12xTCx-bTx,
we can obtain
(22)φ(t+Δt)-φ(t)=12η2uTCu-ηuTr=-(uTr)22uTCu.
It implies
(23)φ(t+Δt)=φ(t)-(uTr)22uTCu<φ(t),
and the functional φ is strictly decreasing. Because the minimal point x*=C-1b 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 x(t) 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 Lemma1 we have proven that the new scheme (16) is convergent, but it cannot tell us what descent direction -u is the best. Although the SDM uses the steepest-descent direction -r, 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
(24)xk+1=xk-(Δrk-1)TΔxk-1∥Δrk-1∥2rk,
where Δrk-1=rk-rk-1 and Δxk-1=xk-xk-1, with initial guesses r0=0 and x0=0. 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)
Selectm and give an initial x0.(ii)
Fork=0,1,…, we repeat the following computations:
(25)rk=b1-Bxk,Arnoldiproceduretosetupvjk,j=1,…,m,(fromv1k=rk∥rk∥),Vk=[v1k,…,vmk],Solve(H¯kTH¯k)yk=∥rk∥H¯kTe1,obtainingyk,zk=Vkyk,xk+1=xk+zk.If xk+1 converges according to a given stopping criterion ∥rk+1∥<ɛ, then stop; otherwise, go to step (ii). In the above, Vk is an n  ×  m matrix, while H¯k is an augmented Heissenberg upper triangular matrix with dimension (m  +  1)  ×  m, and e1 is the first column of Im+1.


### 3.1. A Generalized SDM

We suppose that the original steepest-descent direction-r in the SDM is replaced by a new descent vector -u:
(26)u=r+∑k=1mαkvk,
which is an optimal combination of r and the m-vector vk,k=1,…,m. The coefficients αk will be optimized in Section 3.2.Now we describe how to set up them-vector vk,k=1,…,m, by the Krylov subspace method. Suppose that we have an m-dimensional Krylov subspace generated by the coefficient matrix C from the steepest-descent vector r:
(27)𝒦m≔span{Cr,…,Cmr}.
Then the Arnoldi process is used to normalize and orthogonalize the Krylov vectors Cjr,j=1,…,m, such that the resultant vectors vi, i=1,…,m, satisfy vi·vj=δij,i,j=1,…,m, where δij is the Kronecker delta tensor.The other way to obtainvi,i=1,…,m, is simply by using the unit vectors in the Euclidean space, where vi is the ith column vector of the identity matrix In. The resultant m-dimensional space is named the unit subspace.


### 3.2. Optimization of Algorithm

LetJ be an n×m matrix, which consists of
(28)J≔[v1,…,vm].
That is, the jth column of J is the vector vj. Because v1,…,vm are linearly independent vectors and m<n, the expansion matrix J has a full rank with rank(J)=m. Then, (26) can be written as
(29)u=r+Jα,
where α≔(α1,…,αm)T.Now, in view of (16) the new steplength becomes
(30)η=rTuuTCu.
When u=r, 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:
(31)max{η=rTuuTCu}.
Here, we use the steplength η as a merit function to be optimized. Hence, we can prove the following result.Theorem 2.
Forrank(J)=m and u∈r+𝒦m, the descent direction -u which maximizes (31) is given by
(32)u=r-ECr+λEr,A=JTCJ,(33)E=JA-1JT,b0=rT(CEC-C)r<0,(34)b1=rT(In-CE)r>0,(35)b2=rTEr>0,(36)λ=b12-b0b2-b1b2.
Moreover, the steplength η defined in (30) is positive; that is,
(37)η=12λ=b22b12-b0b2-2b1>0.Proof.
We divide this proof into three parts.(A)
A Minimization Problem. Inserting (29) into (31) we can obtain the following minimization problem:
(38)minα1,…,αmu·(Cu)r·u,where
(39)r·u=∥r∥2+rTJα,(40)u·(Cu)=rTCr+2rTCJα+αTAα,
in which
(41)A≔JTCJ
is an m×m positive definite matrix. Because the matrix C has a full rank with rank(C)=n and positive definite and J has a full rank with rank(J)=m, from the above equation it follows that A has a full rank with rank(A)=m. Then, we confirm that A is an m×m positive definite matrix by (41).
From (39) and (40) we can derive
(42)y1≔∇αr·u=JTr,(43)y2≔∇αu·(Cu)=2JTCr+2Aα,
where ∇α denotes the gradient with respect to α.
By using(44)∇αu·(Cu)r·u=0⟹r·u∇αu·(Cu)-u·(Cu)∇αr·u=0,
we can derive the following governing equation of vector form to solve the optimal parameter α:
(45)r·uy2-u·(Cu)y1=0.
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 Jα:
(46)Jα=2r·uA+y2y1T-y1y2T.
It is an m×m positive definite matrix, because for any z≠0∈ℝm we have
(47)zTJαz=2r·uzTAz>0,
where we have used the positive definite property of A, and r·u>0. Because η defined by (20) is a steplength, we demand it to be positive; that is, r·u>0 (see the proof given below for r·u>0).(B)
A Quadratic Equation. Equations (40) and (43) can be written as
(48)u·(Cu)=rTCr+2y0TJα+αTAα,(49)y2=2JTy0+2Aα,where
(50)y0≔Cr.
From (45) we can observe that y2 is proportional to y1, which is supposed to be
(51)y2=u·(Cu)r·uy1=2λy1,
where 2λ is a multiplier to be determined. From (51), two equations follow
(52)y2=2λy1,(53)u·(Cu)=2λr·u.
Then, by (42), (49), and (52) we can solve α by
(54)α=λA-1JTr-A-1JTy0.
Inserting (54) for α into (39) and (48) we have that
(55)r·u=∥r∥2-y0TEr+λrTEr,u·(Cu)=rTCr-y0TEy0+λ2rTEr,
where
(56)E≔JA-1JT
is an n×n positive semi-definite matrix.
Now, from (53) and (55) we can derive a quadratic equation to solve λ:
(57)a2λ2+a1λ+a0=0,
where
(58)a0≔y0TEy0-rTCr,a1≔2∥r∥2-2y0TEr,a2≔rTEr.(C)
A Closed-Form Solution of Optimal Parameterα. Now we prove a0<0. By using (50) we have
(59)a0=rT[CEC-C]r<0,if we can prove
(60)CEC-C<0.
From (56) and (41) it follows that
(61)ECJ=J(JTCJ)-1JTCJ=J,(62)JTCE=JTCJ(JTCJ)-1JT=JT,(63)JT[CEC]J=JT[CJ(JTCJ)-1JTC]J=JTCJ.
They indicate that E is a Penrose pseudoinverse of C in the Krylov subspace; hence, (60) follows.
According to (58), it is obvious that a2>0; hence, by (59) we have
(64)a12-4a0a2>0.
As a consequence, (57) has a closed-form positive real solution:
(65)λ=a12-4a0a2-a12a2.
Inserting it into (54) we can obtain the closed-form solution of α:
(66)α=a12-4a0a2-a12a2A-1JTr-A-1JTy0.
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 r·u>0.Remark 3.
Because ofa2>0 and a0<0 in (57), it has a positive solution as shown in (65) and in addition a negative solution as follows:
(67)λ=-a12-4a0a2-a12a2,
where
(68)a1=2rT(In-CE)r>0,
in view of (58), (50), and (62). Hence, by (30) and (37) we have
(69)r·u=∥r∥2-y0TEr+λrTEr<0,
which contradicts (20), where η is a steplength being positive. Below we will give a numerical example to demonstrate this fact.Corollary 4.
Ifm=n, then one has
(70)α=J-1C-1r-J-1r.Proof.
Ifm=n, then J defined by (28) is an n×n nonsingular matrix. Simultaneously, E defined by (56) is equal to C-1, a0=0, a1=0, a2>0, and meanwhile (66) reduces to (70), where we have taken λ=1 according to the above condition.Remark 5.
Ifm=n, by (29) and (70) we have
(71)u=r+Jα=r+J[J-1C-1r-J-1r]=C-1r,
and inserting it into (31) we can prove η=1. In any way, under this condition we have the fastest convergent behavior, because the solution x=C-1b is already obtained. Usually, we cannot take m=n, which leads to a quite difficult task to find the inverse of A=JTCJ when m=n is a large number.


## 3.1. A Generalized SDM

We suppose that the original steepest-descent direction-r in the SDM is replaced by a new descent vector -u:
(26)u=r+∑k=1mαkvk,
which is an optimal combination of r and the m-vector vk,k=1,…,m. The coefficients αk will be optimized in Section 3.2.Now we describe how to set up them-vector vk,k=1,…,m, by the Krylov subspace method. Suppose that we have an m-dimensional Krylov subspace generated by the coefficient matrix C from the steepest-descent vector r:
(27)𝒦m≔span{Cr,…,Cmr}.
Then the Arnoldi process is used to normalize and orthogonalize the Krylov vectors Cjr,j=1,…,m, such that the resultant vectors vi, i=1,…,m, satisfy vi·vj=δij,i,j=1,…,m, where δij is the Kronecker delta tensor.The other way to obtainvi,i=1,…,m, is simply by using the unit vectors in the Euclidean space, where vi is the ith column vector of the identity matrix In. The resultant m-dimensional space is named the unit subspace.


## 3.2. Optimization of Algorithm

LetJ be an n×m matrix, which consists of
(28)J≔[v1,…,vm].
That is, the jth column of J is the vector vj. Because v1,…,vm are linearly independent vectors and m<n, the expansion matrix J has a full rank with rank(J)=m. Then, (26) can be written as
(29)u=r+Jα,
where α≔(α1,…,αm)T.Now, in view of (16) the new steplength becomes
(30)η=rTuuTCu.
When u=r, 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:
(31)max{η=rTuuTCu}.
Here, we use the steplength η as a merit function to be optimized. Hence, we can prove the following result.Theorem 2.
Forrank(J)=m and u∈r+𝒦m, the descent direction -u which maximizes (31) is given by
(32)u=r-ECr+λEr,A=JTCJ,(33)E=JA-1JT,b0=rT(CEC-C)r<0,(34)b1=rT(In-CE)r>0,(35)b2=rTEr>0,(36)λ=b12-b0b2-b1b2.
Moreover, the steplength η defined in (30) is positive; that is,
(37)η=12λ=b22b12-b0b2-2b1>0.Proof.
We divide this proof into three parts.(A)
A Minimization Problem. Inserting (29) into (31) we can obtain the following minimization problem:
(38)minα1,…,αmu·(Cu)r·u,where
(39)r·u=∥r∥2+rTJα,(40)u·(Cu)=rTCr+2rTCJα+αTAα,
in which
(41)A≔JTCJ
is an m×m positive definite matrix. Because the matrix C has a full rank with rank(C)=n and positive definite and J has a full rank with rank(J)=m, from the above equation it follows that A has a full rank with rank(A)=m. Then, we confirm that A is an m×m positive definite matrix by (41).
From (39) and (40) we can derive
(42)y1≔∇αr·u=JTr,(43)y2≔∇αu·(Cu)=2JTCr+2Aα,
where ∇α denotes the gradient with respect to α.
By using(44)∇αu·(Cu)r·u=0⟹r·u∇αu·(Cu)-u·(Cu)∇αr·u=0,
we can derive the following governing equation of vector form to solve the optimal parameter α:
(45)r·uy2-u·(Cu)y1=0.
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 Jα:
(46)Jα=2r·uA+y2y1T-y1y2T.
It is an m×m positive definite matrix, because for any z≠0∈ℝm we have
(47)zTJαz=2r·uzTAz>0,
where we have used the positive definite property of A, and r·u>0. Because η defined by (20) is a steplength, we demand it to be positive; that is, r·u>0 (see the proof given below for r·u>0).(B)
A Quadratic Equation. Equations (40) and (43) can be written as
(48)u·(Cu)=rTCr+2y0TJα+αTAα,(49)y2=2JTy0+2Aα,where
(50)y0≔Cr.
From (45) we can observe that y2 is proportional to y1, which is supposed to be
(51)y2=u·(Cu)r·uy1=2λy1,
where 2λ is a multiplier to be determined. From (51), two equations follow
(52)y2=2λy1,(53)u·(Cu)=2λr·u.
Then, by (42), (49), and (52) we can solve α by
(54)α=λA-1JTr-A-1JTy0.
Inserting (54) for α into (39) and (48) we have that
(55)r·u=∥r∥2-y0TEr+λrTEr,u·(Cu)=rTCr-y0TEy0+λ2rTEr,
where
(56)E≔JA-1JT
is an n×n positive semi-definite matrix.
Now, from (53) and (55) we can derive a quadratic equation to solve λ:
(57)a2λ2+a1λ+a0=0,
where
(58)a0≔y0TEy0-rTCr,a1≔2∥r∥2-2y0TEr,a2≔rTEr.(C)
A Closed-Form Solution of Optimal Parameterα. Now we prove a0<0. By using (50) we have
(59)a0=rT[CEC-C]r<0,if we can prove
(60)CEC-C<0.
From (56) and (41) it follows that
(61)ECJ=J(JTCJ)-1JTCJ=J,(62)JTCE=JTCJ(JTCJ)-1JT=JT,(63)JT[CEC]J=JT[CJ(JTCJ)-1JTC]J=JTCJ.
They indicate that E is a Penrose pseudoinverse of C in the Krylov subspace; hence, (60) follows.
According to (58), it is obvious that a2>0; hence, by (59) we have
(64)a12-4a0a2>0.
As a consequence, (57) has a closed-form positive real solution:
(65)λ=a12-4a0a2-a12a2.
Inserting it into (54) we can obtain the closed-form solution of α:
(66)α=a12-4a0a2-a12a2A-1JTr-A-1JTy0.
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 r·u>0.Remark 3.
Because ofa2>0 and a0<0 in (57), it has a positive solution as shown in (65) and in addition a negative solution as follows:
(67)λ=-a12-4a0a2-a12a2,
where
(68)a1=2rT(In-CE)r>0,
in view of (58), (50), and (62). Hence, by (30) and (37) we have
(69)r·u=∥r∥2-y0TEr+λrTEr<0,
which contradicts (20), where η is a steplength being positive. Below we will give a numerical example to demonstrate this fact.Corollary 4.
Ifm=n, then one has
(70)α=J-1C-1r-J-1r.Proof.
Ifm=n, then J defined by (28) is an n×n nonsingular matrix. Simultaneously, E defined by (56) is equal to C-1, a0=0, a1=0, a2>0, and meanwhile (66) reduces to (70), where we have taken λ=1 according to the above condition.Remark 5.
Ifm=n, by (29) and (70) we have
(71)u=r+Jα=r+J[J-1C-1r-J-1r]=C-1r,
and inserting it into (31) we can prove η=1. In any way, under this condition we have the fastest convergent behavior, because the solution x=C-1b is already obtained. Usually, we cannot take m=n, which leads to a quite difficult task to find the inverse of A=JTCJ when m=n is a large number.


## 4. An Optimal Generalized Steepest-Descent Algorithm

Thus, we arrive at the followingoptimally generalized steepest-descent algorithm (OGSDA).(i)
Selectm and 0≤γ<1, and give an initial value of x0.(ii)
Fork=0,1,…, we repeat the following computations:
(72)rk=Cxk-b,Arnoldiproceduretosetupvjk,j=1,…,m,(fromv1k=Crk∥Crk∥),Jk=[v1k,…,vmk],Ak=JkTCJk,Ek=JkAk-1JkT,Compute(36),obtainingλk,uk=rk-EkCrk+λkEkrk,ηk=12λk,xk+1=xk-(1-γ)ηkuk.If xk+1 converges according to a given stopping criterion ∥rk+1∥<ɛ, then stop; otherwise, go to step (ii).Remark 6.
Depending on the different choice ofvj,j=1,…,m, the optimal algorithm is further classified into two types. When vj,j=1,…,m, 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 vj is simply the jth column of the unit matrix In, 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 J 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 A=JTCJ to obtain A-1 and E. On the other hand, in the unit subspace method, the expansion matrix J is fixed, and we only need to calculate A=JTCJ, A-1 and E 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 solutionxk by OGSDA with k steps. We can estimate the cost of the computation of the scalar λ, which at each iterative step requires an inversion of m×m matrix, the matrix-vector multiplications of n×n-matrix and n-vector three times, and the inner products of two n-vectors five times. This portion is very time saving with totally 3kn+km2, because of m≪n. The computation of the Arnoldi vectors needs totally k(k+1)n+kn2 multiplications. In the computation of u it requires 2kn multiplications. The k steps OGSDA therefore requires totally k(k+6)n+kn2+km2 multiplications. Dividing by the total number of steps k we can obtain that each step requires (k+6)n+n2+m2 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:(73)min{f=12xTPx+qTx},(74)Qx=b0,
where P is an n1×n1 matrix, Q is an m1×n1 matrix, and b0 is an m1-vector, which means that (74) provides m1 linear constraints. According to the Lagrange theory we need to solve (1) with the following b and B:
(75)b1:=[-qb0],B=[PQTQ0m1×m1].For definite we taken1=3 and m1=2 with
(76)min{f=x12+2x22+x32-2x1x2+x3},x1+x2+x3=4,2x1-x2+x3=2.
Under the following parameters m=2, γ=0.2, and ɛ=10-5 and starting from the initial values of x1=x2=x3=1 and two Lagrange multipliers λ1=-2 and λ2=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 α1 and α2 and the values of ηk and a0k. It can be seen that a0k∈[-22.16,-5.72×10-13], as proven by (59) a0k<0. Now we raise m to m=n=5 and also the convergence criterion to ɛ=10-5 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 minf=3.977273 at the point (x1,x2,x3)=(1.909090909092,1.954545454479,0.13636363636).Example 1 showing (a) residual, ((b) and (c)) optimal parameters, and (d) steplength anda0.
(a)
(b)
(c)
(d)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 m=2 and γ=0.2. First the algorithm does not converge within 500 steps under a weak convergence criterion with ɛ=10-2 as shown in Figure 2(a). Then the algorithm gives negative steplength of η as shown in Figure 2(b).For Example 1, using (67) the algorithm does not converge and gives a negative steplength.
(a)
(b)


### 5.2. Example 2

Finding ann-order polynomial function p(x)=a0+a1x+⋯+anxn to best match a continuous function f(x) in the interval of x∈[0,1](77)mindeg(p)≤n∫01[f(x)-p(x)]2dx,
leads to a problem governed by (3), where C is the (n+1)×(n+1) Hilbert matrix defined by
(78)Cij=1i+j-1,x is composed of the n+1 coefficients a0,a1,…,an that appeared in p(x), and
(79)b=[∫01f(x)dx∫01xf(x)dx⋮∫01xnf(x)dx]
is uniquely determined by the function f(x).The Hilbert matrix is a notorious example of highly ill-conditioned matrices. Equation (3) with the matrix C 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 C given by (78). The ill-posedness of (3) increases fast with n. We consider an exact solution with xj=1 and bi is given by
(80)bi=∑j=1n‍1i+j-1+σR(i),
with R(i) being random numbers between [-1,1].Liu [22] has applied the relaxed steepest-descent method (RSDM) with γ=0.06, starting from x1=⋯=x9=0.5, to solve the Hilbert problem with n=9 with a stopping criterion ɛ=10-8 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-3. When we apply the OGSDA with the unit subspace method under m=n=9 and γ=10-5 to solve this problem it converges only three steps and with the maximum error being 4.68×10-6. It is very time saving because we only need to calculate A-1 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 ɛ1=10-5 to find A-1. 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.
Solutions
x
1
x
2
x
3
x
4
x
5
x
6
x
7
x
8
x
9
Exact
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
SVD
0.9999
1.008
0.985
0.995
1.007
1.012
1.009
0.999
0.984
SDM
1.00001
0.99980
1.00091
0.99910
0.99928
1.00035
1.00103
1.00062
0.99898
CGM
1.00000
1.00001
0.99924
1.00019
0.99988
0.99985
1.00008
1.00019
0.99987
RSDM
1.00000
0.99986
1.00087
0.99888
0.99928
1.00058
1.00128
1.00066
0.99856
OGSDA
0.9999999999
1.000000005
0.9999999185
1.000000569
0.9999979488
1.000004126
0.9999953195
1.000002797
0.9999993151As shown in Corollary4 when m is taken to be m=n we have the fastest convergent solution. Usually, we can take a smaller m for the easy solution of A-1. For example, if we take m=5 in the OGSDA with the Krylov subspace method, we can obatin the solution with the maximum error being 4.45×10-4 through only four steps. Under the same conditions, the CGM is convergent with six steps with the maximum error being 7.857×10-3.Next we consider a more difficult case withn=300. In the computation a random noise in the level of σ=10-6 is imposed on the right-hand side, and the convergence criteria used in CGM and OGSDA are both ɛ=10-2 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 m=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 m=9 converges with only two steps, CGM with four steps, and the BBM with ten steps.Example 2 of the Hilbert problem withn=300, showing (a) relative residuals, (b) steplength of OGSDA, and (c) comparing the numerical errors of CGM, GMRES, OGSDA, and BBM.
(a)
(b)
(c)


### 5.3. Example 3

When we apply a central difference scheme to the following two-point boundary value problem:(81)-u′′(x)=f(x),0<x<1,u(0)=a,u(1)=b,
we can derive a linear equations system
(82)Cu=[2-1-12-1·········-12][u1u2⋮un]=[(Δx)2f(Δx)+a(Δx)2f(2Δx)⋮(Δx)2f((n-1)Δx)(Δx)2f(nΔx)+b],
where Δx=1/(n+1) is the spatial length and ui=u(iΔx),i=1,…,n, are unknown values of u(x) at the grid points xi=iΔx. u0=a and un+1=b are the given boundary conditions. The above matrix C is known as a central difference matrix.The eigenvalues ofC are found to be
(83)4sin2kπ2(n+1),k=1,2,…,n,
which together with the symmetry of C indicates that C is positive definite, and
(84)Cond(C)=sin2(nπ/2(n+1))sin2(π/2(n+1))
is a large number when the grid number n is very large.In this numerical test we fixn=200 and thus the condition number of C is 16210.7. Let us consider the boundary value problem in (81) with f(x)=sinπx. The exact solution is
(85)u(x)=a+(b-a)x+1π2sinπx,
where we fix a=1 and b=2.A relative random noise with intensityσ=0.01 is added by f^(iΔx)=f(iΔx)[1+σR(i)] into the right-hand side data of (82). Under a moderate convergence criterion with ɛ=10-7 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-5. Then we apply the OGSDA with γ=0.25 and m=30 to this problem, whose number of iterations is 66, and the maximum error is 1.9×10-5. The number of iterations of GMRES with m=10 is 811, and the maximum error is 4.21×10-5, while the number of iterations of BBM is 1708, and the maximum error is 1.06×10-5.Example 3, comparing (a) relative residuals, (b) convergence rates, and (c) numerical errors of CGM, OGSDA, GMRES, and BBM.
(a)
(b)
(c)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
(86)ConvergenceRateOGSDA≔1sk,sk=1-(1-γ2)[(Cxk-b)·Cuk]2∥Cxk-b∥2∥Cuk∥2,ConvergenceRateCGM≔1sk,sk=1-ηk∥Bxk-b1∥2,
where ηk=∥rk∥2/pkTCpk, in which pk 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 of0<x<ℓ by subjecting to the boundary conditions at two ends of a slab:
(87)ut(x,t)=κuxx(x,t),0<t<T,0<x<ℓ,(88)u(0,t)=u0(t),u(ℓ,t)=uℓ(t),
we solve u under a final time condition:
(89)u(x,T)=uT(x).The fundamental solution of (87) is given by
(90)K(x,t)=H(t)2κπtexp(-x24κt),
where H(t) 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 ofu at the field point z=(x,t) can be expressed as a linear combination of the fundamental solutions U(z,sj):
(91)u(z)=∑j=1ncjU(z,sj),sj=(ηj,τj)∈Ωc,
where n is the number of source points, cj are unknown coefficients, and sj are source points being located in the complement Ωc of Ω=[0,ℓ]×[0,T]. For the heat conduction equation we have the basis functions:
(92)U(z,sj)=K(x-ηj,t-τj).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 thet-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)Bx=b1,
where
(94)Bij=U(zi,sj),x=(c1,…,cn)T,b1=(uℓ(ti),i=1,…,m1;uT(xj),j=1,…,m2;u0(tk),k=m1,…,1(uℓ(ti),i=1,…,m1;uT(xj),)T,
and n=2m1+m2.Since the BHCP is highly ill-posed, the ill-condition of the coefficient matrixB 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:
(95)u(x,t)=cos(πx)exp(-π2t).
For the case with T=1 the value of final time data is in the order of 10-4, which is small by comparing with the value of the initial temperature f(x)=u0(x)=cos(πx) to be retrieved, which is O(1). First we impose a relative random noise with an intensity σ=10% being imposed on the final time data. Under the following parameters m1=11, m2=6, m=10, γ=0.25, and ɛ=ɛ1=10-5, 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-4, which is better than 1.247×10-3 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 m=10 to solve the BHCP, of which both algorithms cannot converge to satisfy the convergence criterion with ɛ=10-5, 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.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.
(a)
(b)
(c)
(d)


### 5.5. Example 5

Let us consider the inverse Cauchy problem for the Laplace equation:(96)Δu=urr+1rur+1r2uθθ=0,(97)u(ρ,θ)=h(θ),0≤θ≤π,(98)un(ρ,θ)=g(θ),0≤θ≤π,
where h(θ) and g(θ) are given function. The inverse Cauchy problem is specified as follows.To seek an unknown boundary functionf(θ) on the part Γ2:={(r,θ)∣r=ρ(θ),π<θ<2π} of the boundary under (96)–(98) with the overspecified data being given on Γ1≔{(r,θ)∣r=ρ(θ),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 ofu at the field point x=(rcosθ,rsinθ) can be expressed as a linear combination of fundamental solutions U(x,sj):
(99)u(x)=∑j=1ncjU(x,sj),sj∈Ωc.
For the Laplace equation (96) we have the fundamental solutions:
(100)U(x,sj)=lnrj,rj=∥x-sj∥.In the practical application of MFS, by imposing the boundary conditions (97) and (98) at N points on (99) we can obtain a linear equations system:
(101)Bc=b1,
where
(102)xi=(xi1,xi2)=(ρ(θi)cosθi,ρ(θi)sinθi),sj=(sj1,sj2)=(R(θj)cosθj,R(θj)sinθj),Bij=ln∥xi-sj∥,ifiisodd,Bij=η(θi)∥xi-sj∥2(ρ(θi)-sj1cosθi-sj2sinθi-ρ′(θi)ρ(θi)cccccccccccccc×[sj1sinθi-sj2cosθi](ρ(θi)-sj1cosθi-sj2sinθi-ρ′(θi)ρ(θi)),ifiiseven,c=(c1,…,cn)T,b1=(h(θ1),g(θ1),…,h(θN),g(θN))T,
in which n=2N, and
(103)η(θ)=ρ(θ)ρ2(θ)+[ρ′(θ)]2.
The above R(θ)=ρ(θ)+D with an offset D can be used to locate the source points along a contour with a radius R(θ). 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:(104)u(x,y)=cosxcoshy+sinxsinhy,
defined in a domain with a complex amoeba-like irregular shape as a boundary:
(105)ρ(θ)=exp(sinθ)sin2(2θ)+exp(cosθ)cos2(2θ).
After imposing the boundary conditions (97) and (98) at N points on (99) we can obtain a linear equations system. The noise being imposed on the measured data h and g is σ=0.01.We solve this problem by the OGSDA withm=10 and γ=0.2. Through 15 iterations the residual is smaller than the convergence criterion ɛ=10-5 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.Example 5 showing (a) residual and steplength of OGSDA, (b) comparing numerical and exact solutions, and (c) numerical errors of OGSDA, GMRES, and BBM.
(a)
(b)
(c)We also apply the GMRES withm=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 functionH(x) for a one-dimensional heat conduction equation:
(106)ut(x,t)=uxx(x,t)+H(x),0<x<ℓ,0<t<tf,(107)u(0,t)=u0(t),u(ℓ,t)=uℓ(t),(108)u(x,0)=f(x).
In order to identify H(x) we impose a Neumann type boundary condition:
(109)ux(0,t)=q(t).We propose a numerical differential method by lettingv=ut. Taking the differentials of (106), (107), and (109) with respect to t and letting v=ut we can derive
(110)vt(x,t)=vxx(x,t),0<x<ℓ,0<t<tf,(111)v(0,t)=u˙0(t),(112)v(ℓ,t)=u˙ℓ(t),(113)vx(0,t)=q˙(t).
This is an inverse heat conduction problem (IHCP) for v(x,t) without using the initial condition.Therefore as being a numerical method, we can first solve the above IHCP forv(x,t) 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 u(x,t) by
(114)u(x,t)=∫0tv(x,ξ)dξ+f(x),
which automatically satisfies the initial condition in (108).From (114) it follows that
(115)uxx(x,t)=∫0tvxx(x,ξ)dξ+f′′(x),
which together with ut=v being inserted into (106) leads to
(116)v(x,t)=∫0tvxx(x,ξ)dξ+f′′(x)+H(x).
Inserting (110) for vxx=vt into the above equation and integrating it we can derive the following equation to recover H(x):
(117)H(x)=v(x,0)-f′′(x).
This approach exhibits threefold ill-posedness: one is the use of the IHCP to solve v(x,t) which is used to provide the data of v(x,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 f(x) in (117).We consider(118)u(x,t)=x2+2xt+sin(2πx),H(x)=2x-2+4π2sin(2πx).
In (117) we disregard the ill-posedness of f′′(x) and suppose that the data f′′(x) are given exactly. We solve this problem by the OGSDA with unit subspace method, where we fix m=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 q˙(t). 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).Example 6 showing (a) steplength of OGSDA, (b) comparing numerical and exact solutions, and (c) numerical error of OGSDA.
(a)
(b)
(c)


## 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:(73)min{f=12xTPx+qTx},(74)Qx=b0,
where P is an n1×n1 matrix, Q is an m1×n1 matrix, and b0 is an m1-vector, which means that (74) provides m1 linear constraints. According to the Lagrange theory we need to solve (1) with the following b and B:
(75)b1:=[-qb0],B=[PQTQ0m1×m1].For definite we taken1=3 and m1=2 with
(76)min{f=x12+2x22+x32-2x1x2+x3},x1+x2+x3=4,2x1-x2+x3=2.
Under the following parameters m=2, γ=0.2, and ɛ=10-5 and starting from the initial values of x1=x2=x3=1 and two Lagrange multipliers λ1=-2 and λ2=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 α1 and α2 and the values of ηk and a0k. It can be seen that a0k∈[-22.16,-5.72×10-13], as proven by (59) a0k<0. Now we raise m to m=n=5 and also the convergence criterion to ɛ=10-5 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 minf=3.977273 at the point (x1,x2,x3)=(1.909090909092,1.954545454479,0.13636363636).Example 1 showing (a) residual, ((b) and (c)) optimal parameters, and (d) steplength anda0.
(a)
(b)
(c)
(d)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 m=2 and γ=0.2. First the algorithm does not converge within 500 steps under a weak convergence criterion with ɛ=10-2 as shown in Figure 2(a). Then the algorithm gives negative steplength of η as shown in Figure 2(b).For Example 1, using (67) the algorithm does not converge and gives a negative steplength.
(a)
(b)


## 5.2. Example 2

Finding ann-order polynomial function p(x)=a0+a1x+⋯+anxn to best match a continuous function f(x) in the interval of x∈[0,1](77)mindeg(p)≤n∫01[f(x)-p(x)]2dx,
leads to a problem governed by (3), where C is the (n+1)×(n+1) Hilbert matrix defined by
(78)Cij=1i+j-1,x is composed of the n+1 coefficients a0,a1,…,an that appeared in p(x), and
(79)b=[∫01f(x)dx∫01xf(x)dx⋮∫01xnf(x)dx]
is uniquely determined by the function f(x).The Hilbert matrix is a notorious example of highly ill-conditioned matrices. Equation (3) with the matrix C 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 C given by (78). The ill-posedness of (3) increases fast with n. We consider an exact solution with xj=1 and bi is given by
(80)bi=∑j=1n‍1i+j-1+σR(i),
with R(i) being random numbers between [-1,1].Liu [22] has applied the relaxed steepest-descent method (RSDM) with γ=0.06, starting from x1=⋯=x9=0.5, to solve the Hilbert problem with n=9 with a stopping criterion ɛ=10-8 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-3. When we apply the OGSDA with the unit subspace method under m=n=9 and γ=10-5 to solve this problem it converges only three steps and with the maximum error being 4.68×10-6. It is very time saving because we only need to calculate A-1 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 ɛ1=10-5 to find A-1. 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.
Solutions
x
1
x
2
x
3
x
4
x
5
x
6
x
7
x
8
x
9
Exact
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
SVD
0.9999
1.008
0.985
0.995
1.007
1.012
1.009
0.999
0.984
SDM
1.00001
0.99980
1.00091
0.99910
0.99928
1.00035
1.00103
1.00062
0.99898
CGM
1.00000
1.00001
0.99924
1.00019
0.99988
0.99985
1.00008
1.00019
0.99987
RSDM
1.00000
0.99986
1.00087
0.99888
0.99928
1.00058
1.00128
1.00066
0.99856
OGSDA
0.9999999999
1.000000005
0.9999999185
1.000000569
0.9999979488
1.000004126
0.9999953195
1.000002797
0.9999993151As shown in Corollary4 when m is taken to be m=n we have the fastest convergent solution. Usually, we can take a smaller m for the easy solution of A-1. For example, if we take m=5 in the OGSDA with the Krylov subspace method, we can obatin the solution with the maximum error being 4.45×10-4 through only four steps. Under the same conditions, the CGM is convergent with six steps with the maximum error being 7.857×10-3.Next we consider a more difficult case withn=300. In the computation a random noise in the level of σ=10-6 is imposed on the right-hand side, and the convergence criteria used in CGM and OGSDA are both ɛ=10-2 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 m=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 m=9 converges with only two steps, CGM with four steps, and the BBM with ten steps.Example 2 of the Hilbert problem withn=300, showing (a) relative residuals, (b) steplength of OGSDA, and (c) comparing the numerical errors of CGM, GMRES, OGSDA, and BBM.
(a)
(b)
(c)


## 5.3. Example 3

When we apply a central difference scheme to the following two-point boundary value problem:(81)-u′′(x)=f(x),0<x<1,u(0)=a,u(1)=b,
we can derive a linear equations system
(82)Cu=[2-1-12-1·········-12][u1u2⋮un]=[(Δx)2f(Δx)+a(Δx)2f(2Δx)⋮(Δx)2f((n-1)Δx)(Δx)2f(nΔx)+b],
where Δx=1/(n+1) is the spatial length and ui=u(iΔx),i=1,…,n, are unknown values of u(x) at the grid points xi=iΔx. u0=a and un+1=b are the given boundary conditions. The above matrix C is known as a central difference matrix.The eigenvalues ofC are found to be
(83)4sin2kπ2(n+1),k=1,2,…,n,
which together with the symmetry of C indicates that C is positive definite, and
(84)Cond(C)=sin2(nπ/2(n+1))sin2(π/2(n+1))
is a large number when the grid number n is very large.In this numerical test we fixn=200 and thus the condition number of C is 16210.7. Let us consider the boundary value problem in (81) with f(x)=sinπx. The exact solution is
(85)u(x)=a+(b-a)x+1π2sinπx,
where we fix a=1 and b=2.A relative random noise with intensityσ=0.01 is added by f^(iΔx)=f(iΔx)[1+σR(i)] into the right-hand side data of (82). Under a moderate convergence criterion with ɛ=10-7 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-5. Then we apply the OGSDA with γ=0.25 and m=30 to this problem, whose number of iterations is 66, and the maximum error is 1.9×10-5. The number of iterations of GMRES with m=10 is 811, and the maximum error is 4.21×10-5, while the number of iterations of BBM is 1708, and the maximum error is 1.06×10-5.Example 3, comparing (a) relative residuals, (b) convergence rates, and (c) numerical errors of CGM, OGSDA, GMRES, and BBM.
(a)
(b)
(c)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
(86)ConvergenceRateOGSDA≔1sk,sk=1-(1-γ2)[(Cxk-b)·Cuk]2∥Cxk-b∥2∥Cuk∥2,ConvergenceRateCGM≔1sk,sk=1-ηk∥Bxk-b1∥2,
where ηk=∥rk∥2/pkTCpk, in which pk 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 of0<x<ℓ by subjecting to the boundary conditions at two ends of a slab:
(87)ut(x,t)=κuxx(x,t),0<t<T,0<x<ℓ,(88)u(0,t)=u0(t),u(ℓ,t)=uℓ(t),
we solve u under a final time condition:
(89)u(x,T)=uT(x).The fundamental solution of (87) is given by
(90)K(x,t)=H(t)2κπtexp(-x24κt),
where H(t) 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 ofu at the field point z=(x,t) can be expressed as a linear combination of the fundamental solutions U(z,sj):
(91)u(z)=∑j=1ncjU(z,sj),sj=(ηj,τj)∈Ωc,
where n is the number of source points, cj are unknown coefficients, and sj are source points being located in the complement Ωc of Ω=[0,ℓ]×[0,T]. For the heat conduction equation we have the basis functions:
(92)U(z,sj)=K(x-ηj,t-τj).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 thet-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)Bx=b1,
where
(94)Bij=U(zi,sj),x=(c1,…,cn)T,b1=(uℓ(ti),i=1,…,m1;uT(xj),j=1,…,m2;u0(tk),k=m1,…,1(uℓ(ti),i=1,…,m1;uT(xj),)T,
and n=2m1+m2.Since the BHCP is highly ill-posed, the ill-condition of the coefficient matrixB 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:
(95)u(x,t)=cos(πx)exp(-π2t).
For the case with T=1 the value of final time data is in the order of 10-4, which is small by comparing with the value of the initial temperature f(x)=u0(x)=cos(πx) to be retrieved, which is O(1). First we impose a relative random noise with an intensity σ=10% being imposed on the final time data. Under the following parameters m1=11, m2=6, m=10, γ=0.25, and ɛ=ɛ1=10-5, 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-4, which is better than 1.247×10-3 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 m=10 to solve the BHCP, of which both algorithms cannot converge to satisfy the convergence criterion with ɛ=10-5, 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.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.
(a)
(b)
(c)
(d)


## 5.5. Example 5

Let us consider the inverse Cauchy problem for the Laplace equation:(96)Δu=urr+1rur+1r2uθθ=0,(97)u(ρ,θ)=h(θ),0≤θ≤π,(98)un(ρ,θ)=g(θ),0≤θ≤π,
where h(θ) and g(θ) are given function. The inverse Cauchy problem is specified as follows.To seek an unknown boundary functionf(θ) on the part Γ2:={(r,θ)∣r=ρ(θ),π<θ<2π} of the boundary under (96)–(98) with the overspecified data being given on Γ1≔{(r,θ)∣r=ρ(θ),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 ofu at the field point x=(rcosθ,rsinθ) can be expressed as a linear combination of fundamental solutions U(x,sj):
(99)u(x)=∑j=1ncjU(x,sj),sj∈Ωc.
For the Laplace equation (96) we have the fundamental solutions:
(100)U(x,sj)=lnrj,rj=∥x-sj∥.In the practical application of MFS, by imposing the boundary conditions (97) and (98) at N points on (99) we can obtain a linear equations system:
(101)Bc=b1,
where
(102)xi=(xi1,xi2)=(ρ(θi)cosθi,ρ(θi)sinθi),sj=(sj1,sj2)=(R(θj)cosθj,R(θj)sinθj),Bij=ln∥xi-sj∥,ifiisodd,Bij=η(θi)∥xi-sj∥2(ρ(θi)-sj1cosθi-sj2sinθi-ρ′(θi)ρ(θi)cccccccccccccc×[sj1sinθi-sj2cosθi](ρ(θi)-sj1cosθi-sj2sinθi-ρ′(θi)ρ(θi)),ifiiseven,c=(c1,…,cn)T,b1=(h(θ1),g(θ1),…,h(θN),g(θN))T,
in which n=2N, and
(103)η(θ)=ρ(θ)ρ2(θ)+[ρ′(θ)]2.
The above R(θ)=ρ(θ)+D with an offset D can be used to locate the source points along a contour with a radius R(θ). 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:(104)u(x,y)=cosxcoshy+sinxsinhy,
defined in a domain with a complex amoeba-like irregular shape as a boundary:
(105)ρ(θ)=exp(sinθ)sin2(2θ)+exp(cosθ)cos2(2θ).
After imposing the boundary conditions (97) and (98) at N points on (99) we can obtain a linear equations system. The noise being imposed on the measured data h and g is σ=0.01.We solve this problem by the OGSDA withm=10 and γ=0.2. Through 15 iterations the residual is smaller than the convergence criterion ɛ=10-5 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.Example 5 showing (a) residual and steplength of OGSDA, (b) comparing numerical and exact solutions, and (c) numerical errors of OGSDA, GMRES, and BBM.
(a)
(b)
(c)We also apply the GMRES withm=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 functionH(x) for a one-dimensional heat conduction equation:
(106)ut(x,t)=uxx(x,t)+H(x),0<x<ℓ,0<t<tf,(107)u(0,t)=u0(t),u(ℓ,t)=uℓ(t),(108)u(x,0)=f(x).
In order to identify H(x) we impose a Neumann type boundary condition:
(109)ux(0,t)=q(t).We propose a numerical differential method by lettingv=ut. Taking the differentials of (106), (107), and (109) with respect to t and letting v=ut we can derive
(110)vt(x,t)=vxx(x,t),0<x<ℓ,0<t<tf,(111)v(0,t)=u˙0(t),(112)v(ℓ,t)=u˙ℓ(t),(113)vx(0,t)=q˙(t).
This is an inverse heat conduction problem (IHCP) for v(x,t) without using the initial condition.Therefore as being a numerical method, we can first solve the above IHCP forv(x,t) 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 u(x,t) by
(114)u(x,t)=∫0tv(x,ξ)dξ+f(x),
which automatically satisfies the initial condition in (108).From (114) it follows that
(115)uxx(x,t)=∫0tvxx(x,ξ)dξ+f′′(x),
which together with ut=v being inserted into (106) leads to
(116)v(x,t)=∫0tvxx(x,ξ)dξ+f′′(x)+H(x).
Inserting (110) for vxx=vt into the above equation and integrating it we can derive the following equation to recover H(x):
(117)H(x)=v(x,0)-f′′(x).
This approach exhibits threefold ill-posedness: one is the use of the IHCP to solve v(x,t) which is used to provide the data of v(x,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 f(x) in (117).We consider(118)u(x,t)=x2+2xt+sin(2πx),H(x)=2x-2+4π2sin(2πx).
In (117) we disregard the ill-posedness of f′′(x) and suppose that the data f′′(x) are given exactly. We solve this problem by the OGSDA with unit subspace method, where we fix m=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 q˙(t). 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).Example 6 showing (a) steplength of OGSDA, (b) comparing numerical and exact solutions, and (c) numerical error of OGSDA.
(a)
(b)
(c)


## 6. Conclusions

In the present paper, we have derived an optimal algorithm including anm-vector optimal search direction in an m-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 n=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 m of the Krylov subspace dimension is better than other algorithms investigated in this paper.



---

*Source: 154358-2013-12-17.xml*