An algorithm for the least-squares solution of rank-deficient linear - - PowerPoint PPT Presentation

An algorithm for the least-squares solution of rank-deficient linear systems

• G. Rodriguez

Dipartimento di Matematica e Informatica, Universit` a di Cagliari viale Merello 92, 09123 Cagliari, Italy in collab. with Antonio Aric`

• (Dip. di Matematica e Informatica, Cagliari)

Applied Inverse Problems 2009 Computational Aspects and Emerging Applications Vienna, July 20–24 2009

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Displacement structure

The idea was originally conceived for Toeplitz matrix ∆ξ,η(T) := ZξT − TZη = GTH∗

T

T ∈ Cm×n, G ∈ Cm×δ, H ∈ Cn×δ, and Zφ =        · · · φ 1 · · · 1 · · · . . . . . . . . . . . . · · · 1        . [Kailath et al. 1979]

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Displacement structure

The idea was originally conceived for Toeplitz matrix ∆ξ,η(T) := ZξT − TZη = GTH∗

T

Toeplitz-like matrices form an algebra any Schur complement inherits the displacement structure it is possible to perform an LU factorization operating only on the generators with complexity O(n2) (generalized Schur algorithm)

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Cauchy-like matrices & pivoting

It was later noted that Cauchy-like structure is pivoting-invariant various classical structures (Toeplitz, Hankel, Vandermonde) can be efficiently converted into Cauchy-like [Heinig 1995], [Gohberg, Kailath, Olshevsky 1995] Cij = φ∗

i · ψj

ti − sj , DtC − CDs = GCH∗

C

Dt = diag(t1, . . . , tm), G ∗

C =

• φ1

· · · φm

• This leads to a fast and more stable algorithm, known as GKO,

with complexity O(n2) and storage O(n2).

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Schur complements of augmented matrices

Various matrix problems can be reduced to the computation of Schur complements of suitable augmented matrices [Kailath, Chun 1994]. Tikhonov regularization J(λ, x) = Ax − b2 + λ2Hx2, λ > 0     Im A Ip λH A∗ λH∗ A∗b In    

Sm+n+p

− − − − → (A∗A + λ2H∗H)−1A∗b = xλ

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Schur complements of augmented matrices

Generalized Cross Validation (GCV) V (λ) =

1 m b−Axλ2

[ 1

m trace(I−A(λ))] 2 ,

A(λ) = A(A∗A + λ2H∗H)−1A∗     Im A Ip λH A∗ λH∗ A∗ A    

Sm+n+p

− − − − → A(A∗A + λ2H∗H)−1A∗

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Schur complements of augmented matrices

Multi-parameter regularization J(λ, x) = 2

• Ax − b2 + λ2Hx2 + µ2Kx2

[Brezinski, Redivo-Zaglia, R, Seatzu Numer. Math.]       Im A Ip λH Iq µK A∗ λH∗ µK ∗ A∗b In      

S

− → (A∗A + λ2H∗H + µ2K ∗K)−1A∗b

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Square linear systems

Let Cx = b, with C ∈ Cn×n and b ∈ Cn. Then, AC,b = C b −I O

• ⇒

Sn(AC,b) = C −1b. In [Aric`

• , R 2009] this method has been implemented in

Matlab/C-MEX with O(n2) complexity and O(n) storage. This suite of programs allows to convert various structures into Cauchy-like; works naturally with complex matrices; heavily based on the BLAS library; various pivoting strategies (partial, total, S&B’s, Gu’s); Toeplitz-like: no assumptions; Cauchy-like: reconstructability, can deal with multiple nodes.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Execution time

128 256 512 1024 2048 4096 8192 16384 10

−5

10

−4

10

−3

10

−2

10

−1

10 10

1

10

2

10

3

partial Gu toms729 LU α n2

Random Toeplitz matrix - n = 2k, k = 7, . . . , 14

toms729 only works with real matrices drsolve always uses complex arrays & arithmetics

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Errors w.r.to a prescribed solution

128 256 512 1024 2048 4096 8192 16384 10

−15

10

−14

10

−13

10

−12

10

−11

10

−10

10

−9

10

−8

no piv. partial total Gu toms729 LU

Random Toeplitz matrix - n = 2k, k = 7, . . . , 14

toms729 only works with real matrices drsolve always uses complex arrays & arithmetics

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Least squares by Gaussian elimination

Ax = b, A ∈ Rm×n, rank(A) = n ≤ m Peters-Wilkinson method [1970] PAQ = L1 L2

• U,

L1, U ∈ Rn×n LTLy = LT(Pb), U(QTx) = y [Bj¨

• rck et al. 1988] contains an application to sparse matrices.
• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Least squares by Gaussian elimination

Ax = b, A ∈ Rm×n, rank(A) = n ≤ m Peters-Wilkinson method [1970] Augmented matrix method [Siegel 1965, Bartel et al. 1970, Bj¨

• rck 1992]

I A A∗ r x

• =

b

• ,

r = b − Ax Successfully applied in [Arioli et al. 1989] and [Duff et al. 1990] to sparse matrices, with Bunch-Kaufman pivoting.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Least squares via Schur complement

This method was originally proposed for Toeplitz matrices in [Kailath et al. 1994]. It was further studied (and implemented) in [R 2006] for Toeplitz- and Cauchy-like matrices. Ax = b − → MA =   Im A A∗ A∗b In   Sm+n(MA) = (A∗A)−1A∗b = xLS Since MA has displacement structure, it is possible to apply the generalized Schur algorithm. Moreover, the pseudoinverse matrix can be constructed b = Im ⇒ Sm+n(MA) = A†.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Computation of Schur complement of MA

Computation scheme the structure is converted from Toeplitz to generalized Cauchy, via fast transforms the generalized Schur algorithm is applied

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Computation of Schur complement of MA

Computation scheme the structure is converted from Toeplitz to generalized Cauchy, via fast transforms the generalized Schur algorithm is applied Advantages MA, has displacement structure, but it’s not Toeplitz-like, while MC is Cauchy-like Cauchy-like structure allows to apply partial pivoting fast algorithm O(27r(m + n)2) flops (r = rank∆(MA)) small storage required O(3r(m + n)) easily applicable to other kinds of displacement structure

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting restriction

While computing the Schur complement Sr(M), the action of partial (or total) pivoting must be restricted to the first r rows (and columns) of M. In fact, if M = M11 M12 M21 M22

• ,

P = Pr Im−r

• e

Q = Qr In−r

• then

PMQ = PrM11Qr PrM12 M21Qr M22

• and

Sr(PMQ) = Sr(M).

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Partially reconstructable blocks

Cij = φ∗

i · ψj

ti − sj DtC − CDs = GCH∗

C

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Partially reconstructable blocks

Cij = φ∗

i · ψj

ti − sj DtC − CDs = GCH∗

C

⇒ DLMC − MCDR = GMC H∗

MC

MC =   Im C C ∗ C ∗ In   , DL = Dt ⊕ Ds ⊕ Ds, DR = Dt ⊕ Ds ⊕ Dt,

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Partially reconstructable blocks

Cij = φ∗

i · ψj

ti − sj DtC − CDs = GCH∗

C

⇒ DLMC − MCDR = GMC H∗

MC

MC =   Im C C ∗ C ∗ In   , DL = Dt ⊕ Ds ⊕ Ds, DR = Dt ⊕ Ds ⊕ Dt, Some blocks are partially reconstructable: block (3, 2) − → DsIn − InDs = 0 It is necessary to store separately non recostructable entries and keep pivoting into account.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η

Zξ,mA − AZη,n = GAH∗

A,

ξ, η ∈ C \ {0}

Zφ,k =        · · · φ 1 · · · 1 · · · . . . . . . . . . . . . · · · 1        ∈ Ck×k.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η

Zξ,mA − AZη,n = GAH∗

A,

ξ, η ∈ C \ {0} DtC − CDs = GCH∗

C

Cij = φ∗

i · ψj

ti − sj , tm

i

= ξ, sn

j = η

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η

Zξ,mA − AZη,n = GAH∗

A,

ξ, η ∈ C \ {0} Cij = φ∗

i · ψj

ti − sj , tm

i

= ξ, sn

j = η

Setting ξ = 1 e η = eiπϕ, let ϕ∗ = arg max

ϕ min i,j |ti − sj| .

Then, ϕ∗ = gcd(m, n) m , min

i,j |ti − sj| = 2 sin πϕ∗

2n .

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Optimal values of ξ and η

0.3333 0.6667 10

−9

10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

φ m=3000, n=1000, φ*=0.3333 TLLS TLLS no piv.

KMS matrix (aij = ρ|i−j|) with (m, n) = (3000, 1000), ρ = .99, κ = 3.7 · 104

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Numerical results

1000 2000 3000 4000 5000 10

−16

10

−14

10

−12

ρ=0.5 2 4 6 8 10 10

−16

10

−14

10

−12

1000 2000 3000 4000 5000 10

−10

10

−8

10

−6

10

−4

n ρ=0.99 2 4 6 8 10 10

−10

10

−8

10

−6

10

−4

α TLLS TLLS no piv. QR Cholesky

KMS matrix, m = 2n (left), α = m

n e n = 2000 (right)

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Influence of conditioning

Conditioning of xLS with respect of perturbations in A κLS = κ(A) + κ(A)2 tan θ η with κ(A) = A · A†, η = A·x

Ax

∈ [1, κ(A)], and θ = arccos Ax b = ∡(b, Ax) ∈

• 0, π

2

• QR solution is influenced by κLS

Normal eqs. + Cholesky depend on κ(A)2, and not on η, θ TLLS ?

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Influence of conditioning

20 40 60 80 10

−10

10

−8

10

−6

10

−4

10

−2

θ (degrees) TLLS QR Cholesky 20 40 60 80 10

−10

10

−9

10

−8

10

−7

10

−6

θ (degrees) TLLS QR Cholesky

KMS matrix, (m, n) = (1500, 500), ρ = .99, .999 κ(A) = 3.7 · 104, 1.1 · 106, η ≃ 1

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Computational complexity - unstructured methods

500 1000 1500 2000 10

3

10

6

10

9

10

12

n flops complexity 500 1000 1500 2000 10

−4

10

−2

10 10

2

n seconds computation time TLLS Householder QR Cholesky TLLS Householder QR Cholesky

Comparision between Matlab qr and chol, and the MEX version of TLLS

AMD Athlon 64 3200+, GNU/Linux x86 64, Matlab 64bit KMS matrix (ρ = .99), m = 2n

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Computational complexity - superfast method

5 7 9 11 13 10

−15

10

−10

10

−5

log2 n relative error TLLS TLLS no piv. superfast 5 7 9 11 13 10

−2

10 10

2

10

4

log2 n time in seconds TLLS TLLS no piv. superfast

Comparision between the Matlab versions of TLLS and a superfast algorithm [Van Barel, Heinig, Kravanja 2003]

rand-Toeplitz matrix, m = 2n

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Complexity: crossover w.r. to unstructured methods

The dependance on m is linear for unstructured methods and quadratic for TLLS. Cholesky − → n2(m + 1 3n) Householder QR − → 2n2(m − 1 3n) TLLS − → 143m2 + 284mn + 133n2 This means that for each n, when m grows, there is a point after which unstructured metods become more efficient.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Complexity: crossover w.r. to unstructured methods

1 20 40 60 80 100 1 2 3 4 5 6 2 4 6 8 10

α log10 n

log(complexity/n2) for Householder QR (grid) and TLLS (surf)

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Rank-deficient linear systems

Let us drop the assumption that the matrix A is

• verdetermined and full-rank

Ax = b, A ∈ Rm×n,

• m ⋚ n

rank(A) ≤ min(m, n) If we apply the TLLS algorithm MA =   Im A A∗ A∗b In   − → Sm+n(MA) =? In the full-rank case we have xLS = Sm+n(MA) = (A∗A)−1A∗b

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Phase 1

We perform Gauss elimination to the first m columns. Pivoting is not essential to perform the computation.   Im A A∗ A∗b In   − →   Im A −A∗A A∗b In   This phase is equivalent to the construction of normal equations, unless pivoting is applied (more on this later). A∗A is singular in general

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Phase 2

Let us first assume that A∗A is well ordered: −A∗A = G11 G12 G21 G22

• with rank(A) = k and G11 k × k non singular.

We perform k steps of Gauss algorithm to the submatrix   Im A −A∗A A∗b In   − → −A∗A A∗b In

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Phase 2

Let us first assume that A∗A is well ordered: −A∗A = G11 G12 G21 G22

• with rank(A) = k and G11 k × k non singular.

We perform k steps of Gauss algorithm to the submatrix −A∗A A∗b In

• =

    G11 G12 c1 G21 G22 c1 Ik In−k    

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Phase 2

Let us first assume that A∗A is well ordered: −A∗A = G11 G12 G21 G22

• with rank(A) = k and G11 k × k non singular.

We perform k steps of Gauss algorithm     G11 G12 c1 G21 G22 c1 Ik In−k     − →     Uk Z

• c1

On−k

• c2

Y1 x1 Y2 x2    

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

What the algorithm produces

    Uk Z

• c1

On−k

• c2

Y1 x1 Y2 x2     basic solution of Ax = b xB = x1 x2

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

What the algorithm produces

    Uk Z

• c1

On−k

• c2

Y1 x1 Y2 x2     basic solution of Ax = b xB = x1 x2

• basis for the kernel of A

Y = Y1 Y2

• =
• y1

· · · yn−k

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

What the algorithm produces

    Uk Z

• c1

On−k

• c2

Y1 x1 Y2 x2     basic solution of Ax = b xB = x1 x2

• =
• −G −1

11 c1

• basis for the kernel of A

Y = Y1 Y2

• =
• y1

· · · yn−k

• =
• −G −1

11 G12

In−k

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting in phase 2

In phase 2 pivoting (on both rows/columns) is essential Assume to apply total pivoting −A∗A A∗b In

• −

→ −PA∗AQ PA∗b Q

• and let

G = −PA∗AQ, Q =

• Q1

Q2

• ,

c = PA∗b. Then,   G11 G12 c1 G21 G22 c1 Q1 Q2   − →   Uk L−1G12 L−1c1 On−k c2 − G21G −1

11 c1

Q2 − Q1G −1

11 G12

−Q1G −1

11 c1

 

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting in phase 2

  Uk L−1G12 L−1c1 On−k c2 − G21G −1

11 c1

Q2 − Q1G −1

11 G12

−Q1G −1

11 c1

  basic solution of Ax = b xB = −Q1G −1

11 c1 = Q

• −G −1

11 c1

• ,
• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting in phase 2

  Uk L−1G12 L−1c1 On−k c2 − G21G −1

11 c1

Q2 − Q1G −1

11 G12

−Q1G −1

11 c1

  basic solution of Ax = b xB = −Q1G −1

11 c1 = Q

• −G −1

11 c1

• ,

basis for the kernel of A Y = Q2 − Q1G −1

11 G12 = Q

• −G −1

11 G12

In−k

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting strategies

total pivoting (expensive) partial pivoting + “empirical rules” Bunch-Kaufman pivoting Gu’s pivoting

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Rank revealing?

Numerical rank rǫ = min

E2≤ǫ rank(A + E)

⇔ σrǫ > ǫ ≥ σrǫ+1 The algorithms stops when |Mrr| < τ What is the relation between ǫ and τ?

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting in phase 1

In phase 1 pivoting is not indispensable M =   Im A A∗ A∗b In  

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting in phase 1

In phase 1 pivoting is not indispensable, but it’s essential for the stability of the algorithm. Mγ =   γIm A A∗

1 γ A∗b

In   The scaling parameter γ influences the action of pivoting.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Pivoting in phase 1

In phase 1 pivoting is not indispensable, but it’s essential for the stability of the algorithm. Mγ =   γIm A A∗

1 γ A∗b

In   The scaling parameter γ influences the action of pivoting. [Bj¨

• rck 1992] proved that the optimal value of γ for the method of

the augmented matrix with scaling is γ = 2−1/2σn(A) since it minimizes κ2 γIm A

A∗ 0

• (min ≃

√ 2κ2(A)), and hence the estimate of the forward error on the solution.

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

LDLT factorization

It preserves the symmetry and it is connected to Cholesky and QR.   Im A A∗ A∗b In  

1

→   Im −A∗A A∗b In  

2

→   Im D xLS   Symmetric row/column pivoting if at step k the maximum is at position (ℓ, ℓ) pivot = Mℓℓ if the maximum is at position (r, s) pivot = Mss Msr Mrs Mrr

• (tile pivot)
• G. Rodriguez

Least-squares solution of rank-deficient linear systems

Yet to be done. . .

investigate the role of the scaling parameter γ investigate the effect of various pivoting strategies determine an effective stopping criterion normal solution, up/down-dating, subset selection. . . sparse matrices? implementation in the structured case extensive numerical tests, in comparison with other methods (structured and unstructured)

• G. Rodriguez

Least-squares solution of rank-deficient linear systems

The end

Thanks!

• G. Rodriguez

Least-squares solution of rank-deficient linear systems