Extrapolation techniques and multiparameter treatment for Tikhonov regularization A review Michela Redivo Zaglia University of Padua - Italy Joint work with C. Brezinski (Lille-France) G. Rodriguez, S. Seatzu (Cagliari-Italy) ➜ Tikhonov regularization ➜ Extrapolation techniques ➜ Multiparameter regularization 1
A N 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 − b � 2 + λ � Hx � 2 over all vectors x , where λ ≥ 0 is a parameter, and H a given q × p ( q ≤ p ) matrix. � · � denotes the Euclidean norm. This vector x λ is the solution of the system ( C + λE ) x λ = A T b, where C = A T A and E = H T H . A N APPLICATION TO REGULARIZATION 2
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 λ . A N APPLICATION TO REGULARIZATION 3
Several methods have been proposed to obtain an effective choice of λ (i.e. L -curve, Morozov discrepancy principle, GCV , . . . ). But each of these methods can fail. 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 λ . A N APPLICATION TO REGULARIZATION 4
We assume that H is a p × p non singular matrix, and we set y = Hx . Hence J ( λ, x ) = � AH − 1 y − b � 2 + λ � y � 2 Using the SVD of AH − 1 = U Σ V T , it can be proved that x λ = H − 1 y λ with p � σ i γ i y λ = i + λv i , and γ i = ( u i , b ) . σ 2 i =1 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 v i are known or not . We assume, without loss of generality, that H = I (i.e. x λ = y λ ). A N APPLICATION TO REGULARIZATION 5
E XTRAPOLATION TECHNIQUES R ESTRICTED CASE : If the vectors v 0 , . . . , v k are known , we interpolate by the rational function k � a i R k ( λ ) = b i + λ v i , i =1 where the a i ’s and the b i ’s are 2 k unknown scalars . For determining the parameters a i and b i , we impose the interpolation conditions ( x n = x λ n ) k � a i x n = v i b i + λ n i =1 k � a i x n +1 = v i . b i + λ n +1 i =1 E XTRAPOLATION TECHNIQUES 6
Then, extrapolating at λ = 0 will give k � a j v j = y ( n ) x ≃ R k (0) = k . b j j =1 We assume that vectors w 1 , . . . , w k such that ( v i , w j ) = δ ij are known , and we obtain k � λ n +1 − λ n y ( n ) = ( x n , w j )( x n +1 , w j ) λ n +1 ( x n +1 , w j ) − λ n ( x n , w j ) v j . k j =1 Increasing the value of k , for n fixed we also have + a k +1 y ( n ) k +1 = y ( n ) v k +1 , k = 0 , 1 , . . . , p − 1 . k b k +1 E XTRAPOLATION TECHNIQUES 7
If ( v i , v j ) = δ ij and w j = v j , the process is equivalent to the Truncated SVD (TSVD) , that is k � γ i y ( n ) = v i , γ i = ( u i , b ) . k σ i i =1 In this case, we can drop the superscript n since y ( n ) does not k depend on n . We set e k = x − y k and r k = b − Ay k . It holds � e k � 2 − γ 2 � e k +1 � 2 k +1 /σ 2 = k +1 � r k � 2 − γ 2 � r k +1 � 2 = k +1 . E XTRAPOLATION TECHNIQUES 8
We have also the following identities, ∀ k , k � γ 2 � y k � 2 i = σ 2 i i =1 � y k � 2 + γ 2 k +1 � y k +1 � 2 = σ 2 k +1 p � γ 2 � e k � 2 i = σ 2 i i = k +1 ( y k , e k ) = ( e k − 1 − e k , e k ) = 0 � y k � 2 + � e k � 2 � x � 2 = p � � r k � 2 γ 2 = i . i = k +1 E XTRAPOLATION TECHNIQUES 9
E XTRAPOLATION TECHNIQUES F ULL CASE : If the vectors v i are unknown , we base the extrapolation process of a rational function of the form k � 1 R k ( λ ) = b i + λ w i , k ≤ p, i =1 where the b i ’s are unknown numbers and the w i are unknown vectors . They are determined, as before, by imposing that x n = R k ( λ n ) for some values of n . Then we extrapolate at the point λ = 0 . E XTRAPOLATION TECHNIQUES 10
It corresponds to computing k � 1 x ≃ y k = R k (0) = w i . b i i =1 R k can be written as R k ( λ ) = P k − 1 ( λ ) /Q k ( λ ) with α 0 + · · · + α k − 1 λ k − 1 , α i ∈ R p P k − 1 ( λ ) = k � ( b i + λ ) = β 0 + · · · + β k − 1 λ k − 1 + λ k , Q k ( λ ) = β i ∈ R . i =1 We gave 6 different algorithms for determining the two unknowns needed α 0 and β 0 ( R k (0) = α 0 /β 0 ). E XTRAPOLATION TECHNIQUES 11
Let us describe one of them (the most satisfying one). We have to solve the interpolation problem Q k ( λ i ) x i = P k − 1 ( λ i ) , for i = 0 , . . . , k − 1 . Since Q k and P k − 1 are polynomials, we have, by Lagrange’s formula k � Q k ( λ ) = L i ( λ ) Q k ( λ i ) i =0 k − 1 � � P k − 1 ( λ ) = L i ( λ ) P k − 1 ( λ i ) i =0 k − 1 � � = L i ( λ ) Q k ( λ i ) x i i =0 E XTRAPOLATION TECHNIQUES 12
with k k − 1 � � λ − λ j λ − λ j � L i ( λ ) = and L i ( λ ) = . λ i − λ j λ i − λ j j =0 j =0 j � = i j � = i Let λ k � = λ j , for j = 0 , . . . , k − 1 . We have k − 1 � � L i ( λ k ) Q k ( λ i ) x i = Q k ( λ k ) x k . i =0 Let s 1 , . . . , s p be linearly independent vectors . Setting c i = Q k ( λ i ) /Q k ( λ k ) and multiplying scalarly the preceding equation by s j , for j = 1 , . . . , p , leads to the following linear system k − 1 � � L i ( λ k )( x i , s j ) c i = ( x k , s j ) , j = 1 , . . . , p. i =0 Solving this system in the least squares sense gives c 0 , . . . , c k − 1 . E XTRAPOLATION TECHNIQUES 13
Since the polynomial Q k ( λ ) = � k i =0 L i ( λ ) Q k ( λ i ) is monic and c k = 1 , we have a supplementary condition which gives the value Q k ( λ k ) . Thus Q k ( λ i ) = c i Q k ( λ k ) , and, finally, β 0 = Q k (0) is given by k � β 0 = L i (0) Q k ( λ i ) . i =0 From what precedes, we see that k − 1 � � α 0 = P k − 1 (0) = L i (0) Q k ( λ i ) x i i =0 and it follows that y k = R k (0) = P k − 1 (0) = β 0 . Q k (0) α 0 E XTRAPOLATION TECHNIQUES 14
N UMERICAL EXAMPLES : A wide numerical experimentation has been performed. We used several kind of matrices A , heat, ilaplace, shaw, spikes Hansen Matlab Regularization Toolbox hilbert, lotkin, moler Matlab (gallery function) different solutions x t ( t means true solution ), given defined as in the Regularization Toolbox ones x i = 1 lin x i = i � � p �� 2 , quad x i = i − i = 1 , . . . , p � � 2 sin x i = sin 2 π ( i − 1) /p i = 1 , . . . , p various matrices H , I Identity matrix D 1 , D 2 , D 3 discrete approximations of the first, second and third derivative and we also try the case of a noised data vector. E XTRAPOLATION TECHNIQUES 15
The tests show the effectiveness of the procedures , but that the best approximation , denoted by x opt , depends on the values of λ n chosen for interpolating and that the norm of the error � x t − x opt � can be strongly influenced by the choice of the regularizating matrix H . E XTRAPOLATION TECHNIQUES 16
M ULTIPARAMETER 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 � � k � � Ax − b � 2 + λ i � H i x � 2 J ( λ, x ) = k , i =1 with λ = ( λ 1 , . . . , λ k ) T . M ULTIPARAMETER REGULARIZATION 17
It is also the solution of the system � � k � x λ = A T b, C + λ i E i i =1 with C = A T A and E i = H T i H i . We have the following relation between x λ and x � � k � λ i C − 1 E i I + x λ = x. i =1 We note that x λ is also the vector minimizing k � � � Ax − b � 2 + kλ i � H i x � 2 � J ( λ, x ) = . i =1 M ULTIPARAMETER REGULARIZATION 18
Hence, we consider the k vectors x λ i solving the minimization problems � � Ax − b � 2 + kλ i � H i x � 2 � min , i = 1 , . . . , k. x These vectors satisfy the linear systems ( C + kλ i E i ) x λ i = A T 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 k � x λ ( α ) = � α i x λ i , i =1 where α = ( α 1 , . . . , α k ) T and � k i =1 α i = 1 . M ULTIPARAMETER REGULARIZATION 19
How to choose the vector α ? The following relation between x λ and � x λ ( α ) holds � � − 1 k � x λ = � x λ ( α ) + C + λ i E i ρ λ ( α ) , i =1 where k k � � kλ i E i − x λ i . ρ λ ( α ) = α i λ j E j i =1 j =1 The vector α is chosen to minimize ρ λ ( α ) . It is the solution of an overdetermined system which is solved in the least-squares sense . M ULTIPARAMETER REGULARIZATION 20
Recommend
More recommend
