A simple factorization is given of an arbitrary hermitian, positive definite matrix in which the factors are well-conditioned, hermitian, and positive definite. In fact, given knowledge of the extreme eigenvalues of the original matrix A, an optimal improvement can be achieved, making the condition numbers of each of the two factors equal to the square root of the condition number of A. This technique is to applied to the solution of hermitian, positive definite Toeplitz systems. Large linear systems with hermitian, positive definite Toeplitz matrices arise in some signal processing applications. A stable fast algorithm is given for solving these systems that is based on the preconditioned conjugate gradient method. The algorithm exploits Toeplitz structure to reduce the cost of an iteration to O(n log n) by applying the fast Fourier Transform to compute matrix-vector products. Matrix factorization is used as a preconditioner

Pan, Victor

Schrieber, Robert

English

NASA Technical Reports Server

d T_2
A Fast, Preconditioned Conjugate
Gradient Toeplitz Solver
Robert Schrieber
M,_Lrch, 1989
Research Institute for Advanced Computer Science
NASA Ames Research Center
RIACS Technical Report 89.14
NASA Cooperative Agreement Number NCC 2-387
(NA_A-C_-l. GO367) A FAST, PRECONDITIONED
CONJUGAT_ GRAOIFNT TOrPLITZ S_LVER
{Rese]rch Inst. for Advanced Computer
science) 16 p CSCL f_9_
G3/o]_
f_90-2_T3[
unclas
0230o73
Research Institute for Advanced Computer Science
https://ntrs.nasa.gov/search.jsp?R=19900015415 2020-03-19T22:02:00+00:00Z

A Fast, Preconditioned Conjugate
Gradient Toeplitz Solver
Victor Pan
Robert Schrieber
March, 1989
Research Institute for Advanced Computer Science
NASA Ames Research Center
RIACS Technical Report 89.14
NASA Cooperative Agreement Number NCC 2-387

A Fast, Preconditioned Conjugate
Gradient Toeplitz Solver
Victor Pan*
Computer Science Department, SUNY, All)any, NY 12222
ait(t
Mathematics and Computer Sci_mce Department,
Lehman College, CUNY, Bronx, NY 10468
Robert Schreiber _
Research Institute for Advanced Computer Science
Moffett Field, CA 94035
March 15, 1989
Abstract
We give a, simple factorization of an arbitrary hermitian, positive deti,_it,.
matrix in which the factors are well-conditioned, hermitian, and positive
definite. In fact, given knowl(xlge of the exl, rcme eigcnva:lues of the original
matrix A, we can achieve an optimal iml)rovement , making the condition
numbers of each of the two factors eqmd to the square root of the conditiol,
number of A.
_vVe apply !_his technique to the solution of hermitian, positive definite
Toeplitz systems. Large linear systems with hermitian, positive definite
*Supported by NSF Grant CCR-8805782 and PSC-CUNY A_vard 668541
_Supported by ONR Contract number N00014 86-K-0610, ARO Grant DAAL03-86-
K-0112, and by Cooperative Agreement NCC2-a87 between NASA and the Universi',i<:
Space Research Association.
t On leave from l{e.nsselaer Polyt.ech n ic I nstit.u tc
Toeplitz matrices arise in some signal processing applications. We give a
stable fast algorithm for solving these systems that is based on the precondi-
tioned conjugate gradient method. The algorithm exploits Tocplitz structure
to reduce the cost of an iteration to O(n log n) by applying the fasl Fourier
Transform to compute matrix-vector products. We use our matrix factoriza-
tlon as a preconditioner,
1 Introduction
We give a simple factorization of an arbitrary hermitian, positive definite
matrix A in which the factors are hermitian, positive definite, and are sub-
stantially better conditioned than the original matrix A. In fact, given knowl-
edge of the extreme eigenvalues of A, we can achieve an optimal improve-
ment, making tile condition numbers of each of the two factors equal to the
square root of the condition number of A. The factorization is of the form
A = (A + #I)(I - It(A + ltl)-_). We discuss the optimum choice of t_ in
Section 3.
Consider the linear system Ax = b where A is an n x n hermitian, positive
definite Toeplitz matrix, A = [a_j] = [aj_] = [al__Jl 1. Several direct methods
for solving such a system using O(n 2) arithmetic operations are known [18],
[14], [12], [3]. In addition, some newer methods that take O(n log 2 n) opera-
tions have been developed [4], [5], [13], ill], [2]. The method of Gragg and
Ammar, for example, requires 8n log 2 n real arithmetic operations for a real
Toeplitz system. Stability of these methods is discussed by Bunch [6].
Recently, some new attention has b¢_n given to the preconditioned con-
jugate gradient n mthod as a Toeplitz solver. The motivation is that a single
iteration of this method can be implemented, using the fast Fourier Trans-
form, at a cost of O(n log n) arithmetic operations. The hope is that with
suitable preconditioners, the number of conjugate gradient iterations can be
made small enough to make the method practical. One approach has been
to use circulant approximations to A as preeonditioners [16], [7]. These work
remarkably well for some Toeplitz matrices (finite sections of a singly in-
finite matrix Aoo = [a__,] for which Ek=0_ ak < oo) [8], [17]. As we shall
show in Section 4, however, they are not always effective. Here we use in-
stead the general preconditioning strategy discussed above. This strategy is
most useful when A is Toeplitz, for in this case, the factor, (A -4-ttI), is still
2
Toeplitz. Conjugategradient iterations are thereforeinexpensive.Moreover
its inverse has a representationas the differencebetween two products of
triangular Toeplitz matrices (the Gohberg-Semencul formula) or a triangu-
lar Toeplitz matrix and a circulant matrix (the Ammar-Gader formula); this
allows iteration with the other factor, (I - #(A +/tI)-Z) to be carried out
cheaply.
In the next section we give a general outline of our method. In Section 3
we give the details of this method, providing an optimum shift parameter FL.
We compare it with several competitors in Section 4.
1.1 Notation
For z C C,z" denotes the complex conjugate of z. For a complex matrix
A, A H denotes the conjugate transpose of A, and A(A) denotes the set of
eigenvalues of the square matrix A. A matrix A is hermitian, positive definite
(hereafter h.p.d.) ifA = A n and for alla E A(A),a > 0. Let A beh.p.d.
and let al be the largest and a,, the smallest of A's eigenvalues. Then
a(A) = al/an is called the (spectral) condition number of A.
2 A condition-improving matrix factoriza-
tion
Let A be a given h.p.d, matrix. Let p E R. Let
B - A + _I (1)
and
C = I- #B-'. (2)
Lemma 1 Let A be a given h.p.d, matrix. Let B and C be given by ([) aud
(2). Then A = BC = CB. If -t_ is not an eigenvalue of A then both B and
C have inverses and A -1 = C-aB -I = B-1C -a.
Proof:
obvious.
]emma.
A = B - pI = B(I - #B -1) = (1 - #B-_)B. Invertibility of/3 is
And C = B-1A must therefore be invertible too. This proves the
3
Let the eigenvalues of A, B, and C be given by
{_, <....<_ _,} = A(A),
{_. _<... _<_,} _(B),
{7. _<"-< _,} = a(C).
By (1) and (2) we havc
(a)
(4)
Lemma 2 Lct B and C be given by (1) and (2).
orB and C arc given bg
,_(_) -
an + I*
and
Thus, .for all tt > O,
Then the condition numbers
(5)
cq(a,,+ St) (6)
,_(C)- _,(_, + St).
n(A)--- x(B)_(C). (r)
Proof: Immediate.
Lemma 3 Let # = _v/-dTdT.Then _( B) = n(C) = _.
Proof: Immediate.
In most cases, this factorization does not lead to a reduction in the cost of
inverting A or of solving Ax = b. For if we employ an algorithm for inverting
B and C, such as Gaussian elimination, that is insensitive to condition num-
ber, we might as well have used it on A; if the cost of this algorithm depends
on log _, as it does for Newton's method (see, for example, [15]) then we
have gained nothing. If the cost exceeds log n, as it does for most iterative
methods, we may have gained something. But iteration to solve Cx = z is
expensive because we need to solve a system with B at each iteration. Some
form of inner/outer iteration method could be used to reduce this cost.
There is a case,however,in which we cansolvemany systemswith B in
little more time than it takes to solve one. When A is an h.p.d. Toeplitz
matrix, then so is B. Then we may use the Gohberg-Semcncul formula, or
a recent variant, the Ammar-Gader formula, to represent B -1. This reduces
the cost of one iteration of CG for C to O(n log n).
3 A Fast Toeplitz Solver
Let us apply the matrix factorization of the previous section to the solution
of a linear system
Ax = t, (8)
where A = [a_j] is an n x n h.p.d. Toeplitz matrix, aij = aji = alj_i j. We
suppose that the extremal eigenvalues of A have been estimated (more on
this later on, in Section 5). In this case, the matrix B is also h.p.d, and
Toeplitz, is completely defined by its first column, and its inverse can be
represented as
B -1 = LIL_ - L_L_ t
where L1 and L2 are lower triangular Toeplitz matrices whose elements de-
pend only on the first column of B -1 ([9], [18]). Recently, an alternative rep-
resentation by Toeplitz-circulant products was given by .A,mmar and Cadet
[1]:
13 -1 = L1E T - L2E (9)
where Ll and L2 are again lower triangular Toeplitz, and E is circulant.
These factors are again functions only of the first column of B-1.
This leads to the following algorithm:
Algorithm 1:
Input: An h.p.d. Toeplitz matrix A, a vector b, and a shift parameter p.
Output: A-lb.
Method:
(Step 1) Solve By = e,, where el = (1,0,...,0) H. Construct L1, L2, E
satisfying (9).
(Step 2) Solve Bz = b.
(Step 3) Solve Cx = z. Return x.
5
We solve tile linear systems at Steps 1 and 3 of Algorithm 1 by the
conjugate gradient (CG) method [10]. At Step 2, we use the representation
(9) constructed at Step l; there is no iteration and only FFTs of n-vectors
are required.
In the CG method, each iteration costs 5n flops and one matrix-vector
product. For the Toeplitz matrix B in Steps 1 and 2, we compute the matrix
vector product Bu by a standard technique of embedding B in a circulant
matrix of order 2n and appending n zeros to u; a circulant matrix times a
vector is a convolution, for which the fast Fourier Transform is used:
;, = 01);
_5 = _.,/);
,, =
Bu = v(1 :n);
Here 0 represents the zero vector of order n; Fk and Wk are the forward and
inverse discrete Fourier Transform operators on C k, the operator .* represents
elementwise multiplication of two vectors, v(1 : n) is the first n components
of v, and
-_- /7_2n([bo, bll. ,bii'l,O,b;_l,...,bl] ).
Note that/) is real and this reduces the cost of the elementwise multiplica-
tion. The cost of a CG iteration with B is therefore
Cost(B) : 4@(n)
where ¢(n) is tile cost of an n-point FFT. (This does not include the cost of
computing/), which is computed once and for all.)
For Step 3 above, we use CG, too. We require tile product Cu and hence
also B-Xu at every iteration. Using the Ammar-Gader formula, this cost may
be reduced to
Cost(C) = 7(I)(n)
(see [1]).
The alternative of using B as a preconditioner in the preconditioned CG
method is less attractive than Algorithm 1. We would in effect be solving
Cx = B-lb
by the conjugate gradient method. Thus, we have an equivalent of Steps
2 and 3 above. But at each iteration, we would have to form Cp for some
vector p as B-lAp. The cost of each iteration would therefore be ll_(r_)
rather than 7_(n) as it is in our implementation.
3.1 The Optimal Shift
The choice # = x/_a_ is not optimal. Since each iteration of CG with the
matrix C costs roughly 1.75 times as much as a CG step with B, we would
be better off to choose # to make C better conditioned at tile expense of
worsening the condition number of B. It would be reasonable to choose # lo
minimize the estimate of the total work given by
+ 7nc)¢(n),
where nB is the number of CG iterations at Step 1 above, nc tile number of
CG iterations at Step 3 above. According to [10], we may model these by
and
ns = F_/_ (lo)
nc = f @_(C) (11)
with F a constant. Then, by Lemma 2,
n_nc = F2_-_ - M = const.
Therefore, for any constant K (and motivated by the estimate above we may
think of K = 7/4),
f(nc) -- nB + Knc
M
= _ + Knc
rtc
is minimized at
nc
nB
M
nc
(12)
Now wewant to choosethe best shift parmeter #. In view of (10), (11), and
(12) wechoose# to solve
_-(B)= i_%(c), (13)
which becomes, using (5) and (6),
_.(_[ + 2._1 + v')= K=-,(-'. + 2.4. + #9.
Now change variables from t, to m, where # = m _x/-&'7"_'_;then m satisfies
../2( K2t.2¢, __ C_rt ) JF ?Tt(2( K 2 -- 1) _X/'_'_n ) -1-- (I(20grt --O_1) = 0
or, dividing by an > 0, with _ =- n(A) = alia,,,
rn2(I(2g - 1) 4 m(2(K 2- 1)v/-_) + (K 2 - x) =0.
The roots are
-(K = - 1)v/-_ + K(,_- 1) (14)
no,+ = K2 e; - 1
Note that rn+ _> 0only when __> K 2. Ifx < K 2 there is no possibi]ityof
satisfying (13) anyway, in view of (7). Thus, we assume that x >_ 1( 2 and
that # = m+ _x/'ffT"_ > 0. Fol _ > K 2 the root m+ is monotone increasing
with g and is asymptotic to 1/K.
Lemma 4 Let/, be given by rn_ where m = rn+ (see (14) above). Then
t-(U) : K_, (1.5)
_(C) = K-'_. (16)
Proof: According to the derivation of (14), the relation (13) holds. Then,
by (7),
= _(t_).(c)
= (_(R)/I_-)_,
whence (15) follows. And (16) follows from (13) and (15). This proves the
lemma.
Now,assume(10) and (11) and choose/_= m+v'_C_ so that (13) holds.
In our application, K = t.75. The total cost is then
= 4(nB + Knc)_(n)
= 4(2nB)@(n)
=
= 8vqT '/4F (n)
=
(4rib + 7nc)_(n)
(17)
For comparison, let n,co be the number of iterations required by CG for A.
We assume the model
Coat(CG) = 4nca¢(n) = 4F¢(n)v/_A) • (18)
Compared with (18) we can see that the new method is an improvement for
_; > 49.
3.2 Recursiv:e Preconditioning
The factorization A = BC can be used recursively. In particular, let us
consider solving the equation By = el at Step I of Algorithm I by using this
preconditioned conjugate gradient solver. We now solve for two optimum
shift paramemters; g_, which we use in solving for y in Step 1, is given by
the theory above in which we substitute ¢_, 5,,, and _(B) for C_l,C_, and
t_(A). The other,/l, which we use to define B, is now chosen to minimize the
sum of the cost of Step 1, which by (17) is
4v/qn( B) '14 F_(n) ,
and the cost of Step 3, which is
7n(C)'/2F@(n)
subject to (7). The solution is to take
K(C) = (16h:(A)49 ]'_1/3
and
=
Making the substitution It = m _/-d-]'_ as before leads to a cubic in m that
has a positive root for all _¢> 49/16. The overall work becomes
9 (_) 2/3 F_(n)_(A)1/6 .
This is better that 4F(P(n)tc(A) 1/2 (compare (18)) for s:(A) > 140 arid better
than 4v/-']FC,(n)(K(.4)) '/4 (compare (17)) for _(A) > 3222.
4 Numerical Results and Conclusions
First, we used Algorithm 1 to solve the system (8) where A was as follows.
We generated a real time series
s_. = 2sin(x_) + sin(4xk + _r/6) + u_
where tile n sample points xk were uniformly spaced in [0, 2_-) and vk was a
Gaussian random variable (white noise) with mean zero and variance (power)
0"2. Next, we took aj to be the autocorrelation of {sk} with lag j, that is,
n-1
aj _--- _ 5kNk+ j
k=0 = :
where the index k + j was taken modulo n.
To choose/z, we used tt,_ Lanczos algorithm to approximately tridiago-
nalize A, then con_l)uted the extreme eigenvalues ri and rm of the resulting
m x m symmetric tridiagonal matrix. We used the approximations
Tables 1 and 2 give experimental results for a = 10 and a = 1. In each
table we list the number of sample points, n; the number of CG iterations
at Step 1 of Algorithm 1, nR; the number of CG iterations at Step 3 of
10
n16
64
256
1024
4096
nB
f7
39
50
70
59
nc
7
13
18
23
38
7 nnB+_ C
29
62
82
110
126
nCG
17
46
73
12,1
162
Ilax - yli_
1.8(-12)
9.7(-12)
1.2(-11)
5.3(-11)
1.2(-10)
Ilax_c-vll
1.2613)
6.9(-12)
8.1(-12)
2.o(-_1)
6.1(-11)
Table 1: Results for _7 = 10.
rt
f6
64
256
1024
4096
rt B
17
41
49
50
35
nc
7
np+ -_nc nCG IIA_- Vll, llAxc_-,_ll
12
32
47
75
217
38
97
131
i81
415
18
80
156
253
540
2.s(-12)
6.6(-12)
2.s(-_a)
s.7(-_1)
1.7(-9)
1.9(-13)
7.7(-12)
1.a(-t t )
2.7(-_1)
Table 2: Results for cr = I .
11
7t
16
64
256
i024
4096
4096
nB
18
51
64
67
46
109
_c
12
24
34
54
J69
73
7
nB + -_nc
39
93
124
162
342
237
nCG
18
80
156
253
540
540
IIAx- YII 
7.7(-13)
9.0(-12)
2.9(-11)
7.9(-11)
1.9(-10)
4.0(-9)
5
5
5
5
5
10
Table 3: Results for o = 1. Shift # = rm/£
Algorithm 1, nc; the equivalent number of CG iterations taken by Algorithm
T
1, nB+ _nc; the number of iterations taken by unpreconditioned CG, nca;
the residuals achieved by CG and by Algorithm 1. The number of Lanczos
iterations that were used to estimate eigenvalues varied from 10 to 40.
In our experience, TI is extremely accurate, but rm may be several times
larger than (_,_. We also found that an underestimate of a,,_, which results
in a smaller # than we would take given perfect knowIedge and hence a bias
in favor of more iterations aI, Step 1 and fewer at Step 3, is preferable to an
overestimate. Thus, we replaced our estimate of c_,, by
rm/5
where 6 > 1 is a parameter of the algorithm. Results are given in Table 3 for
the same (:lass of matrices as above, with a = 1. Note that with this better
estimate of the optimal shift, we are doing far better than with unprecondi-
tioned CG.
Strang [16] and Chan [7] have advocated circulant preconditioners. Strang
takes the circulant C = [c(j_ Omo_j,,] with ck = ak, k = 0,...,rn where
m = n/2. The other diagonals of C are determined by symmetry and the
requirement that C be circulant. Thus C coincides with A in the central
half of the diagonals. Chan's choice is to take C to minimize IIA - CHF , the
Frobenius norm of the difference, among all circulaj)t C. For this, one takes
ck = (ka,__k + (n - k)ak)/n. The cost of each iteration is 1.5 times greater
than the cost of a CG iteration, since solving the circu]ant preconditioning
system at each step requires both a forward and an inverse FFT of length n.
12
n 1.5 ns 1.5 nc
16 33 30
64 122 114
256 662 293
Table 4: Results for Chan and Strang preconditioners, o" =- 1.
We used these techniques for the model problem above, with c_ = 1. q'he
results, in Table 4, show that neither technique is competitive with ours
for this problem. For other problems, in particular when ak decays with
increasing k, we have found them to be superior, as has been predicted by
Strang, Chan, and Edelman ([17]).
5 Conclusions
We have demonstrated a preconditioner for Toeplitz systems that achieves a
significant speedup when compared with unpreconditioned conjugate gradi-
ent method and, for certain problems, is considerably better than other previ-
ously proposed preconditioners. The asymptotic complexity is O(n log n_ 1/4)
where _; is the condition number of the problem the exponent I/4 may be
made smaller at the expense of an increase in the constant factor by applying
our technique recursively. Compared with the complexity (8n log 2 n) of other
superfast methods, we can see that for large, well-conditioned problems the
present technique may be quite useful.
6 Acknowledgement
We would like to thank Joowan Chun and Tom t(ailath for pointing oul
Reference [1].
References
[1] Gregory Ammar and Paul Gader. A variant of the Gohberg-Semencul
formula involving circulant matrices. Subrnitted to SIAM Journa/ on
Matrix Analysis and Applications.
13
[2] GregoryS.Ammar and William I3.Gragg.Superfastsolution of real pos-
itive definite Toeplitz systems.SIAM J. Matrix Anal. Appl. 9:1 (1988),
pp. 61-76.
[3] E. H. Bareiss.Numerical solution of linear equationswith Toeplitz and
vector Toeplitz matrices.Numer. Math. 13 (1969),pp. 404-424.
[4] R. R. Bitrnead and B. D. O. Anderson.Asymptotically fast solution of
Toeplitz and related systemsof linear equations.Linear Algebra Appl.
34 (1980), pp. 103-116.
[5] R. P. Brent, F. G. Gustavson and D. Y. Y. Yun. Fast solution of Toeplitz
systems of equations and computation of Pad6 approximants. J. Algo-
rithms 1 (1980), pp. 259-295.
[6] James R. Bunch. Stability of methods for solving Toeplitz systems of
equations. SIAM J. Scient. and Star. Comput. 6:2 (1985), pp. 365-375.
[7] Tony F. Chan. An optimal circulant preeonditioner for Toeplitz systems.
SIAM J. Scient. and Star. Comp. 9:4 (1988), pp. 766-771.
[8] Raymond H. Chan and Gilbert Strang. Toeplitz equations by conju-
gate gradients with circulant preconditioner. SIAM J. Scient. and Star.
Comput. 10:1 (1989), pp. 104-119.
[9] I.C. Gohberg and A.A. Semencul. On the Inversion of Finite To_,.plitz
Matrices and Their Continuous Analogs. Mat. Issled. (in Russian) 2:1
(1972), pp. 201-2.33.
[10] G.H. Golub and C.F. van Loan . Matrix Computations. Johns IIopkins
Univ. Press, Baltimore, Maryland, 1983
[11] F. de Hoog. A new algorithm for solving Toeplitz systems of equations.
Linear Algebra Appl. 88/89 (1987), pp. 122-138.
[12] J. R. Jain. An efl'icient algorithm for a large Toeplitz set of linear equa-
tions. IEEE Trans. Acoust. Speech and Signal Processing 27 (1979), pp.
612-615.
14


A fast, preconditioned conjugate gradient Toeplitz solver

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server