AbstractA variant of the preconditioned conjugate gradient method to solve generalized least squares problems is presented. If the problem is min (Ax − b)TW−1(Ax − b) with A ∈ Rm×n and W ∈ Rm×m symmetric and positive definite, the method needs only a preconditioner A1 ∈ Rn×n, but not the inverse of matrix W or of any of its submatrices. Freund's comparison result for regular least squares problems is extended to generalized least squares problems. An error bound is also given

Iusem, A.N.

Yuan, J.Y.

Elsevier - Publisher Connector 

Preconditioned conjugate gradient method for generalized least squares problems 

JOURNAL OF 
COMPUTATIONAL AND 
APPLIED MATHEMATICS 
ELSEVIER Journal of Computational nd Applied Mathematics 71 (1996) 287-297 
Preconditioned conjugate gradient method for generalized least 
squares problems 
J .Y .  Yuan  a,., 1, A .N .  Iusem b,2 
a Departamento de Matemrtiea, Universidade Federal do Paran6, Centro Politkenieo, CP: 19.081, Curitiba, CEP: 
81531-990, Brazil 
b Instituto de Matemhtica Pura e Aplieada, Estrada Dona Castorina 110, Jardim Botfmieo, Rio de Janeiro, CEP: 
22460-320, Brazil 
Received 10 March 1995; revised 27 September 1995 
Abstract 
A variant of the preconditioned conjugate gradient method to solve generalized least squares problems is presented. If
the problem is 
rain (Ax - b)T W - 1 (Ax - b ) 
with A E ~,,,x,, and W E •mx,, symmetric and positive definite, the method needs only a preconditioner A1 E R nxn, but 
not the inverse of matrix W or of any of its submatrices. Freund's comparison result for regular least squares problems 
is extended to generalized least squares problems. An error bound is also given. 
Keywords." Preconditioned conjugate gradient method; Least squares; Generalized least squares problems 
AMS classifications: 65F20 
1. Introduction 
Paige [7, 8] transformed the generalized least squares problem 
min (Ax  - b)  T W+(Ax  - b), 
x E R" 
(1.1) 
* Corresponding author. E-mail: yuan@inf.ufpr.br., jin@mat.ufpr.br, 
1 The work of this author was supported by CNPq grant 301035/93-8. 
2 The work of this author was supported by CNPq grant 301280/86. 
0377-0427/96/$15.00 (~)1996 Elsevier Science B.V. All rights reserved 
SSDI  0377-0427(95)00239-1  
288 J. Y. Yuan, A.N. Iuseml Journal of Computational and Applied Mathematics 71 (1996) 287-297 
where A E Rmx”, WE lRmx” is symmetric and positive semi-definite, and m > II, into an equality- 
constrained least squares problem 
min 
xEn”,&-@ II412, (l-2) 
s.t. Ax + Bu = b, (1.3) 
where W = BBT, B E R” xk and k = rank( W). Then he presented direct method to solve problem 
(1.2) and (1.3). He also provided an error analysis for his method. However, Paige’s method is a 
direct method which is not appropriate to handle large sparse problems. 
In terms of (1.2) and (1.3), problem (1.1) can be solved by implicit null space iterative methods 
given by James in [5, 61, or by the BNP algorithm [2]. But, as in Paige’s method, the BNP algorithm 
needs to decompose the weighted matrix W, and find the inverse of B or of a submatrix of (A, B). 
Also the BNP algorithm does not work when W is semi-definite because (A, B) is not full rank. 
Even when W is positive definite, the BNP algorithm needs to invert an m x m matrix to be used 
as a preconditioner. In order to overcome these limitations, Yuan [9] recently proposed the 2-cyclic 
SOR method for problem (1 . 1 ), which we describe next. 
Algorithm 1.1 (2-cyclic SOR Algorithm). 
(1) Factorize Al and Wz2, set initial vectors X(O) = 0, Y(O) = 0; 
(2) Select a relaxation parameter 0; 
(3) Iterate for k = 0, 1, . . . , until “convergence” 
xCk+‘) = (1 - o)_x(~) + WA,’ [b, - (W,, - W,,PT)$‘], 
#kfl) = (1 - w)yp 
2 + coW2;‘(bz - W,$-ik’ - A2dk+‘)), 
rjk+l) = (1 - o)ri (k) _ opTr(k+') 2 3 
where 
w= 3 
and 
P = A2A;‘. 
This method is always convergent (see [9]), but it requires estimation of the relaxation para- 
meter cc). 
It is well known that the conjugate gradient method does not need any iterative parameter and 
enjoys a finite termination property. In this paper, the conjugate gradient method is extended to 
solve generalized least squares problems by preconditioner techniques. The preconditioned conju- 
gate gradient method for problem (1.1) also has the finite termination property in the absence of 
roundoff errors. For the sake of simplicity, we consider only the case of a symmetric and positive 
definite weight matrix W and a full rank matrix A. In fact, our method can also work in the case 
of semi-definite W with rank ( W) 2 m - n. The numerical experiments show that the method can 
also work for some cases in which rank(W) <m - n. Freund showed comparison result for the 
J.Y. Yuan, A.N. IusemlJournal of Computational and Applied Mathematics 71(1996) 287-297 289 
conjugate gradient method and the SOR method for regular least squares problems in [4]. We es- 
tablish here a similar comparison result between the preconditioned conjugate gradient method and 
the 2-cyclic SOR method for problem (1.1), so extending Freund's result [4] to generalized least 
squares problems. Our analysis follows the analysis in [2]. 
In Section 2, a symmetric and positive definite subsystem is derived from the normal equation 
of problem (1.1) by preconditioner techniques. Then the conjugate gradient method is applied to 
solve the subsystem. The error bound of the method for problem (1.1) is given in Section 3. 
Some theoretical comparison results among SOR-type methods and preconditioned conjugate gradient 
method for problem (1.1) are obtained in Section 4. It follows from the comparison results that 
the conjugate gradient method is better than SOR-type methods for problem (1.1) in the same 
Krylov subspace. Finally, some remarks on theoretical results and numerical experiments are given in 
Section 5. 
We define the elliptic norm I1" liD as Ilxllo -- ~ where D is symmetric and positive definite. 
2. Preconditioned conjugate gradient algorithm 
Suppose that A is split in the following way: 
X = (A12), 
where A1 E •nxn is nonsingular and A2 E R (m-")×n. 
(WlI ~r12~ (bl~, 
W--- W r W22) and b= \b2)  
(2.1) 
Then W and b have corresponding splittings: 
(2.2) 
where Wn C ~nxn and W22 E []~(m--n)×(m--n) are symmetric and positive definite, because W is sym- 
metric and positive definite, and W12 C ~,x(m--n). The normal equation of (1.1) is 
W22 W~] r2 = , (2.3) 
A T A T ,] rl 
(rl) where r = = W- l (b -Ax)  is the weighted residual vector corresponding to the splitting form 
r2 
(2.1) of A. We will use preconditioner matrices D1 and D2 defined respectively as 
D1 = I 0 and D2 = I , 
0 AT _pT 
where P ---- AzA~ 1. 
AoI W12 
D71 W22 
It follows from (x) 
W~ D2D~ 1 r2 
A T rl 
= D~ -1 (2.4) 
290 J.Y. Yuan, A.N. IusemlJournal of Computational and Applied Mathematics 71 (1996) 287-297 
that 
i All(/'V12 - I'VllpT) 
W22 - PWIE  - -  ( W T - PWll )pT 
0 
Therefore it holds that 
AlX = bl - (W12 - WllpT)r2, 
(P, - I )w(P ; ) r2=b2-Pbb 
and 
x A~lbl 
W~ -- PWll r2 = b2 Pbl 
r~ + pTr2 ; " 
(2.5) 
(2.6) 
(2.7) 
rl + pTr2 = O. 
It is evident hat the matrix in the left-hand side of (2.7) is symmetric and positive definite, since 
W is symmetric and positive definite. Hence the conjugate gradient method can be applied to the 
reduced system (2.7), and the method converges to the solution of (2.7) in at most m - n steps, in 
the absence of rounding errors. Next, we present he conjugate gradient algorithm for generalized 
least squares problems. 
Algorithm 2.1. 
(1) Factorize A1, and set initial values r~ °) = 0,v (°) = bE - -Pbb p(0) = v(0); 
(2) Iterate for k = 0, 1,... ,  until v (k+l) = 0 (or Ilvek+ )ll tolerance), 
2k-  (p(k),q)' 
D (k+l )  ~--- V (k )  - -  2kq, 
rE k+~ = rE k + 2kp (k), 
IIv¢ +*)ll @ 
p(k+l) = v(k+l) + ~k+lP(k); 
(3) Solve the extra subsystem 
AlX = bl + (Wll PT - W12)r~ l),
where r~ ° is the solution obtained in step 2. 
J.Y. Yuan, A.N. IusemlJournal of Computational nd Applied Mathematics 71 (1996) 287-297 291 
3. Error bound 
In this section, we estimate the error bound for Algorithm 2.1. We need three preparatory results. 
Lemma 3.1. Let S be an arbitrary m x n matrix, and Q an m x m symmetric matrix. Then 
x(Sr QS) <<. x(Q)x(ST S), (3.1) 
where x(B) is the spectral condition number of matrix B. 
Proofi Since STQS and Q are symmetric, there exist orthogonal matrices U, V E R 'n×'~ such that 
uTsTQsu = X = diag (al, a2,. . . ,  ak, 0 , . . . ,  0), (3.2) 
and 
VTQV = A = diag (21,.. .  , 21, 0 , . . . ,  0), 
where ~1 ~>a21> "'" ~>ak > 0 and 21 >~ 22 i> . "  ~> At > 0. Hence 
al = max xTXx = max (VSUx)TA(VSUx) ~< 21 max (Ux)rSrS(Ux), 
Ilxll=l Ilxll=l Ilxll=a 
and 
ak = min xrXx = min (VSUx)rA(VSUx) t> At min (ux)rsTs(Ux).  
Ilxll=l Ilxll=l Ilxll=l 
It follows from (3.4) and (3.5) that 
x(STQS) = cr_La <<x(Q)x(sTs). [] 
ffk 
Lemma 3.2. 
x (( P - I  ) w (P~ ) ) ~< (l + llPll2)x(w), 
where W ~ R mxr" is symmetric, and P c R ("-")×". 
Proof. Apply Lemma 3.1 with S = - I  ' and note that x(SrS)~<l + IlPll . 
Lemma 3.3. Let the function f (x)  be defined by 
E(I+ 1 x- -1  f (x)  = x/x)2 
[ 1 x- -1  1+ (1 + V,~)2 
(3.3) 
(3.4) 
(3.5) 
(3.6) 
[] 
(3.7) 
292 J.Y. Yuan, A.N. IusemlJournal of Computational nd Applied Mathematics 71 (1996) 287-297 
for x >>. 1. Then f (x )  is an increasing function of  x for all x >1 1 and any k >~ O. 
Proof. Taking the derivative in (3.7) we obtain 
l - y2  k (v /x - l )  k-1 
(1 + fY)2 x /~-  ~ ~/~- i  >/0, f ' (x )  - -  
where 
[ x -1  l' [~- l l  k 
y(x) = [(1 ~_-~/~)2j = [1 + v,~j ' 
for all x >/1 and all k/> 0, i.e., f (x )  is an increasing function of x for all x i> 1 and any k >f 0. [] 
Now we can establish the desired error bound. 
Theorem 3.4. The standard error bound based on the Chebyshev polynomials for the CG method 
applied to the generalized linear least squares problem is given by 
HA(x* - x(k))ll ~<2 
I IA(x* - - ,  
1+ 
( l+~2) f l _  1 ]k 
(1 + V/(1 + ~2)fl)2 
(1 + ~2)fl_ 1 ]2k, 
(1 + X/(1 + ~2)fl)2 J
(3 .8 )  
where 0~ = Ilell= = IIAEA11 [[2, f l :  x (w)  = #m~x(W)/bt~n(W) is the spectral condition number of  the 
symmetric and positive definite matrix W, x (°) is a vector corresponding to an arbitrary initial 
vector r~ °) and x* is the solution of  the problem (1.1). 
Proof. It is well known [1] that the conjugate gradient method applied to system (2.7) computes 
from any starting vector r~ °) E ~m-,, with residual r0 = b2 - Pbl - Er~ °), where 
E=(P  - I )W - I  ' 
and further iterates r~ k) E R m-" with the following minimization property in terms of the norm Ilr211e: 
[Ir6 k) - -  r lle = min IIr2 - r~ne, (3.9) 
rEEr~°) +Kk(ro',E) 
r2 - r~ °) ~ Kk(ro;E), 
where r~ is the exact solution of the reduced system (2.7) and 
Kk ( ro ; E ) = Span { ro, Ero, . . . , Ek -  l ro } 
is the kth Krylov subspace of R m-~ generated by r0 and E. 
(3.10) 
(3.11) 
J.Y. Yuan, A.N. IusemlJournal of Computational nd Applied Mathematics 71 (1996) 287-297 293 
It is well known that the error bound of the regular conjugate gradient method applied to a 
symmetric and positive definite system Sx = d, is given 
IIx k) - x*l12 2 L~J  
IIx °)- x*ll: [ X/~ - 1] 2k' 
1+ Lv +l j 
where x is the spectral condition number of S [1]. It follows from Lemma 3.2, and Lemma 3.3 and 
the bound above it 
[ ( l+~2) f l -1  ] k 
IIr  - r~llE (1 + v/(1 + ~2)fl)2 
~< 2 (3.12) 
(1 + ~2)fl_ 1 ]2k. 
[[r~°)-r~lle 1 + (1 +~(1-~) - f l )2 j  
Since 
IIr2 - .2  (P ; )  * r2 lie = (r: - r~)T (p,  - I  ) W (r2 - r 2 ) (3.13) 
and 
(P ; )  * (PT( r2 - r~) )  (3.14, 
(r2 - r2 ) --- - ( r2 - r~ ) ' 
(rl) 
it follows from rl = --pTr2 and r = that 
r2 
- r2 lie = ( r  - r*)TW(r -- r*) = (x --x*)TATW-1A(x --X*) (3.15) lit: * 
using the fact that r - r* = W-1A(x - x*). The result follows replacing (3.15) in (3.12). 
In Algorithm 2.1, we do not compute effectively the matrix P = A2A? ~. Instead of computing it, 
we solve the subsystem A~y = c and perform the matrix-vector product z = .42y or ATy = ATc by 
some direct method. Observe that for many problems n is much smaller than m, so that the size of 
A~ is also much smaller than the size W. In these cases, it is much easier to factorize Al than W 
(or B), as is needed when using BNP algorithm [2], Freund's method [4], or James' method [5]. 
How to get a preconditioner A~ from a general sparse matrix A will be discussed in [3]. 
4. Comparison results 
In the previous section, we presented the conjugate gradient method for problem (1.1). In this 
section, we will compare this algorithm with the 2-cyclic SOR algorithm for problem (1.1) presented 
in Section 1. We start with some known results. 
Lemma 4.1. Consider the linear system 
( I  - B )x  = c.  
294 J.Y. Yuan, A.N. Iusem/Journal of Computational nd Applied Mathematics 71 (1996) 287-297 
Let B be a weakly 2-cyclic matrix. Suppose that the initial vector x (°~ is partitioned in the same 
way as B so that 
Then the SOR method for (I - B)x -= c generates iterates x (k) with 
(a) x} 0 = x} °~, x~ 1~ = x~ °l - cou, 
and for k = 2, 3,... 
(b) xl k) E x} O) -~ Kk_l(BlU;Bl,2); 
E °) + Kk(u; 
where 
(0  B1) (4.1) 
B= B2 0 ' 
B1,2 = B1B2 and B2,1 = B2B~, (4.2) 
and Kk(c;B) = Span{c, Bc, B2c .. . .  ,Bk-lc} is the kth Krylov subspace generated by the vector c 
and the matrix B. 
Proof. See [4]. [] 
For the system corresponding to Algorithm 1.1, the matrix B of the form is 
(0  B1) (4.3) 
B= Ba 0 ' 
where 
( -A l l (W12-Wl lPT)  ) (4.4) 
B1 = __pT 
and 
BE = -- ( W~1A2 W~IW T ) (4.5) 
with the notation of Lemma 4.1. Hence we can extend Freund's result to our problem. 
Theorem 4.2. Suppose that 2-cyclic SOR algorithm is started with initial vectors x (°), r} °) and r~ °). 
Then the method generates iterates x (k), rl k) and r~ k) with 
x (1) = x(°), rl 1) = rl °), r2 (1) = r~ °) - coy(°), 
where 
v (°) = W~lb2 - r~ °) - W~(Azx (°) + W~rl°)), (4.6) 
J.Y. Yuan, A.N. Iusem/Journal of Computational and Applied Mathematics 71 (1996) 287-297 295 
and for k = 2, 3,... 
X (k) E x (°) A-AI-I(wI2- WIlpT)Kk_I(v(O);E), (4.7) 
rl k) E r~ O) + pTKk_l(V(°); E), (4.8) 
r~ k) E r~ °) + K~(v(°);E), (4.9) 
where P = AEA11, and 
E = W~l[P(W12 - Wll PT) + wTpT]. (4.10) 
Proof. By Lemma 4.1 with 81 and BE as in (4.4) and (4.5) respectively, we have 
r~ k) E r~ °) + Kk(v(°);Bzl), (4.11) 
and 
x k) f x (°) 
r~k) C ~ r~O) ) -}- Kk_l(BlV(°);B12 ). (4.12) 
Now (4.9) follows from 
B21 = W~a[P(W12 - WlIP T) + W~PT]. (4.13) 
Since 
fAll(W12 - WllpT)w~21A2 All(W12 -- WI1pT)w~Iw T)  
B12---- ~k pTW221A2 pTW~21W~2 (4.14) 
and for any index v >/0 
(B12fB1 = Bl(B21 f ,  
we get 
Kk-  1 (81V(°); B12 ) = Span{B1 v(°), B1281 v(°), • • •, (812)k-2B1V (0) } 
---- B1Span{v  (°), 821V(0), • • •, (B21)k-2/)(0) 
~- B1Kk-I(V(°); B21 ). 
The proof is complete. [] 
Corollary 4.3. I f  the 2-cyclic SOR algorithm 1.1 and the conjugate gradient algorithm 2.1 are 
both started with the same vector x (°) E •", then 
lib (k) ~..(k+l) -Axc~llw-,<<.[[b-..~SOR rv-,, k = 0,1,... (4.15) 
where ~(k) and ~(k) - c~ SOR are generated by the conjugate gradient algorithm 2.1 and the 2-cyclic SOR 
algorithm 1.1 respectively. 
296 J.Y. Yuan, A.N. IusemlJournal of Computational nd Applied Mathematics 71 (1996) 287-297 
Proof. In the conjugate gradient algorithm 2.1, r~ k~ has the following property 
r~ k) = arg minr2er~0)+/q(~(0); -F)1[r2 - -  r~ lIE 
= arg mint2 e~0,+K,(~,0); F>II"2 - r~llE, 
where 
F = W~21[p(w12 -ml lP  T) + wTpT], 
and E = I - F. Since 
IIr~ k~-  r~llE = lib - Ax(k)llw-,, 
the result follows immediately from Theorem 4.2. [] 
(4.16) 
5. Remarks and conclusions 
From the analysis of the convergence behavior of the conjugate gradient method, we know that 
the error bound given in (3.8) still works in the presence of roundoff errors in which case we may 
have r (n) # 0. 
If W =1, i.e., for regular linear least squares problems, and x(W)= 1, it follows from (3.6) that 
~(E).< 1 + IIPIIL Therefore (3.8) reduces to 
E  ]2k IIA(x* - x~k~)ll= ~ 2 1 + 
IIA(x* -x~°~)ll= 1 + 1 + 
(5.1) 
which coincides with the result given by Freund in [4]. 
Algorithm 2.1 obtains an (m-  n)-dimensional vector 1"2 by the conjugate gradient method, and 
then solves the extra subsystem (2.7). The algorithm requires only AI 1 and W but not W -1. So this 
algorithm is more efficient in actual applications than Freund's algorithm, which requires W -1. In 
the case of W12 --- 0, we can apply the conjugate gradient method to the reduced systems given in 
[9] and get another preconditioned conjugate gradient algorithm which will reduce to other known 
preconditioned conjugate gradient algorithms. Theorem 4.2 and Corollary 4.3 prove that Freund's 
result [4] is true also for the generalized least squares problem. The result in Corollary 4.3 reduces 
to Freund's result and James' result [5] respectively for the regular least squares problem and the 
equality constrained least squares problem. 
Our results how that the preconditioned conjugate gradient algorithm is preferable to SOR method 
for solving generalized least squares problems. 
Numerical results presented in [9] show that though the preconditioned conjugate gradient method 
is numerically less stable than Paige's method, it is much faster than this procedure, especially for 
sparse problems. Also the method works for semi-definite W, even when rank(W) < m-  n, in 
which case the BNP algorithm does not work. 
J.Y. Yuan, A.N. IusemlJournal of Computational nd Applied Mathematics 71 (1996) 287-297 297 
Acknowledgements 
We thank Dr. ]k. Bj6rck, who conjectured the result of Corollary 4.3 and encouraged us to study 
the issue. 
References 
[1] O. Axelsson and G. Lindskog, On the rate of convergence of the preconditioned conjugate gradient method, Numer. 
Math. 48 (1986) 499-523. 
[2] J.L. Barlow, N.K. Nichols and R.J. Plemmons, Iterative methods for equality constrained least squares problems, 
S lAM J. Sci. Statist. Comput. 9 (1988) 892-906. 
[3] /~. Bjfrck and J.Y. Yuan, Preconditioner for least squares problems by LU factorization, to be published. 
[4] R. Freund, A note on two block SOR methods for sparse least-squares problems, Linear Algebra and its Appl. 88/89 
(1987) 211-221. 
[5] D.J. James, Conjugate gradient methods for constrained least squares problems, Ph.D Dissertation, Dept. of Math., 
North Carolina State University, Raleigh, 1990. 
[6] D.J. James, Implicit nullspace iterative methods for constrained least squares problems, S lAM J. Matrix Anal. Appl. 
13 (1992) 962-978. 
[7] C.C. Paige, Fast numerically stable computations for generalized least-squares problems, S lAM J. Numer.Anal. 16 
(1979) 165-171. 
[8] C.C. Paige, Computer solution and perturbation analysis of generalized linear least-squares problems, Math. Comp. 
33 (1979) 171-184. 
[9] J.Y. Yuan, Iterative methods for generalized least squares problems, Ph.D. Thesis, IMPA, Brazil, 1993. 


file:///data/remote/core/dit/data/elsevier/pdf/b2e/aHR0cDovL2FwaS5lbHNldmllci5jb20vY29udGVudC9hcnRpY2xlL3BpaS8wMzc3MDQyNzk1MDAyMzkx.pdf

Preconditioned conjugate gradient method for generalized least squares problems

Abstract

Similar works

Full text

Available Versions

Elsevier - Publisher Connector

Elsevier - Publisher Connector