Preconditioners for the Least Squares Problems
Based on the Incomplete Givens Orthogonalization
Jun-Feng Yin
National Institute of Informatics

Nov. 16, 2006
yinjf@nii.ac.jp

Joint work with
♦ Zhong-Zhi Bai
Institute of Computational Mathematics
Chinese Academy of China, China♦ Ken Hayami
National Institute of Informatics, Japan
Research subsidized byFirst NII International Internship Program(2006, 1 – 2006, 3)

Outline

IntroductionIncomplete Givens Orthogonalization (IGO) MethodPreconditioned Iterative Method• CGLS• LSQR• BA-GMRES

Numerical ResultsConclusion

Problems
Consider the solution of linear least squares problems
min kb − Axk,

x∈RA ∈ Rm×n.

where A is large and sparse, which arise from many scientific applications,• Statistical model• Linear and nonlinear programming• Signal and image processing

Normal Equation
The iterative method solves the normal equation
AT Ax = AT b.
The coefficient matrix of the equation becomes:• Symmetric positive definite
suitable for conjugate gradient method.
• More larger condition number
require the efficient preconditioners.

Recently Contributions♦ Iterative Method:• CGLS M.R. Hestenes and E. Stiefel,
1952.
• LSQR C.C. Paige and M.A. Saunders,1982.
• BA(AB)-GMRES K. Hayami and T. ITO, 2005.
♦ Preconditioning techniques:• CIMGS X. Wang and et al., 1997.
• IGO Z.Z. Bai, I.S. Duff and A.J. Wathen,• RIF M. Benzi and M. Tuma, 2003.
General View

Å. Björck’s Book,

1996.

2001.

Review of Givens Rotation
A Givens rotation (or plane rotation) G(i, j, θ) ∈ Rn×n is
equal to the identity matrix except thatc s,G([i, j], [i, j]) =−s c
for k 6= i, j,xk ,Algebraically, yk = cxi + sxj ,for k = i, 1 ≤ k ≤ n,−sxi + cxj , for k = j,and so, yj = 0 if
xjxi
s=q,
c=q.
x2i + x2j
x2i + x2j

The Process of the IGO Method  ∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗
→ ∗ ∗ ∗ ↓ ∗ ∗ ∗ ↑ ∗ ∗ ∗ or 
 or 

→ → ∗ ∗ ↓ ↓ ∗ ∗ ↑ ↑ ∗ ∗→ → → ∗↓ ↓ ↓ ∗↑ ↑ ↑ ∗Choice of the order depends on• Storage format• Computational implementation

Row IGO Method
Method : The General row IGO Method
1. For i = 2, . . . , m Do:
2. For j = 1, . . . , min{i − 1, n} and aij 6= 0 Do:
3. If |aij | < τ , Then
4. Set aij = 0
5. Cycle to Next j
6. EndIfq
7. Compute ρ := a2jj + a2ij
8. Compute c := ajj /ρ
9. Compute s := aij /ρ
10. Set ajj = ρ , aij = 0
11. Store c, sRow IGO Method(continue)
12. For k = min{j + 1, n}, . . . , n Do:
13. Compute tempj := cajk + saik
14. Compute temp i := −sajk + caik
15. If |ajk | > τ , Set ajk := tempj
16. Else, Set ajk := 0
17. If |aik | > τ , Set aik := temp i
18. Else, Set aik := 0
19. EndDo
20. EndDo
21. EndDo
22. For i = 1, . . . , n
23. Keep only the p largest entries in the i-th row
24. Set rij := aij
25. EndDo

Theorem of the IGO-method
Theorem (Bai/Duff/Wathen BIT(2001))
Let A ∈ Rn×n be nonsingular, and Q, R ∈ Rn×n be
the incomplete orthogonal and upper triangular matrices,
respectively, produced by the IGO-method. Then(i) R is sparse and nonsingular, and its diagonal entries
are positive except possibly for the last one;(ii) Q = GT is orthogonal, provided PQ = Pn, where1Q
G = Gn−1Gn−2 · · · G1 ≡Gj .
j=n−1

Theorem of the Row IGO Method
Theorem Let A ∈ Rm×n be full column rank, and Q, R ∈Rn×n be the incomplete orthogonal and upper triangular matrices, respectively, produced by the IGO-method.
Then(i) R ∈ Rn×n is sparse and nonsingular, and its diagonal
entries are positive except possibly for the last one;(ii) Q = GT ∈ Rm×m is orthogonal, where G =1Q
Gn−1Gn−2 · · · G1 ≡Gj .
j=n−1

More Consideration for the IGO Method
For the squares matrices,Gn . . . G1 A ≈ R
For the rectangular matrices R,Gn . . . G1 A ≈0

 RA≈Q, M.
0More details on♦ Sparse structure and storage scheme♦ The magnitude of the fill-ins♦ How to solve the equation M z = v

Sparse Structure
By define
PA,L = {(i, j) | aij 6= 0,PA,R = {(i, j) | aij 6= 0,i > j,i ≤ j,

1 ≤ i, j ≤ n},
1 ≤ i, j ≤ n},

we have the dropping rule
PR = PA,Rand

PG = PA,L

By introducing a parameter l, 0 ≤ l ≤ n, we have
PA,L,l = PA,L ∪ {(i, j) | i ≤ j ≤ i + l, 1 ≤ i, j ≤ n},m
PA,R,l = PA,R ∪ {(i, j) | i − d el ≤ j ≤ i, 1 ≤ i, j ≤ n},nGeneral Remarks• The IGO method will never breakdown and produce an orthogonal matrix Q and an upper triangular matrix R,• The upper triangular matrix R is nonsingular in case of theoriginal matrix is full column rank ,• All the factorization processes are numerically stable,• Need more operations than ILU method.

Preconditioned Iterative Method

Preconditioned CGLS Method
M = RTR,R−TATAR−1u = R−TATb,

u = Rx.

• Q should not be stored.
• The preconditioner M is symmetric positive definite.
• The additional computation is only solving the two triangleequations:
RTw = v, Rz = w.

Preconditioned LSQR Method
The preconditioned LSQR method solves the problem
min kb − AR−1xk,
instead ofmin kb − Axk,
which imply
θ1v1 = ATb, ρ1p1 = AR−1v1.
T −Tθi+1vi+1 = A R pi − ρivi
, i = 1, 2, . . .
−1
ρi+1pi+1 = AR vi − θi+1pi

Preconditioned BA-GMRES Method
The method solve
BAx = Bbinstead of the normal equation, where B should satisfied:
R(A) = R(B T),

R(AT) = R(B)

Theorem (Hayami/ITO (2005)) If R(A) = R(B T) andR(AT) = R(B), the BA-GMRES method determines a leastsquares solution of minx∈Rn kb − Axk2 for all b ∈ Rm andfor all x0 ∈ Rn without breakdown.

Choice of the mapping matrix B
When A is full column rank, R is nonsingular, there are atleast three choices:
B = AT(1)

B = [R−1 0]G(2)T

B = (R R)
R(A) = R(B T),−1

AT

R(AT) = R(B)(3)

Theoretical Analysis
R(A) = R(B T) R(AT) = R(B)
√√
T
B=A√
−1
B = [R 0]G?
√√
−1 TTB = (R R) A

Numerical ResultsComparison between B(1) and B(3)CGLS LSQR B(1) PCGLS PLSQR B(3)RANDS2408411 162
3535 360.25
1.47 1.69 1.11
0.26 0.28
RANDS4 1891 1926 239
8687 73
6.81 7.84 2.00
0.53
0.56 0.56
RANDS6 9871 10148 258
325337 130
35.61 42.10 2.23
1.80
2.03 1.09
RANDS8 46520 44634 266
1035
1051 153
168.54 185.11 2.38
5.49
6.11 1.34m = 1000, n = 320, den(A) = 4.9%, CPUP = 0.08 sec

Numerical Comparison of B(2) and B(3)lPCGLS PLSQR B(2) B(3)438431967 256
100 5.04 17.76 17.59 37.94 13.72
387379961 244200 7.12 15.58 15.95 41.91 13.48
333
324837 220
400 12.62 14.39 14.29 42.28 12.23
256248512 188
600 17.75 11.20 11.13 28.14 10.21
188
179
160 142800 22.72 8.38
8.14 9.40 7.45
111
11000 24.91 0.05
0.05 0.07 0.05
Tr
m = 10000, n = 1000, κ(A) = 80, den(A) = 1.5%, A< 10−6AT b

010PCGLSPLSQR
BA−GMRES−2
BA−GMRES−3−5

relative residual

10−1010−15100100200300400500
iteration

600m = 10000, n = 1000, l = 700, κ(A) = 104

700

800

900

010PCGLSPLSQR
BA−GMRES−2BA−GMRES−3−5relative residual

10−1010−1510050100150200250
iteration

300350m = 10000, n = 1000, l = 800, κ(A) = 104400450500

010

PCGLS
PLSQR
BA−GMRES−2BA−GMRES−3−210−410−6relative residual

10−810−1010−1210−1410−1610−18100102030
iteration

40m = 10000, n = 1000, l = 900, κ(A) = 104

5060

Numerical Results when l = 200
PCGLS
PLSQRBAGMRESITCPUITCPU IT CPU
RAND1 387 15.58 379 15.95 244 13.48
RAND2 914 36.81 917 37.74 420 26.94
RAND3 1932 78.16 1942 80.02 598 44.53
RAND4 5708 230.38 5897 243.06 721 58.81
RAND5 15211 613.45 16561 682.27 770 65.17
RAND6 44023 1777.38 50540 2082.52 786 67.06
RAND7 90389 3671.70 ––802 69.61
RAND8 ––––776 66.06
MAXIT=100000,m = 10000, n = 1000, CPUP = 7.12 sec

Conclusion• Introduce a class of preconditioner base on incomplete
Givens orthogonal factorization.
• Analysis their convergence performance of the preconditioned iterative methods for the least square problems.
• Examine their properties numerically, especially for theBA-GMRES method.

Thank You!