[PPT] - Rank Revealing QR factorization F. Guyomarch, D. Mezher and B. PowerPoint Presentation

SLIDE 1

Rank Revealing QR factorization

F. Guyomarc’h, D. Mezher and B. Philippe

1

SLIDE 2

Outline

Introduction
Classical Algorithms

⋆ Full matrices ⋆ Sparse matrices

Rank-Revealing QR
Conclusion

CSDA 2005, Cyprus 2

SLIDE 3

Situation of the problem

See Bjorck SIAM96 : Numerical Methods for Least Squares Problems.

Data : A ∈ Rm×n, b ∈ Rm. Problem P : Find x ∈ S such that x minimum where S = {x ∈ Rn | b − Ax minimum }. (. ≡ .2) The solution is unique : ˆ x = A+b. Property : x ∈ S iff AT(b − Ax) = 0, iff ATAx = ATb. The last equation is consistent since R(ATA) = R(AT).

CSDA 2005, Cyprus 3

SLIDE 4

Rank-Revealing QR factorization

Theorem : A ∈ Rm×n with rank(A) = r (r ≤ min(m, n)). There exist Q, E, R11 and R12 such that

Q ∈ Rm×m is orthogonal,
E ∈ Rn×n is a permutation,
R11 ∈ Rr×r is upper-triangular with positive diagonal entries,
R12 ∈ Rr×(n−r),

and AE = Q

R11

R12

.

RRQR The factorization is not unique. Let RRQR be any of them.

CSDA 2005, Cyprus 4

SLIDE 5

Actually, if we consider a complete orthogonal decomposition A = Q

T
V T

where Q and V are orthogonal and T triangular with positive diagonal entries, we have A+ = V

T −1
QT.

Such a factorization can be obtained from the previous one by per- forming a QR factorization of

RT

11

RT

12

CSDA 2005, Cyprus

5

SLIDE 6

Column pivoting strategy

Householder QR factorization using Householder reflections :

A = v = A1..12,1 + ||A1..12,1||e1 H1 = I12 − 2vvt

vtv

H1A = v = A2..12,2 + ||A2..12,2||e1 H2 = I11 − 2vvt

vtv

CSDA 2005, Cyprus 6

SLIDE 7

Sparse QR factorization

Factorizing a sparse matrix implies fill-in in the factors. The situation is worse with QR than with LU since when updating the tailing matrix :

LU : the elementary transformation x −

→ y = (I − αveT

k )x keeps x

invariant when xk = 0.

QR : the elementary transformation x −

→ y = (I − αvvT)x keeps x invariant when x⊥v. Since ATA = RTR, the QR factorization A is related to the Cholesky factorization of ATA. It is known that a symmetric permutation on a sparse s.p.d. matrix changes the level of fill-in.

CSDA 2005, Cyprus 7

SLIDE 8

Therefore, a permutation of the columns of A changes the fill-in of R. ⇒ there is conflict between pivoting to minimize fill-ins and pivoting associated with numerical properties. We choose to decouple the sparse factorization phase and the rank- revealing phase For a standard QR factorization of a sparse matrix : Multi-frontal QR method [Amestoy-Duff-Puglisi ’94] Routine MA49AD in Library HSL.

CSDA 2005, Cyprus 8

SLIDE 9

Column pivoting strategy

Householder QR factorization using Householder reflections :

A = v = A1..12,1 + ||A1..12,1||e1 H1 = I12 − 2vvt

vtv

H1A = v = A2..12,2 + ||A2..12,2||e1 H2 = I11 − 2vvt

vtv

Householder QR with Column Pivoting (Businger and Golub) : At step k, the column jk ≥ k defining the Householder transformation is chosen so that A(k : m, jk) ≥ A(k : m, j) for j ≥ k.

CSDA 2005, Cyprus 9

SLIDE 10

Some properties on R11

At step k : A(k) ≡ Hk · · · H1AE(k) =

  R(k)

11

R(k)

12

A(k)

22

  .

Any column a of

  R(k)

12 (k, :)

A(k)

22

 

satisfies R(k)

11 (k, k) ≥ a.

Moreover A ≥ R(k)

11 ≥ R(k) 11 (1, 1)

R(k)

11 (k, k) ≥ σmin(R(k) 11 ) ≥ σmin(A),

This implies that : cond(A) ≥ cond(R(k)

11 ) ≥ R(k)

11 (1,1)

R(k)

11 (k,k). CSDA 2005, Cyprus 10

SLIDE 11

Bad new

The quantity R(k)

11 (1,1)

R(k)

11 (k,k) cannot be considered as an estimator of cond(R(k)

11 )

since it can be of different order of magnitude : Example (Kahan’s matrix of order n) : for θ ∈ (0, π

2), c = cos θ and s = arcsin θ:

MATLAB : d=c .ˆ [0:n-1] ; M=diag(d)(eye(n)-striu(ones(n),1)); n=100 ; θ = arcsin(0.2) : σmin(R11) = 3.7e − 9 ≪ R11(100, 100) = 1.3e − 1 Therefore a better estimate of cond(R(k)

11 ) is needed.

CSDA 2005, Cyprus 11

SLIDE 12

Incremental Condition Estimator (ICE) [Bischof 90] which is implemented in LAPACK : it estimates cond(R(k)

11 ) from an estimation of cond(R(k−1) 11

). However, the strategy which consists in stopping the factorization when cond(R(k)

11 ) < 1

ǫ ≤ cond(R(k+1)

11

) may fail in very special situations : Counterexample : M = Kahan’s matrix of order 30 with θ = arccos(0.2) cond(R(16)

11

) < 1

ǫ < cond(R(17) 11

) indicates a numerical rank equal to 16 although the numerical rank of M computed by SVD is 22.

CSDA 2005, Cyprus 12

SLIDE 13

Pivoting strategies

Convert R to reveal the rank by pushing the singularities towards the right end

Apply an Incremental Condition Estimator to evaluate σmax and

σmin of R1→i,1→i,

if σmax/σmin > τ, move column i towards the right end and re-
rthogonalise R using Givens rotations

CSDA 2005, Cyprus 13

SLIDE 14

Pivoting strategies

RRQR uses a set of thresholds to overcome numerical singularities,
A predifined set of thresholds might fail

A =

  

2 10−6 1 1 10−3 −10−3

  

rank=2 if τ = {107}, rank=3 if τ = {104, 107}

RRQR adapts the thresholds to avoid failure

⋆ Upon completion, check if

R22
< σmax

τ

≃ σmin ⋆ On failure, insert new thresholds and restart algorithm

CSDA 2005, Cyprus 14

SLIDE 15

Numerical results

5 10 15 20 25 30 35 40 45 50 10 10

5

10

15

10

20

10

25

10

30

10

35

10

40

10

45

10

50

kahan(50,acos(0.2)), cond=2.4721e+49, rank=20, size=50x50, req=0 condition number smax/smin cond([q,r,p]=qr(A) RRQR, thres=1e3,1e7,1e12,1e25 RRQR, thres=1e3,1e9,1e25

CSDA 2005, Cyprus 15

SLIDE 16

Numerical results

5 10 15 20 25 30 35 40 45 50 10 10

5

10

15

10

20

10

25

10

30

10

35

10

40

10

45

10

50

kahan(50,acos(0.2)), cond=2.4721e+49, rank=20, size=50x50 condition number smax/smin cond([q,r,p]=qr(A) RRQR, thres=1e30 RRQR, thres=1e20 RRQR, thres=1e25, req=2

CSDA 2005, Cyprus 16

SLIDE 17

Second strategy

The k + 1 column is not the worst usually for the condition number. Can we select the worst one for the conditionement of R11 ? Then we reject this worst column to the end.

k+1 j k+1

We need to recompute the singular values and vectors estimates. Reverse ICE

CSDA 2005, Cyprus 17

SLIDE 18

Second strategy

100 200 300 400 500 10 10

5

10

15

10

20

Random matrix, n=500 Exact Estimated tol=1e−12 Estimated tol=1e−7 CSDA 2005, Cyprus 18

Rank Revealing QR factorization

1

Outline

⋆ Full matrices ⋆ Sparse matrices

CSDA 2005, Cyprus 2

Situation of the problem

See Bjorck SIAM96 : Numerical Methods for Least Squares Problems.

Data : A ∈ Rm×n, b ∈ Rm. Problem P : Find x ∈ S such that x minimum where S = {x ∈ Rn | b − Ax minimum }. (. ≡ .2) The solution is unique : ˆ x = A+b. Property : x ∈ S iff AT(b − Ax) = 0, iff ATAx = ATb. The last equation is consistent since R(ATA) = R(AT).

CSDA 2005, Cyprus 3

Rank-Revealing QR factorization

Theorem : A ∈ Rm×n with rank(A) = r (r ≤ min(m, n)). There exist Q, E, R11 and R12 such that

and AE = Q

R12

RRQR The factorization is not unique. Let RRQR be any of them.

CSDA 2005, Cyprus 4

Actually, if we consider a complete orthogonal decomposition A = Q

where Q and V are orthogonal and T triangular with positive diagonal entries, we have A+ = V

Such a factorization can be obtained from the previous one by per- forming a QR factorization of

11

RT

12

5

Column pivoting strategy

Householder QR factorization using Householder reflections :

A = v = A1..12,1 + ||A1..12,1||e1 H1 = I12 − 2vvt

H1A = v = A2..12,2 + ||A2..12,2||e1 H2 = I11 − 2vvt

CSDA 2005, Cyprus 6

Sparse QR factorization

Factorizing a sparse matrix implies fill-in in the factors. The situation is worse with QR than with LU since when updating the tailing matrix :

→ y = (I − αveT

k )x keeps x

invariant when xk = 0.

→ y = (I − αvvT)x keeps x invariant when x⊥v. Since ATA = RTR, the QR factorization A is related to the Cholesky factorization of ATA. It is known that a symmetric permutation on a sparse s.p.d. matrix changes the level of fill-in.

CSDA 2005, Cyprus 7

CSDA 2005, Cyprus 8

Column pivoting strategy

Householder QR factorization using Householder reflections :

A = v = A1..12,1 + ||A1..12,1||e1 H1 = I12 − 2vvt

H1A = v = A2..12,2 + ||A2..12,2||e1 H2 = I11 − 2vvt

Householder QR with Column Pivoting (Businger and Golub) : At step k, the column jk ≥ k defining the Householder transformation is chosen so that A(k : m, jk) ≥ A(k : m, j) for j ≥ k.

CSDA 2005, Cyprus 9

Some properties on R11

At step k : A(k) ≡ Hk · · · H1AE(k) =

  R(k)

11

R(k)

12

A(k)

22

  .

Any column a of

  R(k)

12 (k, :)

A(k)

22

 

satisfies R(k)

11 (k, k) ≥ a.

Moreover A ≥ R(k)

11 ≥ R(k) 11 (1, 1)

R(k)

11 (k, k) ≥ σmin(R(k) 11 ) ≥ σmin(A),

This implies that : cond(A) ≥ cond(R(k)

11 ) ≥ R(k)

11 (1,1)

R(k)

11 (k,k). CSDA 2005, Cyprus 10

Bad new

The quantity R(k)

11 (1,1)

R(k)

11 (k,k) cannot be considered as an estimator of cond(R(k)

11 )

since it can be of different order of magnitude : Example (Kahan’s matrix of order n) : for θ ∈ (0, π

2), c = cos θ and s = arcsin θ:

MATLAB : d=c .ˆ [0:n-1] ; M=diag(d)*(eye(n)-s*triu(ones(n),1)); n=100 ; θ = arcsin(0.2) : σmin(R11) = 3.7e − 9 ≪ R11(100, 100) = 1.3e − 1 Therefore a better estimate of cond(R(k)

11 ) is needed.

CSDA 2005, Cyprus 11

Incremental Condition Estimator (ICE) [Bischof 90] which is implemented in LAPACK : it estimates cond(R(k)

11 ) from an estimation of cond(R(k−1) 11

MATLAB : d=c .ˆ [0:n-1] ; M=diag(d)(eye(n)-striu(ones(n),1)); n=100 ; θ = arcsin(0.2) : σmin(R11) = 3.7e − 9 ≪ R11(100, 100) = 1.3e − 1 Therefore a better estimate of cond(R(k)