Dynamical systems gradient method for solving ill-conditioned linear algebraic systems

0.36 MB·English

by  N. S. Hoang

Dynamical systems gradient method for solving
ill-conditioned linear algebraic systems
9
0 N. S. Hoang ∗ A. G. Ramm ‡
† †
0
2
Mathematics Department, Kansas State University,

n †
Manhattan, KS 66506-2602, USA
a
J

8
2
Abstract

] A version of the Dynamical Systems Method (DSM) for solving ill-conditioned linear
A
algebraicsystems is studied in this paper. An a priori anda posteriori stopping rules
N
are justified. An algorithmfor computing the solutionusingaspectraldecomposition
.
h of the left-handside matrix is proposed. Numericalresults showthat when a spectral
t decompositonoftheleft-handsidematrixisavailableornotcomputationallyexpensive
a
m to obtain the new method can be considered as an alternative to the Variational
Regularization.
[
Keywords. Ill-conditionedlinearalgebraicsystems,DynamicalSystemsMethod

3 (DSM), Variational Regularization
v MSC: 65F10; 65F22
3
3
9
1 Introduction
3
.
2
The Dynamical Systems Method (DSM) was systematically introduced and investigated
0
8 in [11] as a general method for solving operator equations, linear and nonlinear, especially
0
ill-posed operator equations. In several recent publications various versions of the DSM,
:
v proposedin [11], were shownto beas efficient andeconomical as variational regularization
i
X methods. This was demonstrated, for example, for the problems of solving ill-conditioned
r linear algebraic systems [2], and stable numerical differentiation of noisy data [8], [9], [3].
a
Theaim of thispaperis toformulate aversion of theDSMgradient methodfor solving
ill-posed linear equations and to demonstrate numerical efficiency of this method. There
is a large literature on iterative regularization methods. These methods can be derived
from a suitable version of the DSM by a discretization (see [11]). In the Gauss-Newton-
type version of the DSM one has to invert some linear operator, which is an expensive
procedure. Thesame is truefor regularized Newton-type versions of the DSMand of their
iterative counterparts. In contrast, the DSM gradient method we study in this paper does
not require inversion of operators.
∗Email: [email protected]
‡Corresponding author. Email: [email protected]
1
We want to solve equation
Au= f, (1)
where A is a linear bounded operator in a Hilbert space H. We assume that (1) has a
solution, possibly nonunique, and denote by y the unique minimal-norm solution to (1),
y := (A) := u : Au = 0 , Ay = f. We assume that the range of A, R(A),
⊥ N N { }
is not closed, so problem (1) is ill-posed. Let f , f f δ, be the noisy data. We
δ δ
k − k ≤
want to construct a stable approximation of y, given δ,f ,A . There are many methods
δ
{ }
for doing this, see, e.g., [4], [5], [6], [11], [13], to mention a few books, where variational
regularization, quasisolutions, quasiinversion, iterative regularization, and the DSM are
studied.
The DSM version we study in this paper consists of solving the Cauchy problem
du
u˙(t) = A∗(Au(t) f), u(0) = u , u N, u˙ := , (2)
0 0
− − ⊥ dt
whereA∗ istheadjointtooperatorA,andprovingtheexistenceofthelimitlim u(t) =
t→∞
u( ), and the relation u( ) = y, i.e.,
∞ ∞
lim u(t) y = 0. (3)
t→∞k − k
If the noisy data f are given, then we solve the problem
δ
u˙ (t)= A∗(Au (t) f ), u (0) = u , (4)
δ δ δ δ 0
− −
and prove that, for a suitable stopping time t , and u := u (t ), one has
δ δ δ δ
lim u y = 0. (5)
δ
δ→0k − k
In Section 2 these results are formulated precisely and recipes for choosing t are
δ
proposed.
Thenovel results inthis paperincludetheproofof thediscrepancy principle(Theorem
3), an efficient method for computing u (t ) (Section 3), and an a priori stopping rule
δ δ
(Theorem 2).
Our presentation is essentially self-contained.
2 Results
Suppose A : H H is a linear bounded operator in a Hilbert space H. Assume that
→
equation
Au= f (6)
has a solution not necessarily unique. Denote by y the unique minimal-norm solution i.e.,
y := (A). Consider the following Dynamical Systems Method (DSM)
⊥ N N
u˙ = A∗(Au f),
− − (7)
u(0) = u ,
0
2
where u is arbitrary. Denote T := A∗A, Q := AA∗. The unique solution to (7) is
0
⊥ N
t
u(t) = e−tTu +e−tT esTdsA∗f.
0
Z
0
Let us show that any ill-posed linear equation (6) with exact data can be solved by the
DSM.
2.1 Exact data
Theorem 1 Suppose u . Then problem (7) has a unique solution defined on [0, ),
0
⊥ N ∞
and u( ) = y, where u( ) = lim u(t).
t→∞
∞ ∞
Proof. Denote w := u(t) y, w = w(0). Note that w . One has
0 0
− ⊥ N
w˙ = Tw, T = A∗A. (8)
−
The unique solution to (8) is w = e−tTw . Thus,
0
kTk
w 2 = e−2tλd E w ,w .
λ 0 0
k k Z h i
0
where u,v is the inner product in H, and E is the resolution of the identity of the
λ
h i
selfadjoint operator T. Thus,
kTk
w( ) 2 = lim e−2tλd E w ,w = P w 2 = 0,
λ 0 0 N 0
k ∞ k t→∞Z h i k k
0
where P = E E is the orthogonal projector onto . Theorem 1 is proved. 2
N 0 −0
− N
2.2 Noisy data f
δ
Let us solve stably equation (6) assuming that f is not known, but f , the noisy data, are
δ
known, where f f δ. Consider the following DSM
δ
k − k≤
u˙ = A∗(Au f ), u (0) = u .
δ δ δ δ 0
− −
Denote
w := u y, T := A∗A, w (0) = w := u y ⊥.
δ δ δ 0 0
− − ∈N
Let us prove the following result:
Theorem 2 If lim t = , lim t δ = 0, and w , then
δ→0 δ δ→0 δ 0
∞ ⊥ N
lim w (t ) = 0.
δ δ
δ→0k k
3
Proof. One has
w˙ = Tw +η , η = A∗(f f), η A δ. (9)
δ δ δ δ δ δ
− − k k ≤ k k
The unique solution of equation (9) is
t
w (t) = e−tTw (0)+ e−(t−s)Tη ds.
δ δ δ
Z
0
Let us show that lim w (t) = 0. One has
t→∞ δ
k k
t
lim w (t) lim e−tTw (0) + lim e−(t−s)Tη ds . (10)
δ δ δ
t→∞k k ≤ t→∞k k t→∞(cid:13)Z (cid:13)
(cid:13) 0 (cid:13)
(cid:13) (cid:13)
One uses the spectral theorem and gets: (cid:13) (cid:13)
t t kTk
e−(t−s)Tdsη = dE η e−(t−s)λds
δ λ δ
Z Z Z
0 0 0 (11)
kTk etλ 1 kTk 1 e−tλ
= e−tλ − dE η = − dE η .
λ δ λ δ
Z λ Z λ
0 0
Note that
1 e−tλ
0 − t, λ> 0,t 0, (12)
≤ λ ≤ ∀ ≥
since 1 x e−x for x 0. From (11) and (12), one obtains
− ≤ ≥
t 2 kTk 1 e−tλ
e−(t−s)Tdsη = − 2d E η ,η
δ λ δ δ
(cid:13)Z (cid:13) Z λ h i
(cid:13) 0 (cid:13) 0 (cid:12) (cid:12)
(cid:13) (cid:13) kT(cid:12)k (cid:12) (13)
(cid:13) (cid:13) t2 d E η ,η
λ δ δ
≤ Z h i
= t2 η 2.
δ
k k
Since η A δ, from (10) and (13), one gets
δ
k k ≤ k k
lim w (t ) lim e−tδTw (0) +t δ A = 0.
δ δ δ δ
δ→0k k ≤ δ→0(cid:18)k k k k(cid:19)
Here we have used the relation:
lim e−tδTw (0) = P w = 0,
δ N 0
δ→0k k k k
and the last equality holds because w ⊥. Theorem 2 is proved. 2
0
∈N
From Theorem 2, it follows that the relation t = C, γ = const, γ (0,1) and C > 0
δ δγ ∈
is a constant, can be used as an a priori stopping rule, i.e., for such t one has
δ
lim u (t ) y = 0. (14)
δ δ
δ→0k − k
4
2.3 Discrepancy principle
Let us consider equation (6) with noisy data f , and a DSM of the form
δ
u˙ = A∗Au +A∗f , u (0) = u . (15)
δ δ δ δ 0
−
forsolvingthisequation. Equation(15)hasbeenusedinSection2.2. Recallthaty denotes
the minimal-norm solution of equation (6).
Theorem 3 Assume that Au f > Cδ. The solution t to the equation
0 δ δ
k − k
h(t) := Au (t) f = Cδ, C = const, C (1,2), (16)
δ δ
k − k ∈
does exist, is unique, and
lim u (t ) y = 0. (17)
δ δ
δ→0k − k
Proof. Denote
v (t) := Au (t) f , T := A∗A, Q = AA∗, w(t) := u(t) y, w := u y.
δ δ δ 0 0
− − −
One has
d
v (t) 2 = 2Re Au˙ (t),Au (t) f
δ δ δ δ
dtk k h − i
= 2Re A[ A∗(Au (t) f )],Au (t) f (18)
δ δ δ δ
h − − − i
= 2 A∗v (t) 2 0.
δ
− k k ≤
Thus, v (t) is a nonincreasing function. Let us prove that equation (16) has a solution
δ
k k
for C (1,2). Recall the known commutation formulas:
∈
e−sTA∗ = A∗e−sQ, Ae−sT = e−tQA.
Using these formulas and the representation
t
u (t)= e−tTu + e−(t−s)TA∗f ds,
δ 0 δ
Z
0
one gets:
v (t) = Au (t) f
δ δ δ
−
t
= Ae−tTu +A e−(t−s)TA∗f ds f
0 δ δ
Z −
0
t
= e−tQAu +e−tQ esQdsQf f
0 δ δ
Z −
0
= e−tQA(u y)+e−tQf +e−tQ(etQ I)f f
0 δ δ
− − −
= e−tQAw e−tQf +e−tQf.
0 δ
−
5
Note that
lim e−tQAw = lim Ae−tTw = AP w = 0.
0 0 N 0
t→∞ t→∞
Here the continuity of A, and the following relation
kTk
lim e−tTw = lim e−stdE w = (E E )w = P w ,
0 s 0 0 −0 0 N 0
t→∞ t→∞Z −
0
were used. Therefore,
lim v (t) = lim e−tQ(f f ) f f δ, (19)
δ δ δ
t→∞k k t→∞k − k ≤ k − k≤
because e−tQ 1. The function h(t) is continuous on [0, ), h(0) = Au f > Cδ,
0 δ
k k ≤ ∞ k − k
h( ) δ. Thus, equation (16) must have a solution t .
δ
∞ ≤
Let us prove the uniqueness of t . Without loss of generality we can assume that
δ
there exists t > t such that Au (t ) f = Cδ. Since v (t) is nonincreasing and
1 δ δ 1 δ δ
k − k k k
v (t ) = v (t ) , one has
δ δ δ 1
k k k k
v (t) = v (t ) , t [t ,t ].
δ δ δ δ 1
k k k k ∀ ∈
Thus,
d
v (t) 2 = 0, t (t ,t ). (20)
δ δ 1
dtk k ∀ ∈
Using (18) and (20) one obtains
A∗v (t)= A∗(Au (t) f )= 0, t [t ,t ].
δ δ δ δ 1
− ∀ ∈
This and (15) imply
u˙ (t) = 0, t (t ,t ). (21)
δ δ 1
∀ ∈
One has
u˙ (t)= Tu (t)+A∗f
δ δ δ
−
t
= T e−tTu + e−(t−s)TA∗f ds +A∗f
0 δ δ
− (cid:18) Z (cid:19) (22)
0
= Te−tTu (I e−tT)A∗f +A∗f
0 δ δ
− − −
= e−tT(Tu A∗f ).
0 δ
− −
From (22) and (21), one gets Tu A∗f = etTe−tT(Tu A∗f) = 0. Note that the
0 0
− −
operator etT is an isomorphism for any fixed t since T is selfadjoint and bounded. Since
Tu A∗f = 0, by (22) one has u˙ (t) =0, u (t)= u (0), t 0. Consequently,
0 δ δ δ
− ∀ ≥
Cδ < Au (0) f = Au (t ) f = Cδ.
δ δ δ δ δ
k − k k − k
This is a contradiction which proves the uniqueness of t .
δ
6
Let us prove (17). First, we have the following estimate:
Au(t ) f Au(t ) Au (t ) + Au (t ) f + f f
δ δ δ δ δ δ δ δ
k − k≤ k − k k − k k − k
tδ (23)
e−tδQ esQQds f f +Cδ+δ.
δ
≤ (cid:13) Z (cid:13)k − k
(cid:13) 0 (cid:13)
(cid:13) (cid:13)
Let us use the inequality: (cid:13) (cid:13)
tδ
e−tδQ esQQds = I e−tδQ 2,
Z k − k ≤
(cid:13) 0 (cid:13)
(cid:13) (cid:13)
and conclude from (23), that
lim Au(t ) f = 0. (24)
δ
δ→0k − k
Secondly, we claim that
limt = .
δ
δ→0 ∞
Assume the contrary. Then there exist t > 0 and a sequence (t )∞ , t < t , such that
0 δn n=1 δn 0
lim Au(t ) f = 0. (25)
n→∞k δn − k
Analogously to (18), one proves that
d v 2
k k 0,
dt ≤
where v(t) := Au(t) f. Thus, v(t) is nonincreasing. This and (25) imply the relation
− k k
v(t ) = Au(t ) f = 0. Thus,
0 0
k k k − k
0 = v(t )= e−t0QA(u y).
0 0
−
This implies A(u y) = et0Qe−t0QA(u y) = 0, so u y . Since u y ⊥, it
0 0 0 0
− − − ∈ N − ∈ N
follows that u = y. This is a contradiction because
0
Cδ Au f = f f δ, 1< C < 2.
0 δ δ
≤ k − k k − k ≤
Thus, lim t = .
δ→0 δ
∞
Let us continue the proof of (17). Let w (t) := u (t) y. We claim that w (t) is
δ δ δ
− k k
nonincreasing on [0,t ]. One has
δ
d
w (t) 2 = 2Re u˙ (t),u (t) y
δ δ δ
dtk k h − i
= 2Re A∗(Au (t) f ),u (t) y
δ δ δ
h− − − i
= 2Re Au (t) f ,Au (t) f +f Ay
δ δ δ δ δ
− h − − − i
2 Au (t) f Au (t) f f f
δ δ δ δ δ
≤ − k − k(cid:18)k − k−k − k(cid:19)
0.
≤
7
Here we have used the inequalities:
Au (t) f Cδ > f Ay = δ, t [0,t ].
δ δ δ δ
k − k ≥ k − k ∀ ∈
Let ǫ > 0 be arbitrary small. Since lim u(t) = y, there exists t > 0, independent
t→∞ 0
of δ, such that
ǫ
u(t ) y . (26)
0
k − k ≤ 2
Since lim t = , there exists δ such that t > t , δ (0,δ ). Since w (t) is
δ→0 δ 0 δ 0 0 δ
∞ ∀ ∈ k k
nonincreasing on [0,t ] one has
δ
w (t ) w (t ) u (t ) u(t ) + u(t ) y , δ (0,δ ). (27)
δ δ δ 0 δ 0 0 0 0
k k ≤ k k ≤ k − k k − k ∀ ∈
Note that
t0 t0
u (t ) u(t ) = e−t0T esTdsA∗(f f) e−t0T esTdsA∗ δ. (28)
δ 0 0 δ
k − k k Z − k≤ k Z k
0 0
Since e−t0T t0esTdsA∗ is a bounded operator for any fixed t , one concludes from (28)
0 0
that lim Ru (t ) u(t ) = 0. Hence, there exists δ (0,δ ) such that
δ→0 δ 0 0 1 0
k − k ∈
ǫ
u (t ) u(t ) , δ (0,δ ). (29)
δ 0 0 1
k − k ≤ 2 ∀ ∈
From (26)–(29), one obtains
ǫ ǫ
u (t ) y = w (t ) + = ǫ, δ (0,δ ).
δ δ δ δ 1
k − k k k ≤ 2 2 ∀ ∈
This means that lim u (t ) = y. Theorem 3 is proved. 2
δ→0 δ δ
3 Computing u (t )
δ δ
3.1 Systems with known spectral decomposition
OnewaytosolvetheCauchyproblem(15)istouseexplicitEulerorRunge-Kuttamethods
with a constant or adaptive stepsize h. However, stepsize h for solving (15) by explicit
numerical methods is often smaller than 1 and the stopping time t = nh may be large.
δ
Therefore, the computation time, characterized by the number of iterations n, for this
approach may be large. This fact is also reported in [2], where one of the most efficient
numerical methods for solving ordinary differential equations (ODEs), the DOPRI45 (see
[1]), is used for solving a Cauchy problem in a DSM. Indeed, the use of explicit Euler
method leads to a Landweber iteration which is known for slow convergence. Thus, it
may be computationally expensive to compute u (t ) by numerical methods for ODEs.
δ δ
However, when A in (15) is a matrix and a decomposition A= USV∗, where U and V
are unitary matrices and S is a diagonal matrix, is known, it is possible to compute u (t )
δ δ
at a speed comparable to other methods such as the variational regularization (VR) as it
will be shown below.
8
We have
t
u (t)= e−tTu +e−tT esTdsA∗f , T := A∗A. (30)
δ 0 δ
Z
0
Suppose that a decomposition
A= USV∗, (31)
where U and V are unitary matrices and S is a diagonal matrix is known. These matrices
possibly contain complex entries. Thus, T = A∗A = VS¯SV∗ and eT = eVS¯SV∗. Using the
formulaeVS¯SV∗ = VeS¯SV∗,whichisvalidifV isunitaryandS¯S isdiagonal, equation(30)
can be rewritten as
t
u (t) = Ve−tS¯SV∗u +V e(s−t)S¯SdsS¯U∗f . (32)
δ 0 δ
Z
0
Here, the overbar stands for complex conjugation. Choose u = 0. Then
0
t
u (t) = V e(s−t)S¯SdsS¯h , h := U∗f . (33)
δ δ δ δ
Z
0
Let us assume that
A∗f = 0. (34)
δ
6
This is a natural assumption. Indeed, if A∗f = 0, then by the definition of h in (33),
δ δ
relation V∗V = I, and equation (31), one gets
S¯h = S¯U∗f = V∗VS¯U∗f = V∗A∗f = 0. (35)
δ δ δ δ
Equations (35) and (33) imply u (t) 0.
δ
≡
The stopping time t we choose by the following discrepancy principle:
δ
Au (t ) f = tδ e(s−tδ)S¯SdsS¯Sh h = e−tδS¯Sh = Cδ.
δ δ δ δ δ δ
k − k (cid:13)Z − (cid:13) k k
(cid:13) 0 (cid:13)
(cid:13) (cid:13)
where 1 < C < 2. (cid:13) (cid:13)
Let us find t from the equation
δ
φ(t):= ψ(t) Cδ = 0, ψ(t) := e−tS¯Sh . (36)
δ
− k k
The existence and uniqueness of the solution t to equation (36) follow from Theorem 3.
δ
We claim that equation (36) can be solved by using Newton’s iteration (43) for any
initial value t such that φ(t ) > 0.
0 0
Let us prove this claim. It is sufficient to prove that φ(t) is a monotone strictly convex
function. This is proved below.
Without loss of generality, we can assume that h (see (36)) is a vector with real
δ
components. The proof remained essentially the same for h with complex components.
δ
First, we claim that
S¯Sh = 0, and S¯Se−tS¯Sh = 0, (37)
δ δ
6 k k6
p p
9
so ψ(t) > 0.
Indeed, since e−tS¯S is an isomorphism and e−tS¯S commutes with √S¯S one concludes
that √S¯Se−tS¯Sh = 0 iff √S¯Sh = 0. If √S¯Sh = 0 then S¯h = 0, and then, by
δ δ δ δ
k k
equation (35), A∗f = S¯h = 0. This contradicts to the assumption (34).
δ δ
Let us now prove that φ monotonically decays and is strictly convex. Then our claim
will be proved.
One has
d
e−tS¯Sh ,e−tS¯Sh = 2 e−tS¯Sh ,S¯Se−tS¯Sh .
δ δ δ δ
dth i − h i
Thus,
d d e−tS¯Sh 2 e−tS¯Sh ,S¯Se−tS¯Sh
ψ˙(t) = e−tS¯Sh = dtk δk = h δ δi. (38)
dtk δk 2 e−tS¯Sh − e−tS¯Sh
δ δ
k k k k
Equation (38), relation (37), and the fact that e−tS¯Sh ,S¯Se−tS¯Sh = √S¯Se−tS¯Sh 2
δ δ δ
h i k k
imply
ψ˙(t) < 0. (39)
From equation (38) and the definition of ψ in (36), one gets
ψ(t)ψ˙(t) = e−tS¯Sh ,S¯Se−tS¯Sh (40)
δ δ
−h i
Differentiating equation (40) with respect to t, one obtains
ψ(t)ψ¨(t)+ψ˙2(t)= S¯Se−tS¯Sh ,S¯Se−tS¯Sh + e−tS¯Sh ,S¯SS¯Se−tS¯Sh
δ δ δ δ
h i h i
= 2 S¯Se−tS¯Sh 2.
δ
k k
This equation and equation (38) imply
e−tS¯Sh ,S¯Se−tS¯Sh 2
ψ(t)ψ¨(t) =2 S¯Se−tS¯Sh 2 h δ δi S¯Se−tS¯Sh 2 > 0. (41)
k δk − e−tS¯Sh 2 ≥ k δk
δ
k k
Here the inequality: e−tS¯Sh ,S¯Se−tS¯Sh e−tS¯Sh S¯Se−tS¯Sh was used. Since
δ δ δ δ
h i ≤ k kk k
ψ > 0, inequality (41) implies
ψ¨(t) > 0. (42)
It follows from inequalities (39) and (42) that φ(t) is a strictly convex and decreasing
function on (0, ). Therefore, t can be found by Newton’s iterations:
δ
∞
φ(t )
n
t = t
n+1 n− φ˙(t )
n
(43)
= t + ke−tnS¯Shδk−Cδ e−tnS¯Sh , n = 0,1,...,
n S¯Se−tnS¯Shδ,e−tnS¯Shδ k δk
h i
for any initial guess t of t such that φ(t ) > 0. Once t is found, the solution u (t ) is
0 δ 0 δ δ δ
computed by (33).
10