Galerkin Projection Methods for Solving Multiple Linear Systems.
In this paper, we consider using conjugate gradient (CG) methods for solving multiple linear systems $A^{(i)} where the coefficient matrices $A^{(i)}$ and the right-hand sides $b^{(i)}$ are different in general.\ In particular, we focus on the seed projection method which generates a Krylov subspace from a set of direction vectors obtained by solving one of the systems, called the seed system, by the CG method and then projects the residuals of other systems onto the generated Krylov subspace to get the approximate solutions.\ The whole process is repeated until all the systems are solved.\ Most papers in the literature [T.\ F.\ Chan and W.\ L.\ Wan, {\it SIAM J.\ Sci.\ Comput.}, Peterson, and R.\ Mittra, {\it IEEE Trans.\ Antennas and Propagation}, 37 (1989), pp. 1490--1493] considered only the case where the coefficient matrices $A^{(i)}$ are the same but the right-hand sides are different.\ We extend and analyze the method to solve multiple linear systems with varying coefficient matrices and right-hand sides. A theoretical error bound is given for the approximation obtained from a projection process onto a Krylov subspace generated from solving a previous linear system. Finally, numerical results for multiple linear systems arising from image restorations and recursive least squares computations are reported to illustrate the effectiveness of the method.
Introduction
We want to solve, iteratively using Krylov subspace methods, the following linear systems:
A (i) x
where A (i) are real symmetric positive definite matrices of order n, and in general A (i) 6= A (j)
and b (i) 6= b (j) for i 6= j. Unlike for direct methods, if the coefficient matrices and the right-hand
sides are arbitrary, there is nearly no hope to solve them more efficiently than as s completely
un-related systems. Fortunately, in many practical applications, the coefficient matrices and the
right-hand sides are not arbitrary, and often there is information sharable among the coefficient
matrices and the right-hand sides. Such a situation occurs, for instance, in recursive least squares
computations [20], wave scattering problem [14, 4, 9], numerical methods for integral equations
[14] and image restorations [13]. In this paper, our aim is to propose a methodology to solve
these "related" multiple linear systems efficiently.
In [24], Smith et al. proposed and considered using a seed method for solving linear systems
of the same coefficient matrix but different right-hand sides, i.e.,
In the seed method, we select one seed system and solve it by the conjugate gradient method.
Then we perform a Galerkin projection of the residuals onto the Krylov subspace generated
by the seed system to obtain approximate solutions for the unsolved ones. The approximate
solutions are then refined by the conjugate gradient method again. In [24], a very effective
implementation of the Galerkin projection method was developed which uses direction vectors
generated in the conjugate gradient process to perform the projection. In [6], Chan and Wan
observed that the seed method has several nice properties. For instance, the conjugate gradient
method when applied to the successive seed system converges faster than the usual CG process.
Another observation is that if the right-hand sides are closely related, the method automatically
exploits this fact and usually only takes a few restarts to solve all the systems. In [6], a theory
was developed to explain these phenomena. We remark that the seed method can be viewed
as a special implementation of the Galerkin projection method which had been considered and
analyzed earlier for solving linear systems with multiple right-hand sides, see for instance, Parlett
[19], Saad [21], van der Vorst [26], Padrakakis et al. [18], Simoncini and Gallopoulos [22, 23]. A
very different approach based on the Lanczos method with multiple starting vectors have been
recently proposed by Freund and Malhotra [9].
In this paper, we extend the seed method to solve the multiple linear systems (1.1), with
different coefficient matrices different right-hand sides (b (j) 6= b (k) ). We
analyze the seed method and extend the theoretical results given in [6]. We will see that the
theoretical error bounds for the approximation obtained from a projection process depends on
the projection of the eigenvector components of the error onto a Krylov subspace generated from
the previous seed system and how different the system is from the previous one.
Unlike in [6], in the general case here where the coefficient matrices A (i) can be different, it
is not possible to derive very precise error bounds since the A (i) 's have different eigenvectors in
general. Fortunately, in many applications, even though the A (i) 's are indeed different, they may
be related to each other in a structured way which allows a more precise error analysis. Such is
the case in the two applications that we study in this paper, namely, image restorations and recursive
least squares (RLS) computations. More precisely, for the image restoration application,
the eigenvectors of the coefficient matrices are the same, while for the RLS computations, the co-efficient
matrices differ by rank-1 or rank-2 matrices. Numerical examples on these applications
are given to illustrate the effectiveness of the projection method. We will see from the numerical
results that the eigenvector components of the right-hand sides are effectively reduced after
the projection process and the number of iterations required for convergence decreases when
we employ the projected solution as initial guess. Moreover, other examples involving more
general coefficient matrices (for instance, that do not have the same eigenvectors or differ by a
low rank matrix), are also given to test the performance of the projection method. We observe
similar behaviour in the numerical results as in image restoration and RLS computations. These
numerical results demonstrate that the projection method is effective.
The paper is organized as follows. In x2, we first describe and analyze the seed projection
algorithm for general multiple linear systems. In x3, we study multiple linear systems arising
from image restoration and RLS applications. Numerical examples are given in x4 and concluding
remarks are given in x5.
2 Derivation of the Algorithm
Conjugate gradient methods can be seen as iterative solution methods to solve a linear system
of equations by minimizing an associated quadratic functional. For simplicity, we let
be the associated quadratic functional of the linear system A (i) x . The minimizer of f j
is the solution of the linear system A (i) x . The idea of the projection method is that
for each restart, a seed system A (k) x is selected from the unsolved ones which are then
solved by the conjugate gradient method. An approximate solution - x (j) of the non-seed system
A (j) x can be obtained by using search direction p k
i generated from the ith iteration
of the seed system. More precisely, given the ith iterate x j
i of the non-seed system and the
direction vector p k
i , the approximate solution -
x (j) is found by solving the following minimization
problem:
It is easy to check that the minimizer of (2.2) is attained at -
(p k
and r j
After the seed system A (k) x is solved to the desired accuracy, a new seed system is
selected and the whole procedure is repeated. In the following discussion, we call this method
Projection Method I. We note from (2.3) that the matrix-vector multiplication A (j) p k
j is required
for each projection of the non-seed iteration. In general, the cost of the method will be expensive
in the general case where the matrices A (j) and A (k) are different. However, in x3, we will consider
two specific applications where the matrices A (k) and A (j) are structurally related. Therefore,
the matrix-vector products A (j) p k
j can be computed cheaply by using the matrix-vector product
A
j generated from the seed iteration.
In order to reduce the extra cost in Projection Method I in the general case, we propose
using the modified quadratic function ~
~
to compute the approximate solution of the non-seed system. Note that we have used A (k)
instead of A (j) in the above definition. In this case, we determine the next iterate of the non-
seed system by solving the following minimization problem:
min ff
~
The approximate solution -
x (j) of the non-seed system A (j) x is given by
where
(p k
and ~ r j
Now the projection process does not require the matrix-vector product involving the coefficient
matrix A (j) of the non-seed system. Therefore, the method does not increase the dominant cost
(matrix-vector multiplies) of each conjugate gradient iteration. In fact, the extra cost is just
one inner product, two vector additions, two scalar-vector multiplications and one division. We
call this method Projection Method II. Of course, unless A (j) is close to A (k) in some sense, we
do not expect this method to work well because ~
f j is then far from the current f j .
To summarize the above methods, Table 1 lists the algorithms of Projection Methods I
and II. We remark that Krylov subspace methods (for instance conjugate gradient), especially
when combined with preconditioning, are known to be powerful methods for the solution of
linear systems [10]. We can incorporate the preconditioning strategy into the projection method
to speed up its convergence rate. The idea of our approach is to precondition the seed system
A preconditioner C (k) for each restart. Meanwhile, an approximate
solution of the non-seed system A (j) x also obtained from the space of direction
vectors generated by the conjugate gradient iterations of preconditioned seed system. We can
formulate the preconditioned projection method directly produces vectors that approximate the
desired solutions of the non-seed systems. Table 2 lists the preconditioned versions of Projection
Methods I and II.
all the systems are solved
Select the kth system as seed
for iteration
for unsolved systems
if j=k then perform usual CG steps
oe k;k
x k;k
r k;k
else perform Galerkin projection
x k;j
r k;j
A (j) p k;k
end for
end for
end for
all the systems are solved
Select the kth system as seed
for iteration
for unsolved systems
if j=k then perform usual CG steps
oe k;k
x k;k
r k;k
else perform Galerkin projection
x k;j
r k;j
A
end for
end for
end for

Table

1: Projection Methods I (left) and II (right). The kth system is the seed for the
restart. The first and the second superscripts is used to denote the kth restart and the jth
system. The subscripts is used to denote the ith step of the CG method.
all the systems are solved
Select the kth system as seed
for iteration
for unsolved systems
if j=k then perform usual CG steps
oe k;k
x k;k
+oe k;k
r k;k
z k;k
preconditioning
else perform Galerkin projection
x k;j
r k;j
end for
end for
end for
all the systems are solved
Select the kth system as seed
for iteration
for unsolved systems
if j=k then perform usual CG steps
oe k;k
x k;k
+oe k;k
r k;k
z k;k
preconditioning
else perform Galerkin projection
x k;j
r k;j
end for
end for
end for

Table

2: Preconditioned Projection Methods I (left) and II (right)
We emphasize that in [6, 19, 21, 23, 24], the authors only considered using the projection
method for solving linear systems with the same coefficient matrix but different right-hand sides.
In this paper, we use Projection Methods I and II to solve linear systems with different coefficient
matrices and right-hand sides . An important question regarding the approximation obtained
from the above process is its accuracy. For Projection Method I, it is not easy to derive error
bounds since the direction vectors generated for the seed system A (k) x are only A (k) -
orthogonal but are not A (j) -orthogonal in general. In the following discussion, we only analyze
Projection Method II. However, the numerical results in x4 shows that Projection method I is
very efficient for some applications and is generally faster convergent than Projection Method
II.
2.1 Analysis of Projection Method II
For Projection Method II, we have the following Lemma in exact arithmetic.
Assume that a seed system A (k) x has been selected. Using Projection
Method II, the approximate solution of the non-seed system A (j) x (j) at the ith iteration is
given by
x k;j
where x k;j
' is 'th iterate of the non-seed system, V k
i is the Lanczos vectors generated by i steps
of the Lanczos algorithm if the seed system A (k) x solved by the Lanczos algorithm,
Proof: Let the columns of V k
be the orthonormal vectors of the i-dimensional
Krylov subspace generated by i steps of the Lanczos method. Then we have the following
well-known three-term recurrence
where e i is the ith column of the identity matrix and fi k
i+1 is a scalar. From (2.4) (or see [24]), the
approximate solution x k;j
i of the non-seed system is computed in the subspace generated by the
direction vectors fp k;k
generated from the seed iteration. However, this subspace generated by
the direction vectors is exactly the subspace spanned by the columns of V k
Therefore,
we have
x k;j
Moreover, it is easy to check from (2.4) and (2.5) that
(p k;k
It follows that the solution x k;j
i can be obtained by the Galerkin projection onto the Krylov
subspace K (k) generated by the seed system. Equivalently, x k;j
i can be determined by solving
the following problem:
Noting that the solution is
0 ). the result follows.
To analyze the error bound of Projection Method II, without loss of generality, consider only
two symmetric positive definite n-by-n linear systems:
A (1) x
The eigenvalues and normalized eigenvectors of A (i) are denoted by - (i)
k and q (i)
k respectively and
2. The theorem below gives error bounds for Projection
Method II for solving multiple linear systems with different coefficient matrices and right-hand
sides.
Theorem 1 Suppose the first linear system A (1) x is solved to the desired accuracy in
steps. Let x 1;2
0 be the solution of the second system A (2) x obtained from the
projection onto Km generated by the first system, with zero vector as the initial guess of the
second system (x 0;2
the eigen-decomposition of x
0 be expressed as
Then the eigenvector components c k can be bounded by:
where
6 (q (2)
Here Vm is the orthonormal vectors of Km , P ?
A (1) is the A (1) -orthogonal
projection onto Km and
m A (1) Vm is the matrix representation of the projection of A (1)
onto Km .
Proof: By (2.6), we get x 1;2
x
Since Vm is the orthogonal vectors of Km and
we have kVm
It follows that
6 (q (2)
6 (q (2)
Theorem 1 basically states that the size of the eigenvector component c k is bounded by E k
and F . If the Krylov subspace Km generated by the seed system contains the eigenvectors q (2)
well, then the projection process will kill off the eigenvector components of the initial error of
the non-seed system, i.e., E k is very small. On the other hand, F depends essentially on how
different the system A (2) x (2) is from the previous one A (1) x In particular, when
is small, then F is also small.
We remark that when A A (2) and b (1) 6= b (2) , the term F becomes zero, and as q (1)
the
6 (q (1)
Km )k. It is well-known that the Krylov subspace Km
generated by the seed system contains the eigenvectors q (1)
k well. In particular, Chan and Wan
[6] have the following result about the estimate of the bound sin
6 (q (1)
6
\Gamma- (1)
(- (1)
is the Chebyshev polynomial of degree j. Then
sin
6 (q (1)
If we assume that the eigenvalues of A (1) are distinct, then Tm\Gammak (1+2- k ) grows exponentially
as m increases and therefore the magnitude sin
6 (q (1)
very small for sufficiently large m.
It implies that the magnitude E k is very small when m is sufficiently large. Unfortunately, we
cannot have this result in the general case since q (1)
k , except in some special cases that
will be discussed in the next section.
3 Applications of Galerkin Projection Methods
In this section, we consider using the Galerkin projection method for solving multiple linear
systems arising in two particular applications from image restorations and recursive least squares
computations. In these applications, the coefficient matrices differ by a parameterized identity
matrix or a low rank matrix. We note from Theorem 1 that the theoretical error bound of the
projection method depends on E k and F . In general, it is not easy to refine the error bound E k
and F . However, in these cases, the error bound E k and F can be further investigated.
3.1 Tikhonov Regularization in Image Restorations
Image restoration refers to the removal or reduction of degradations (or blur) in an image using
a priori knowledge about the degradation phenomena; see for instance [13]. When the quality
of the images is degraded by blurring and noise, important information remains hidden and
cannot be directly interpreted without numerical processing. In matrix-vector notation, the
linear algebraic form of the image restoration problem for an n-by-n pixel image is given as
follows:
where b, x, and j are n 2 -vectors and A is an n 2 -by-n 2 matrix. Given the observed image b, the
matrix A which represents the degradation, and possibly, the statistics of the noise vector j, the
problem is to compute an approximation to the original signal x.
Because of the ill-conditioning of A, naively solving will lead to extreme instability
with respect to perturbations in b, see [13]. The method of regularization can be used to
achieve stability for these problems [1, 3]. In the classical Tikhonov regularization [12], stability
is attained by introducing a stabilizing operator D (called a regularization operator), which
restricts the set of admissible solutions. Since this causes the regularized solution to be biased,
a scalar -, called a regularization parameter, is introduced to control the degree of bias. More
specifically, the regularized solution is computed as the solution to
min
b#
A
-D
or min
The term kDx(-)k 2
2 is added in order to regularize the solution. Choosing D as a kth order
difference operator matrix forces the solution to have a small kth order derivative. When
the rectangular matrix has full column rank, one can find the solution by solving the normal
equations
The regularization parameter - controls the degree of smoothness (i.e., degree of bias) of
the solution, and is usually small. Choosing - is not a trivial problem. In some cases a priori
information about the signal and the degree of perturbations in b can be used to choose - [1],
or generalized cross-validation techniques may also be used, e.g., [3]. If no a priori information
is known, then it may be necessary to solve (3.10) for several values of -. For example, in the
L-curve method discussed in [7], choosing the parameter - requires solving the linear systems
with different values of -. This gives rise to multiple linear systems which can be solved by our
proposed projection methods.
In some applications [13, 5], the regularization operator D can be chosen to be the identity
matrix. Consider for simplicity two linear systems:
2:
In this case, we can employ Projection Method I to solve these multiple linear systems as the
matrix-vector product (- 2 I in the non-seed iteration can be computed cheaply by
adding (- 1 I +A T A)p generated from the seed iteration and (- together. Moreover, we
can further refine the error bound of Projection Method II in Theorem 1. Now assume that m
steps of the conjugate gradient algorithm have been performed to solve the first system. We
note in this case that the eigenvectors of the first and the second linear systems are the same,
i.e., q (1)
k . Therefore, we can bound sin
6 (q (1)
using Lemma 2. We shall prove that if
the Krylov subspace of the first linear system contains the extreme eigenvectors well, the bound
for the convergence rate is effectively the classical conjugate gradient bound but with a reduced
condition number.
Theorem 2 Let x 1;2
0 be the solution of the second system obtained from the projection onto Km
generated by the first system. The bound for the A (2) -norm of the error vector after i steps of
the conjugate gradient process is given by
A (2) - 4kx
x 1;2
A (2)
x 1;2
i is ith iterate of the CG process for A (2) x with the projection of x
span fq (1)
'+1 is the reduced condition number of
A (2) and
- (2)
Proof: We first expand the eigen-components of x
It is well-known [10] that there exists a polynomial -
of degree at most i and constant term
1 such that
x 1;2
x 1;2
By using properties of the conjugate gradient iteration given in [10], we have
A
A (2)
A (2)
A (2)
x 1;2
Now the term kx
A (2) can be bounded by the classical CG error estimate,
x 1;2
A (2) - 4kx
x 1;2
A (2)
Noting that
using
Theorem 1 and Lemma 2, the result follows by substitution (2.7) into (3.12).
We see that the perturbation term ffi contains two parts. One depends on the ratio - 2 =- 1
of the regularization parameters between two linear systems and the other depends on how well
the Krylov subspace of the seed system contains the extreme eigenvectors. We remark that the
regularization parameter - in practice is always greater than 0 in image restoration applications
because of the ill-conditioning of A. In particular, - 1 6= 0. If the ratio - 2 =- 1 is near to 1, then
the magnitude of this term will be near to zero. On the other hand, according to Lemma 2, the
Galerkin projection will kill off the extreme eigenvector components and therefore the quantity
in (3.11) will be also small for k close to 1. Hence the perturbation term ffi becomes
very small and the CG method, when applied to solve the non-seed system, converges faster
than the usual CG process.
3.2 Recursive Least Squares Computations in Signal Processing
Recursive least squares (RLS) computations are used extensively in many signal processing and
control applications; see Alexander [2]. The standard linear least squares problem can be posed
as follows: Given a real p-by-n matrix X with full column rank n (so that X T X is symmetric
positive definite) and a p-vector b, find the n-vector w that solves
min w
In RLS computations, it is required to recalculate w when observations (i.e., equations) are
successively added to, or deleted from, the problem (3.13). For instance, in many applications
information arrives continuously and must be incorporated into the solution w. This is called
updating. It is sometimes important to delete old observations and have their effect removed
from w. This is called downdating and is associated with a sliding data window. Alternatively,
an exponential forgetting factor fi, with instance [2]), may be incorporated
into the updating computations to exponentially decay the effect of the old data over time. The
use of fi is associated with an exponentially-weighted data window.
3.2.1 Rank-1 Updating and Downdating Sliding Window RLS
At the time step t, the data matrix and the desired response vector are given by
d t\Gammap+1
respectively, where p is the length of sliding window (one always assumes that p - n). We solve
the following least squares problem: min w(t) Now we assume that a row
is added and a row is removed at the step t + 1. The right-hand-side
desired response vector modified in a corresponding fashion. One now seeks to solve
the modified least squares problem min w(t+1) for the updated least
squares estimate vector w(t + 1) at the time step t + 1. We note that its normal equations are
given by
Therefore, the coefficient matrices at the time step t and t differ by a rank-2 matrix.
3.2.2 Exponentially-weighted RLS
For the exponentially-weighted case, the data matrix X(t) and desired response vector d(t) at
the time step t are defined [2] recursively by
and
where fi is the forgetting factor, and x T
. The RLS algorithms recursively solve for the least squares estimator w(t) at time t,
with t - n. The least squares estimator at the time t and t can be found by solving the
corresponding least squares problems and their normal equations are given by
and
respectively. We remark that these two coefficient matrices differ by a rank-1 matrix plus a
scaling.
3.2.3 Multiple Linear Systems in RLS computations
We consider multiple linear systems in RLS computations, i.e., we solve the following least
squares problem successively
where s is an arbitrary block size of RLS computations. The implementation of recursive least
squares estimators have been proposed and used [8]. Their algorithms updates the filter coefficients
by minimizing the average least squares error over a set of data samples. For instance,
the least squares estimates can be computed by modifying the Cholesky factor of the normal
equations with O(n 2 ) operations per adaptive filter input [20]. For our approach, we employ the
Galerkin projection method to solve the multiple linear systems arising from sliding window or
exponentially-weighted RLS computations.
For the sliding window RLS computation with rank-1 updating and downdating, by (3.15),
the multiple linear systems are given by
1st system : X(t) T
\Theta

s
s
s
s
For the exponentially-weighted case, by (3.16), the multiple linear systems are given by
1st system : X(t) T
\Theta

s
s
According to (3.18), the consecutive coefficient matrices only differ by a rank-2 matrix in
the sliding data window case. From (3.19), the consecutive coefficient matrices only differ by a
rank-1 matrix and the scaled coefficient matrix in the exponentially-weighted case. In these RLS
computations, Projection Method I can be used to solve these multiple linear systems as the
matrix-vector product in the non-seed iteration can be computed inexpensively. For instance,
the matrix-vector product for the new system can be computed by
is generated from the seed iteration. The extra cost is some inner
products. We remark that for the other linear systems in (3.18) and (3.19), we need more inner
products because the coefficient matrices X(t) T X(t) and differ by a rank-s
or rank-2s matrices.
We analyze below the error bound given by Projection Method II for the case that the
coefficient matrices differ by a rank-1 matrix, i.e.,
A
where r has unit 2-norm and each component is greater than zero. For the exponentially-
weighted case, we note that
By using the eigenvalue-eigenvector decomposition of A (1) , we obtain
A
with
is a diagonal matrix containing eigenvalues - (1)
i of A (1) and
r. It has been shown in [11] that if - (1)
k for all k, then the eigenvalues - (2)
k of A (2)
can be computed by solving the secular equation
[(q (1)
(- (1)
0:
Moreover, the eigenvectors q (2)
k of A (2) can be calculated by the formula:
q (2)
Theorem 3 Suppose the first linear system A (1) x is solved to the desired accuracy in m
CG steps. Then the eigenvector components c k of the second system are bounded by jc
6 (q (1)
and
(q (1)
(- (1)
(q (1)
where fq (1)
i g is the orthonormal eigenvectors of A (1) and Km is the Krylov subspace generated
for the first system.
Proof: We just note from Theorem 1 that jc k j - j(P ?
By using (3.20), Theorem 1 and Lemma 2, we can analyze the term j(P ?
6 (q (1)
Since jfl i;k j and j sin
6 (q (1)
are less than 1, we have
small and large i
6 (q (1)
remaining i
From Lemma 2, for i close to 1 or n,
6 (q (1)
sufficiently small when m is large.
Moreover, we note that if ae ? 0, then
see [10]. Therefore, if the values (q (1)
are about the same magnitude for each eigenvector q (1)
then the maximum value of jfl i;k j is attained at either may expect that the
second term of the inequality (3.21) is small when k is close to 1 or n. By combining these facts,
we can deduce that E k is also small when k is close to 1 or n. On the other hand, if the scalar ae
is small (i.e., the 2-norm of rank-1 matrix is small), then F is also small. To illustrate the result,
we apply Projection Method II to solve A (1) x
b (1) and b (2) are random vectors with unit 2-norm. Figures 1 and 2
show that some of the extreme eigenvector components of b (2) are killed off by the projection
especially when jaej is small. This property suggests that the projection method is useful to solve
multiple linear systems arising from recursive lease squares computations. Numerical examples
will be given in the next section to illustrate the efficiency of the method.
In this section, we provide experimental results of using Projection Methods I and II to solve
multiple linear systems (1.1). All the experiments are performed in MATLAB with machine
. The stopping criterion is: kr k;j
tol is the tolerance we
used. The first and the second examples are Tikhonov regularization in image restoration and
the recursive least squares estimation, exactly as discussed in x3. The coefficient matrices A (i) 's
have the same eigenvectors in the Example 1. In Example 2, the coefficient matrices A (i) 's differ
component number
log
of
the
component
RHS before projection
component number
log
of
the
component
RHS after projection
(a) (b)
component number
log
of
the
component
RHS before projection
component number
log
of
the
component
RHS after projection
(c) (d)

Figure

1: Size distribution of the components of (a) the original right hand side b (2) , (b) b (2)
after Galerkin projection when ae = 1. Size distribution of the components of (c) the original
right hand side b (2) , (d) b (2) after Galerkin projection when
component number
log
of
the
component
RHS before projection
component number
log
of
the
component
RHS after projection
(a) (b)
component number
log
of
the
component
RHS before projection
component number
log
of
the
component
RHS after projection
(c) (d)

Figure

2: Size distribution of the components of (a) the original right hand side b (2) , (b) b (2)
after Galerkin projection when ae = \Gamma1. Size distribution of the components of (c) the original
right hand side b (2) , (d) b (2) after Galerkin projection when
Linear Systems (1) (2) (3) (4) Total
Starting with Projection Method I 36 37 43
Starting with Projection Method II 36 48 55 76 205
Starting with previous solution 36 54 66 87 243
Starting with random initial guess 38
Starting with Projection Method I 9 9 9 11 38
using preconditioner
Starting with Projection Method II 9 9 11
using preconditioner
Starting with previous solution 9 13 16 23 61
using preconditioner

Table

3: (Example 1) Number of matrix-vector multiplies required for convergence of all the
systems. Regularization parameter
by a rank-1 or rank-2 matrices. We will see that the extremal eigenvector components of the
right-hand sides are effectively reduced after the projection process. Moreover, the number of
iterations required for convergence when we employ the projected solution as initial guess is less
than that required in the usual CG process.
Example We consider a 2-dimensional deconvolution problem arising in
ground-based atmospheric imaging and try to remove the blurring in an image (see Figure
3(a)) resulting from the effects of atmospheric turbulence. The problem consists of a 256-by-
256 image of an ocean reconnaissance satellite observed by a simulated ground-based imaging
system together with a 256-by-256 image of a guide star (Figure 3(b)) observed under similar
circumstances. The data are provided by the Phillips Air Force Laboratory at Kirkland AFB,
NM through Prof. Bob Plemmons at Wake Forest University. We restore the image using the
identity matrix as the regularization operator suggested in [5] and therefore solve the linear
systems (3.10) with different regularization parameters -. We also test the effectiveness of the
preconditioned projection method. The preconditioner we employed here is the block-circulant-
circulant-block matrix proposed in [5].

Table

3 shows the number of matrix-vector multiplies required for the convergence of all
the systems. Using the projection method, we save on number of matrix-vector multiplies
in the iterative process with or without preconditioning. From Table 3, we also see that the
performance of Projection Method I is better than that of Projection Method II. For comparison,
we present the restorations of the images when the regularization parameters are 0.072, 0.036,
and 0.009 in

Figure

3. We see that when the value of - is large, the restored image is very
smooth, while the value of - is small, the noise is amplified in the restored image. By solving
these multiple linear systems successively by projection method, we can select Figure 3(e) that
presents the restored image better than the others.
(a) (b)
(c) (d)

Figure

3: (Example 1) Observed Image (a), guide star image (b), restored images using regularization
parameter
Linear Systems (1) (2) (3) (4) (5) Total
Starting with Projection Method I 45 31 28 25 24 153
Starting with Projection Method II 45 37
Starting with previous solution 45 43 44 42 40 214
(a)
Linear Systems (1) (2) (3) (4) (5) Total
Starting with Projection Method I 68 51 45 36
Starting with Projection Method II 68 55
Starting with previous solution 68 61 59 56 54 308
(b)

Table

4: (Example 2) Number of matrix-vector multiplies required for convergence of all the
systems. (a) Exponentially-weighted RLS computations and (b) Sliding window RLS computation

Example In this example, we test the performance of Projection Methods
I and II in the block (sliding window and exponentially-weighted) RLS computations.
We illustrate the convergence rate of the method by using the adaptive Finite Impulse Response
system identification model, see [15]. The second order autoregressive process
is a white noise process with variance being 1, is used to
construct the data matrix X(t) in x3.2. The reference (unknown) system w(t) is an n-th order
FIR filter. The Gaussian white noise measurement error with variance 0.025 is added into the
desired response d(t) in x3.2. In the tests, the forgetting factor fi is 0.99 and the order n of filter
is 100.
In the case of the exponentially-weighted RLS computations, the consecutive systems differ
by a rank-1 positive definite matrix, whereas in the case of the sliding window computations, the
consecutive systems differ by the sum of a rank-1 positive definite matrix and a rank-1 negative
definite matrix. Table 4 lists the number of matrix-vector multiplies required for the convergence
of all the systems arising from exponentially-weighted and sliding window RLS computations.
We observe that the performance of Projection Method I is better than that of Projection Method
II. The projection method requires less matrix-vector multiplies than that using the previous
solution as an initial guess. We note from Figures 4 and 5 that the eigenvector components of
b (2) are effectively reduced after projection in both cases of exponentially-weighted and sliding
window RLS computations. We see that the decreases of eigenvector components when using
Projection Method I are indeed greater than those when using Projection Method II.
In the next three examples, we consider more general coefficient matrices, i.e., the consecutive
linear systems do not differ by the scaled identity matrix and rank-1 or rank-2 matrices. In these
examples, the matrix-vector products for the non-seed iteration may not be computed cheaply,
component number
log
of
the
component
RHS before Projection
component number
log
of
the
component
RHS using Projection Method
(a) (b)
component number
log
of
the
component
RHS using Projection Method II
component number
log
of
the
component
RHS using the previous solution as initial guess
(c) (d)

Figure

4: (Example 2) Exponentially-weighted RLS computations. Size distribution of the
components of (a) the original right hand side b (2) , (b) b (2) after using Projection Method I, (c)
b (2) after using Projection Method II, (d) b (2) \Gamma A (2) x (1) (using the previous solution as an initial
component number
log
of
the
component
RHS before projection
component number
log
of
the
component
RHS using Projection Method
(a) (b)
component number
log
of
the
component
RHS using Projection Method II
component number
log
of
the
component
RHS using the previous solution as initial guess
(c) (d)

Figure

5: (Example 2) Sliding window RLS computations. Size distribution of the components of
(a) the original right hand side b (2) , (b) b (2) using Projection Method I, (c) b (2) using Projection
we therefore only apply Projection Method II to solve the multiple linear systems. However, the
same phenomena as in Examples 1 and 2 is observed in these three examples as well.
Example 3 In this example, we consider a discrete ill-posed problem, which
is a discretization of a Fredholm integral equation of the first kind
a
The particular integral equation that we shall use is a one dimensional model problem in image
reconstruction [7] where an image is blurred by a known point-spread function. The desired
solution f is given by
while the kernel K is the point spread function of an infinitely long slit given by
ae sin[-(sin s
We use collocation with n (=64) equidistantly spaced points in [\Gamma-=2; -=2] to derive the matrix
A and the exact solution x. Then we compute the exact right-hand sides
perturb it by uncorrelated errors (white noise) normally distributed with zero mean and standard
derivation 10 \Gamma4 . Here we choose a matrix D equal to the second derivative operator (D =
Different regularization parameters - are used to compute the L-curve
(see

Figure

and test the performance of the Projection Method II for solving multiple linear
systems
s:
We emphasize that the consecutive systems do not differ by the scaled identity matrix.

Table

5 shows the number of iterations required for convergence of all 10 systems using
Projection Method II and using the previous solution as initial guess having the same residual
norm. We see that the projection method requires 288 matrix-vector multiplies to solve all
the systems, but the one using the previous solution as initial guess requires 365 matrix-vector
multiplies. In particular, the tenth system can be solved without restarting the conjugate
gradient process after the projection.
Example 4 We consider the integral equation
Z 2-f(t)dt
corresponding to the Dirichlet problem for the Laplace equation in the interior of an ellipse with
semiaxis c - d ? 0. We solve the case where the unique solution and the right-hand side are
given by
Linear Systems (1) (2) (3) (4) (5) (6) (7) (8)
Starting with Projection 79 38 33 25 27 26 23 21 15 1 288
Method II
Starting with previous 79 44 37 34
solution

Table

5: (Example 3) Number of matrix-vector multiplies required for convergence of all the
systems with -
6.4 6.5 6.6 6.7 6.8 6.9 7 7.1 7.2
least squares residual norm
solution
semi-norm

Figure

(Example 3) The Tikhonov L-curve with regularization parameters used in Table 4.
# of matrix-vector multiply
log
of
residual
norm

Figure

7: (Example 4) The convergence behaviour of all the systems (i)
2:3822 and
d)=(c+d). The coefficient matrices A (k) and the right hand sides b
are obtained by discretization of the integral equation (4.22). The size of all systems is 100.
The values of c and d are arbitrary chosen from the intervals [2; 5] and [0; 1] respectively. We
emphasize that in this example, the consecutive discretized systems do not differ by low rank or
small norm matrices.
The convergence behaviour of all the systems is shown in Figure 7. In the plot, each steepest
declining line denotes the convergence of a seed and also for the non-seed in the last restart.
Note that we plot the residual norm against the cost (the number of matrix-vector multiply)
in place of the iteration number so that we may compare the efficiency of these methods. We
remark that the shape of the plot obtained is similar to those numerical results given in [6] for
the Galerkin projection method for solving linear systems with multiple right hand sides. If
we use the solution of the second system as an initial guess for the third system, the number
of iteration required is 13. However, the number of iteration required is just 8 for Projection
Method II to have the same residual norm as that of the previous solution method; see Figure 8.

Figure

9 shows the components of the corresponding right-hand side of the third system before
the Galerkin projection, after the projection and using the previous solution as initial guess.
The figure clearly reveals that the eigenvector components of b (3) are effectively reduced after
the projection.
Example 5 The matrices for the final set of experiments corresponding to the
three-point centered discretization of the operator \Gamma d
dx (a(x) du
dx ) in [0; 1] where the function a(x)
is given by and d are two parameters. The discretization is performed
using a grid size of h = 1=65, yielding matrices of size 64 with different values of c and d. The
right hand sides of these systems are generated randomly with its 2-norm being 1. We remark
that the consecutive linear systems do not differ by low rank or small norm matrices in this
# of CG iteration
log
of
residual
norm
(a) (b)
(c)

Figure

8: (Example 4) The convergence behaviour of the third system, (a) with projected
solution as initial guesses, (b) with previous solution vector as initial guess and (c) with random
vector as initial guess.
Linear Systems (1) (2) (3) (4) (5) (6) (7) (8)
Starting with Projection 83
Method II
Starting with previous
solution

Table

(Example 5) Number of matrix-vector multiplies required for convergence of all the
systems with c
example.

Table

6 shows the number of iterations required for convergence of all the systems using
Projection Method II and using previous solution as initial guess having the same residual
norm. We observe from the results that the one using the projected solution as the initial
guess converges faster than that using the previous solution as initial guess. Figure 10 shows
the components of the corresponding right-hand side of the seveth system before the Galerkin
projection and after the projection. Again, it illustrates that the projection can reduce the
eigenvector components effectively.
e
components
RHS before projection
(a)
component number
log
of
the
component
RHS after using Projection Method II
component number
log
of
the
component
RHS using the previous solution as initial guess
(b) (c)

Figure

9: (Example distribution of the components of (a) the original right hand side
b (3) , (b) b (3) after Galerkin projection, (c) b (3) \Gamma A (3) x (2) (using the previous solution as an initial
RHS before projection
(a)
component number
log
of
the
component
RHS after using Projection Method II
component number
log
of
the
component
RHS using the previous solution as initial guess
(b) (c)

Figure

10: (Example 5) Size distribution of the components of (a) the original right hand side
b (7) , (b) b (7) after Galerkin projection and (c) b (using the previous solution as an
Concluding Remarks
In this paper, we developed Galerkin projection methods for solving multiple linear systems.
Experimental results show that the method is an efficient method. We end with concluding
remarks about the extensions of the Galerkin projection method.
1. A block generalization of the Galerkin projection method can be employed in many appli-
cations. The method is to select more than one system as seed so that the Krylov subspace
generated by the seed is larger and the initial guess obtained from the Galerkin projection
onto this subspace is expected to be better. One drawback of the block method is that
it may break down when singularity of the matrices occurs arising from the conjugate
gradient process. For details about block Galerkin projection methods, we refer to Chan
and Wan [6].
2. The literature for nonsymmetric systems with multiple right-hand sides is vast. Two
methods that have been proposed are block generalizations of solvers for nonsymmetric
systems; the block biconjugate gradient algorithm [17, 16], block GMRES [25], block QMR
[4, 9]. Recently, Simoncini and Gallopoulos [23] proposed a hybrid method by combining
the Galerkin projection process and Rishardson acceleration technique to speed up the
convergence rate of the conjugate gradient process. In the same spirit, we can modify
the above Galerkin projection algorithms to solve nonsymmetric systems with multiple
coefficient matrices and right-hand sides.



--R

regularization and super- resolution
Springer Verlag

A Block QMR Method for Computing Multiple Simultaneous Solutions to Complex Symmetric Systems
Generalization of Strang's Preconditioner with Applications to Toeplitz Least Squares Problems
Analysis of Projection Methods for Solving Linear Systems with Multiple Right-hand Sides
Analysis of Discrete Ill-posed Problems by Means of the L-curve
Block Implementation of Adaptive Digital Filters
A Block-QMR Algorithm for Non-Hermitian Linear Systems with Multiple Right-Hand Sides
Matrix Computations
Some Modified matrix Eigenvalue Problems
The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind
Fundamentals of Digital Image Processing

Fast RLS Adaptive Filtering by FFT-Based Conjugate Gradient Iterations
Variable Block CG Algorithms for Solving Large Sparse Symmetric Positive Definite Linear Systems on Parallel Computers
The block conjugate gradient algorithm and realted methods
A New Implementation of the Lanczos Method in Linear Problems
A New Look at the Lanczos Algorithm for Solving Symmetric Systems of Linear Equations

On the Lanczos Method for Solving Symmetric Linear Systems with Several Right-Hand Sides
A Memory-conserving Hybrid Method for Solving Linear Systems with Multiple Right-hand Sides
An Iterative Method for Nonsymmetric Systems with Multiple Right-hand Sides
A Conjugate Gradient Algorithm for the Treatment of Multiple Incident Electromagnetic Fields
Etude de quelques m'ethodes de r'esolution de probl'emes lin'eaires de grande taille sur multiprocesseur
An Iteration Solution Method for Solving f(A)
--TR
krylov space;conjugate gradient method;multiple linear systems;galerkin projection