[PPT] - Extrapolation techniques and multiparameter treatment for Tikhonov PowerPoint Presentation

SLIDE 1

Extrapolation techniques and multiparameter treatment for Tikhonov regularization A review

Michela Redivo Zaglia University of Padua - Italy Joint work with

SLIDE 2

AN APPLICATION TO REGULARIZATION

When a p × p system Ax = b is ill-conditioned, its solution cannot be computed accurately. Tikhonov’s regularization consists in computing the vector xλ which minimizes the quadratic functional J(λ, x) = Ax − b2 + λHx2

q × p (q ≤ p) matrix. · denotes the Euclidean norm. This vector xλ is the solution of the system (C + λE)xλ = AT b, where C = AT A and E = HT H.

SLIDE 3

If λ is close to zero, then xλ is badly computed while, if λ is far away from zero, xλ is well computed but the norm of the error x − xλ is quite large. For decreasing values of λ, the norm of the error x − xλ first decreases, and then increases when λ approaches 0. Thus the norm of the error, which is the sum of the theoretical error and the error due to the computer’s arithmetic, passes through a minimum corresponding to the optimal choice of the regularization parameter λ.

SLIDE 4

Several methods have been proposed to obtain an effective choice of λ (i.e. L-curve, Morozov discrepancy principle,

Idea: compute xλ for several values of λ, interpolate by a vector function of λ which mimics the exact behaviour, and then extrapolate at λ = 0. The main problem is to use a conveniently chosen vector function. For that purpose, let us study the exact behaviour of xλ with respect to λ.

SLIDE 5

We assume that H is a p × p non singular matrix, and we set y = Hx. Hence J(λ, x) = AH−1y − b2 + λy2 Using the SVD of AH−1 = UΣV T , it can be proved that xλ = H−1yλ with yλ =

σiγi σ2

and γi = (ui, b). So, we will choose an interpolation function of the same form but with a sum running only from i = 1 to k, with k < p. Several extrapolation methods based on this idea exist. They depend whether the vectors vi are known or not. We assume, without loss of generality, that H = I (i.e. xλ = yλ).

SLIDE 6

EXTRAPOLATION TECHNIQUES

RESTRICTED CASE: If the vectors v0, . . . , vk are known, we interpolate by the rational function Rk(λ) =

ai bi + λ vi, where the ai’s and the bi’s are 2k unknown scalars. For determining the parameters ai and bi, we impose the interpolation conditions (xn = xλn) xn =

ai bi + λn vi xn+1 =

ai bi + λn+1 vi.

SLIDE 7

Then, extrapolating at λ = 0 will give x ≃ Rk(0) =

aj bj vj = y(n)

We assume that vectors w1, . . . , wk such that (vi, wj) = δij are known, and we obtain y(n)

=

(xn, wj)(xn+1, wj) λn+1 − λn λn+1(xn+1, wj) − λn(xn, wj) vj. Increasing the value of k, for n fixed we also have y(n)

+ ak+1 bk+1 vk+1, k = 0, 1, . . . , p − 1.

SLIDE 8

If (vi, vj) = δij and wj = vj, the process is equivalent to the Truncated SVD (TSVD), that is y(n)

=

γi σi vi, γi = (ui, b). In this case, we can drop the superscript n since y(n)

does not depend on n. We set ek = x − yk and rk = b − Ayk. It holds ek+12 = ek2 − γ2

rk+12 = rk2 − γ2

SLIDE 9

We have also the following identities, ∀k, yk2 =

γ2

σ2

yk+12 = yk2 + γ2

σ2

ek2 =

γ2

σ2

(yk, ek) = (ek−1 − ek, ek) = 0 x2 = yk2 + ek2 rk2 =

γ2

SLIDE 10

EXTRAPOLATION TECHNIQUES

FULL CASE: If the vectors vi are unknown, we base the extrapolation process of a rational function of the form Rk(λ) =

1 bi + λ wi, k ≤ p, where the bi’s are unknown numbers and the wi are unknown vectors. They are determined, as before, by imposing that xn = Rk(λn) for some values of n. Then we extrapolate at the point λ = 0.

SLIDE 11

It corresponds to computing x ≃ yk = Rk(0) =

1 bi wi. Rk can be written as Rk(λ) = Pk−1(λ)/Qk(λ) with Pk−1(λ) = α0 + · · · + αk−1λk−1, αi ∈ Rp Qk(λ) =

(bi + λ) = β0 + · · · + βk−1λk−1 + λk, βi ∈ R. We gave 6 different algorithms for determining the two unknowns needed α0 and β0 (Rk(0) = α0/β0).

SLIDE 12

Let us describe one of them (the most satisfying one). We have to solve the interpolation problem Qk(λi)xi = Pk−1(λi), for i = 0, . . . , k − 1. Since Qk and Pk−1 are polynomials, we have, by Lagrange’s formula Qk(λ) =

Li(λ)Qk(λi) Pk−1(λ) =

=

SLIDE 13

with Li(λ) =

j=i

λ − λj λi − λj and

j=i

λ − λj λi − λj . Let λk = λj, for j = 0, . . . , k − 1. We have

Let s1, . . . , sp be linearly independent vectors. Setting ci = Qk(λi)/Qk(λk) and multiplying scalarly the preceding equation by sj, for j = 1, . . . , p, leads to the following linear system

j = 1, . . . , p. Solving this system in the least squares sense gives c0, . . . , ck−1.

SLIDE 14

Since the polynomial Qk(λ) = k

ck = 1, we have a supplementary condition which gives the value Qk(λk). Thus Qk(λi) = ciQk(λk), and, finally, β0 = Qk(0) is given by β0 =

Li(0)Qk(λi). From what precedes, we see that α0 = Pk−1(0) =

and it follows that yk = Rk(0) = Pk−1(0) Qk(0) = β0 α0 .

SLIDE 15

NUMERICAL EXAMPLES: A wide numerical experimentation has been performed. We used several kind of matrices A, heat, ilaplace, shaw, spikes

hilbert, lotkin, moler

different solutions xt (t means true solution), given

xi = 1 lin xi = i quad xi =

p

2 , i = 1, . . . , p sin xi = sin

various matrices H, I

D1, D2, D3

and we also try the case of a noised data vector.

SLIDE 16

The tests show the effectiveness of the procedures, but that the best approximation, denoted by xopt, depends on the values of λn chosen for interpolating and that the norm of the error xt − xoptcan be strongly influenced by the choice

SLIDE 17

MULTIPARAMETER REGULARIZATION

A good choice of the matrix H often depends on the mathematical properties of the solution. Using several regularization terms avoids this difficult choice. We are looking for the vector xλ which minimizes the quadratic functional J(λ, x) = k

λiHix2

with λ = (λ1, . . . , λk)T .

SLIDE 18

It is also the solution of the system

λiEi

with C = AT A and Ei = HT

We have the following relation between xλ and x

λiC−1Ei

We note that xλ is also the vector minimizing J(λ, x) =

.

SLIDE 19

Hence, we consider the k vectors xλi solving the minimization problems min

, i = 1, . . . , k. These vectors satisfy the linear systems (C + kλiEi)xλi = AT b, i = 1, . . . , k. Thus we compute k one-parameter regularized solutions, and, after, we consider the approximation of xλ given by the linear combination

αixλi, where α = (α1, . . . , αk)T and k

SLIDE 20

How to choose the vector α? The following relation between xλ and xλ(α) holds xλ = xλ(α) +

λiEi

ρλ(α), where ρλ(α) =

αi  kλiEi −

λjEj   xλi. The vector α is chosen to minimize ρλ(α). It is the solution of an overdetermined system which is solved in the least-squares sense.

SLIDE 21

How to estimate the vector λ? The parameters λ1, . . . , λk can be chosen according to a test based on a modification of the Generalized Cross Validation. It consists in minimizing the function Z(λ) = k

Vi(λi)

with Vi(λ) = Ni(λ) Di(λ), i = 1, . . . , k, and Ni(λ) =  

(γ(i)

2 

, Di(λ) =

kλ (γ(i)

. The c(i)

decompositions.

SLIDE 22

NUMERICAL EXAMPLES: In many tests performed, regularizing with more than one parameter seems to increase the possibilities of computing a good approximation to the solution of an ill-conditioned linear system. Moreover, in general, the error is never worse than the worst

The algorithm seems enough stable, robust and easy to implement (only the solution of k one-parameter regularization problems and one additional linear system). In the sequel some tables reporting the relative errors x − x x corresponding to the solution obtained with the best λ.

SLIDE 23

σ = 0 I D1 D2 MP heat

8.7 · 10−6 1.4 · 10−14 2.1 · 10−15 6.9 · 10−16 lin 3.7 · 10−1 1.8 · 10−2 7.7 · 10−15 2.3 · 10−15 quad 2.8 · 10−5 2.8 · 10−5 9.7 · 10−3 9.4 · 10−3 sin 9.6 · 10−2 6.7 · 10−6 5.0 · 10−7 5.9 · 10−10 ilaplace

6.1 · 10−1 5.3 · 10−15 1.7 · 10−15 9.8 · 10−16 lin 8.4 · 10−1 2.8 · 10−1 1.1 · 10−15 2.8 · 10−15 quad 7.9 · 10−1 7.2 · 10−1 6.7 · 10−1 5.6 · 10−1 sin 7.1 · 10−1 5.2 · 10−1 2.1 · 10−1 3.7 · 10−1 lotkin lin 1.1 · 10−6 1.1 · 10−6 1.6 · 10−13 1.7 · 10−15 quad 4.5 · 10−6 4.5 · 10−6 1.1 · 10−3 2.4 · 10−4 sin 2.2 · 10−4 2.2 · 10−4 2.1 · 10−3 5.5 · 10−4 moler

1.6 · 10−8 3.2 · 10−14 2.1 · 10−14 1.9 · 10−14 quad 5.9 · 10−8 8.5 · 10−4 1.3 · 10−3 9.4 · 10−11 sin 5.7 · 10−4 3.3 · 10−4 2.4 · 10−4 1.3 · 10−7

SLIDE 24

σ = 10−6 I D1 D2 MP heat

1.1 · 10−4 1.1 · 10−4 1.7 · 10−6 5.3 · 10−8 lin 1.9 · 10−4 1.9 · 10−4 2.8 · 10−6 2.8 · 10−6 quad 2.5 · 10−4 2.5 · 10−4 9.7 · 10−3 9.4 · 10−3 sin 9.6 · 10−2 8.7 · 10−2 1.8 · 10−2 1.6 · 10−4 ilaplace

7.2 · 10−1 1.7 · 10−5 4.4 · 10−7 5.7 · 10−7 lin 9.4 · 10−1 4.2 · 10−1 7.3 · 10−7 7.3 · 10−7 quad 7.9 · 10−1 6.0 · 10−1 1.1 1.3 · 10−1 sin 7.1 · 10−1 4.1 · 10−1 5.1 · 10−1 3.4 · 10−1 lotkin lin 2.8 · 10−3 2.8 · 10−3 3.4 · 10−7 3.4 · 10−7 quad 4.1 · 10−2 4.1 · 10−2 2.3 · 10−2 2.3 · 10−2 sin 7.3 · 10−2 7.3 · 10−2 4.7 · 10−2 4.7 · 10−2 moler

6.4 · 10−6 3.4 · 10−7 2.0 · 10−8 2.7 · 10−10 quad 3.5 · 10−6 3.5 · 10−6 8.9 · 10−7 7.6 · 10−7 sin 6.6 · 10−7 2.9 · 10−6 1.1 · 10−1 4.9 · 10−7

SLIDE 25

FUTURE WORKS

SLIDE 26

References

Extrapolation techniques for ill–conditioned linear systems

Multi-parameter regularization techniques for ill-conditioned linear systems