Geoscience Reference
In-Depth Information
Fig. 12.1
The GCR
algorithm for solving
Ax
initialize
x
0
,
;
r
0
=
b
−
Ax
0
;
i
=0;
while
(
r
i
r
i
)
1
/
2
>,
do
i
=
i
+1;
u
i
c
i
D
b
=
r
i−
1
;
=
Au
i
;
for
k
=1
,i−
1
,
do
α
k
=
c
T
c
k
;
=
− α
k
c
k
α
k
u
k
c
i
c
i
;
u
i
u
i
=
−
;
end
;
=
/
(
c
i
c
i
)
1
/
2
;
=
/
(
c
i
c
i
)
1
/
2
;
x
i
=
x
i−
1
r
i
c
i
c
i
u
i
u
i
+
(
c
i
r
i−
1
)
u
i
;
=
r
i−
1
−
(
c
i
r
i−
1
)
c
i
;
end
with the Physical-space Statistical Analysis System (PSAS,
Cohn et al. 1998
),
which employs the conjugate-gradient algorithm applied to (
12.6
)using
R
1=2
as
preconditioner, and it was also displayed in
Zaron
(
2006
) with a non-preconditioned
solver. The non-monotonic reduction in the value of the objective function makes
it problematic to establish an acceptable stopping criteria for the iterative solver. In
spite of the fact that
m<<n
, data sets are frequently large enough that executing
full set of
iterations, the worst-case iteration count for conjugate-gradient-type
linear solvers in exact arithmetic, is prohibitive.
Another issue which arises in practice is that the huge condition number of the
covariance matrices and asymmetry of the linearized model and its approximate
adjoint may cause
R
C
HBH
T
to be non-positive-definite symmetric. Experience
with idealized problems, where the operators can be explicitly constructed as
matrices, shows that the lack of monotonic convergence discussed in the previous
paragraph is exacerbated by symmetry errors and lack of positive-definiteness in the
HBH
T
matrix.
A final consideration in the development of new solvers is the availability of
diagnostic data to assess the progress of the iteration or to evaluate the quality of the
state variable which is obtained.
Recent experience has shown that the generalized conjugate residual (GCR)
method (
de Sturler 1994
,
1996
) addresses all the above-mentioned points. GCR
is a general-purpose Krylov method for solving non-symmetric systems,
Ax
m
D
b
,
R
p
m
such that
AU
D
C
. The columns of both
U
and
C
are in the span of the Krylov subspace
which builds matrices
U
and
C
in
A
p
1
b
g,and
C
is orthogonal, such that
C
T
C
D
I
. The GCR algorithm shown in Fig.
12.1
computes
x
p
2
K
K
D
Span
f
b
;
Ab
;:::;
to minimize k
Ax
p
b
k
2
, which is similar to the
minimum residual
algorithm suggested by
El Akkraoui and Gauthier
(
2010
). Although the GCR
algorithm can fail when either the residual is orthogonal to the Krylov subspace or
when
b
is an eigenvector of
A
p
, neither of these situations has occurred in practice.
Search WWH ::
Custom Search