10−3 10−2 10−1 100 101
100 101 102 103 104 105
||Axλ−b||2
||xλ||2
L Kurve; Hilbertmatrix der Dimension 100; Stoerung 0.1%
α=1e−2 α=1 α=1e−4
α=1e−6 α=1e−8
α=1e−10 α=1e−12
α=1e−14 α=1e−16
TUHH Heinrich Voss Least Squares Problems Valencia 2010 66 / 82
Toy problem
The following table contains the errors for the linear systemAx =bwhereAis the Hilbert matrix, andbis such thatx =ones(n,1)is the solution. The regularization matrix isL=Iand the regularization parameter is determined by the L-curve strategy. The normal equations were solved by the Cholesky factorization, QR factorization and SVD.
n=10 n=20 n=40 Tikhonov Cholesky 1.41 E-3 2.03 E-3 3.51 E-3 Tikhonov QR 3.50 E-6 5.99 E-6 7.54 E-6 Tikhonov SVD 3.43 E-6 6.33 E-6 9.66 E-6 The following table contains the results for the LS problems (m=n+20).
n=10 n=20 n=40 Tikhonov Cholesky 3.85 E-4 1.19 E-3 2.27 E-3
Large problems
Large ill-posed problems
Until now we assumed that the SVD ofAor TSVD of(A,L)is available which is only reasonable for not too large dimensions.
We now assume that the matrixAis so large as to make the decomposition of its singular value decomposition undesirable or infeasible.
The first method introduced by Bj ¨orck (1988) combines the truncated SVD with the bidiagonalization of Golub-Kahan-Lanczos.
With the starting vectorb, put
β1u1=b, α1v1=ATu1, (47)
and fori =1,2, . . . compute
βi+1ui+1 = Avi−αiui
αi+1vi+1 = ATui+1−βi+1vi (48) whereαi ≥0 andβi ≥0,i =0,1,2, . . . are chosen so thatkuik=kvik=1.
TUHH Heinrich Voss Least Squares Problems Valencia 2010 69 / 82
Bidiagonalization & TSVD
With
Uk = (u1, . . . ,uk), Vk = (v1, . . . ,vk), B¯k =
α1 β1 α2
β3 . ..
. .. αk βk
(49)
The recurrence relation can be rewritten as
Uk+1(β1e1) =b, (50) AVk =Uk+1B¯k, ATUk+1=VkB¯kT +αk+1vk+1eTk+1. (51) We are looking for an approximate solutionxk ∈span{Vk}and write
Large problems
Bidiagonalization & TSVD cnt.
From the first equation of (51) it followsAxk =AVkyk =Uk+1B¯kyk, and since (in exact arithmetic)Uk+1can be shown to be orthogonal, it follows that
kb−Axkk=kUk+1(β1e1−B¯kyk)k=kβ1e1−B¯kykk.
Hence,kb−Axkkis minimized over span(Vk)ifyk solves the least squares problem
kβ1e1−B¯kykk=min!. (52) Withdk :=β1e1−B¯kyk it holds
AT(Axk−b) =ATUk+1dk =VkB¯kTdk+αk+1vk+1eTk+1dk, and ifyk solves (52) we get frombTkdk =0
AT(b−Axk) =αk+1vk+1eTk+1dk. (53)
TUHH Heinrich Voss Least Squares Problems Valencia 2010 71 / 82
Bidiagonalization & TSVD cnt.
Assume thatαi 6=0 andβi 6=0 fori =1, . . . ,k. If in the next stepβi+1=0, thenB¯k has rankk and it follows thatdk =0, i.e.Axk =b. Ifβk+16=0 but αk+1=0, then by (53)xk is a least squares solution ofAx =b.
Thus, the recurrence (51) cannot break down before solution of kAx −bk=min!is obtained.
The bidiagonalization algorithm (47), (48) is closely related to the Lanczos process applied to the symmetric matrixATA. If the starting vectorv1is used, then it follows from (48)
(ATA)Vk =ATUk+1B¯k =Vk( ¯BkTB¯k) +αk+1βkvk+1ekT.
Large problems
Bidiagonalization & TSVD cnt.
To obtain a regulrized solution ofkAx−bk=min!consider the (full) SVD of the matrixB¯k ∈Rk+1×k
B¯k =Pk Ωk
0
QkT =
k
X
i=1
ωipiqiT (54)
wherePk andQk are square orthogonal matrices and
ω1≥ω2≥ · · · ≥ωk >0. (55)
Then the solutionyk to the least squares problem (52) can be written yk =β1Qk(Ω−1k , 0)PkTe1=β1
k
X
i=1
ωi−1κ1iqi (56) withκij = (Pk)ij.
TUHH Heinrich Voss Least Squares Problems Valencia 2010 73 / 82
Bidiagonalization & TSVD cnt.
The corresponding residual vector is
dk =β1Pkek+1eTk+1PkTe1=β1κ1,k+1pk+1, from which we get
kb−Axkk=kdkk=β1|κ1,k+1|.
For a given thresholdδ >0 we define regularized solution by the TSVD method:
xk(δ) =Vkyk(δ), yk(δ) =β1X
ωi>δ
ωi−1κ1iqi. (57)
Notice, that the norm ofxk(δ)and the residual norm kxk(δ)k2=β21 X
ωi>δ
(κ1i)2 and krk(δ)kr =β12X
ωi≤δ
κ21i can be obtained for fixedk andδwithout formingx (δ)ory (δ)explicitly.
Large problems
Tikhonov regularization for large problems
Consider the Tikhonov regularization in standard form (L=I) and its normal equations
(ATA+µ−1I)x =ATb (58) where the regularization parameterµis positive and finite.
Usuallyλ=µ−1is chosen, but for the method to develop (58) is more convenient.
The solution of (58) is
xµ= (ATA+µ−1I)−1ATb (59) and the discrepancy
dµ=b−Axµ. (60)
We assume that an estimate for the errorεofbis explicitly known. We seek for a parameterµˆsuch that
ε≤ kdµˆk ≤ηε (61)
where the choice ofηdepends on the accuracy of the estimateε.
TUHH Heinrich Voss Least Squares Problems Valencia 2010 75 / 82
Tikhonov regularization for large problems cnt.
Let
φ(µ) :=kb−Axµk2. (62)
Substituting the expression (59) into (62) and applying the identity I−A(ATA+µ−1I)−1AT = (µAAT +I)−1 one gets that
φ(µ) =bT(µAAT +I)−2b. (63)
Hence,φ0(µ)≤0 andφ00(µ)≥0 forµ≥0, i.e.φis monotonically decreasing and convex, and forτ∈(0,kbk)the equationφ(µ) =τhas a unique solution.
Notice, that the functionν7→φ(1/ν)cannot be guaranteed to be convex. This is the reason why the regularization parameterµis used instead ofλ=1/µ.
Large problems
Tikhonov regularization for large problems cnt.
Newton’s method converges globally to a solution ofφ(µ)−ε2=0. However, for large dimension the evaluation ofφ(µ)andφ0(µ)for fixedµis much too expensive (or even impossible due to the high condition number ofA).
Calvetti and Reichel (2003) took advantage of partial bidiagonalization ofAto construct bounds toφ(µ)and thus determine aµˆsatisfying (61).
Application of`≤min{m,n}bidiagonalization steps (cf. (48)) gives the decomposition
AV`=U`B`+β`+1u`+1eT`, ATU`=V`B`T, b=β1U`e1 (64) whereV`∈Rn×`andU`∈Rn×`have orthonormal columns andU`Tu`+1=0.
β`+1is a nonnegative scalar, andku`+1k=1 whenβ`+1>0.
B`∈R`×`is the bidiagonal matrix with diagonal elementsαj and (positive) subdiagonal elementsβj (B`is the submatrix ofB¯`in (49) containing its first` rows).
TUHH Heinrich Voss Least Squares Problems Valencia 2010 77 / 82
Tikhonov regularization for large problems cnt.
Combining the equations in (64) yields
AATU`=U`B`B`T +α`β`+1u`+1e`T (65) which shows that the columns ofU`are Lanczos vectors, and the matrix T`:=B`BT` is the tridiagonal matrix that would be obtained by applying` Lanczos steps to the symmetric matrixAAT with initial vectorb.
Introduce the functions
φ`(µ) :=kbk2e1T(µB`BT` +I`)−2e1, (66) φ¯`(µ) :=kbk2eT1(µB¯`B¯`T +I`+1)−2e1, (67)
Large problems
Tikhonov regularization for large problems cnt.
Substituting the spectral factorizationAAT =WΛWT with Λ =diag{λ1, . . . , λm}andWTW =Iintoφone obtains
φ(µ) =bTW(µΛ +I)−2WTb=
m
X
j=1
b˜j2 (µλj+1)2 =:
∞
Z
0
1
(µλ+1)2dω(λ) (68) where˜b=WTb, andωchosen to be a piecewise constant function with jump discontinuities of heightb˜2j at the eigenvaluesλj.
Thus, the integral in the right-hand side is a Stieltjes integral defined by the spectral factorization ofAAT and by the vectorb.
TUHH Heinrich Voss Least Squares Problems Valencia 2010 79 / 82
Tikhonov regularization for large problems cnt.
Golub and Meurant (1993) proved thatφ`defined in (66) is the`-point Gauß quadrature rule associated with the distribution functionωapplied to the function
ψ(t) := (µt+1)−2. (69) Similarly, the functionφ¯`in (67) is the(`+1)-point Gauß-Radau quadrature rule with an additional node at the origin associated with the distribution functionωapplied to the functionψ.
Since for any fixed positive value ofµthe derivatives ofψwith respect tot of odd order are strictly negative and the derivatives of even order are strictly positive fort≥0, the remainder terms for Gauß and Gauß-Radau quadrature rules show
φ`(µ)< φ(µ)<φ¯`(µ). (70) Instead of computing a valueµthat satisfies (61), one seeks to determine values`andµsuch that
Large problems
Tikhonov regularization for large problems cnt.
A method for computing a suitableµcan be based on the fact that for every µ >0 it holds that
φk(µ)< φ`(µ)< φ(µ)<φ¯`(µ)<φ¯k(µ) for 1≤k < ` (72) which was proved by Hanke (2003).
In general the value`in pairs(`, µ)that satisfy (71) can be chosen quite small (cf. the examples in the paper of Calvetti and Reichel).
Once a suitable valueµˆof the regularization parameter is available, the regularized solution
xµ,`ˆ =V`yµ,`ˆ (73) wheryµ,`ˆ ∈R`satisfies the Galerkin equation
V`T(ATA+ ˆµ−1I)V`yµ,`ˆ =V`ATb. (74)
TUHH Heinrich Voss Least Squares Problems Valencia 2010 81 / 82