Numerical Methods for Ill-Posed Problems I Lothar Reichel Summer School on Applied Analysis TU Chemnitz, Oct. 4-8, 2010
Outline of Lecture 1: • Inverse and ill-posed problems • The singular value decomposition • Tikhonov regularization
Inverse problems • Inverse problems arise when one seeks to determine the cause of an observed effect. – Helioseismology: Determine the structure of the sun by measurements from earth or space. – Medical imaging, e.g., electrocardiographic imaging, computerized tomography. – Image restoration: Determine the unavailable exact image from an available contaminated version. • Inverse problems often are ill-posed.
Ill-posed problems A problem is said to be ill-posed if it has at least one of the properties: • the problem does not have a solution, • the problem does not have a unique solution, • the solution does not depend continuously on the data.
Linear discrete ill-posed problems Ax = b arise from the discretization of linear ill-posed problems (Fredholm integral equations of the 1st kid) or, naturally, in discrete form (image restoration). • The matrix A is of ill-determined rank, possibly singular. System may be inconsistent. • The right-hand side b represents available data that generally is contaminated by an error.
Available contaminated, possibly inconsistents, linear system Ax = b (1) Unavailable associated consistent linear system with error-free right-hand side Ax = ˆ b (2) Let ˆ x denote the desired solution of (2), e.g., the minimal-norm solution. Task: Determine an approximate solution of (1) that is a good approximation of ˆ x .
Define the error e := ˆ b − b noise and let ǫ := � e � The choice of solution method depends on whether an estimate of ǫ is available. How much damage can a little noise really do? How much noise requires the use of special techniques?
Example 1: Fredholm integral equation of the 1st kind � π 0 exp( − st ) x ( t ) dt = 2sinh( s ) 0 ≤ s ≤ π , 2 . s Determine solution x ( t ) = sin( t ). Discretize integral by Galerkin method using piecewise constant functions. Code baart from Regularization Tools.
This gives a linear system of equations Ax = ˆ A ∈ R 200 × 200 , ˆ b ∈ R 200 . b, A is numerically singular. Let the “noise” vector e in b have normally distributed entries with mean zero and ǫ = � e � = 10 − 3 � b � b := ˆ b + e i.e., 0.1% relative noise
right−hand side 0.26 0.25 0.24 no added noise added noise/signal=1e−3 0.23 0.22 0.21 0.2 0.19 0.18 0.17 0 20 40 60 80 100 120 140 160 180 200
Solution Ax=b: noise level 10 −3 13 x 10 1.5 1 0.5 0 −0.5 −1 −1.5 0 20 40 60 80 100 120 140 160 180 200
solution of Ax=b : no added noise 150 100 50 0 −50 −100 0 20 40 60 80 100 120 140 160 180 200
exact solution 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 20 40 60 80 100 120 140 160 180 200
The singular value decomposition
The SVD of the m × n matrix A , m ≥ n : A = U Σ V T U = [ u 1 , u 2 , . . . , u m ] orthogonal , m × m, V = [ v 1 , v 2 , . . . , v n ] orthogonal , n × n, Σ = diag[ σ 1 , σ 2 , . . . , σ n ] , m × n with σ 1 ≥ σ 2 ≥ . . . ≥ σ n ≥ 0 .
Application: Least-squares approximation Let the matrix A ∈ R m × n , m ≥ n , represent the “model” and the vector b ∈ R m the data. Solve x � Ax − b � 2 = min x � U Σ V T x − b � 2 = min x � Σ V T x − U T b � 2 . min Let y = [ y 1 , y 2 , . . . , y n ] T = V T x and b ′ = [ b ′ m ] T = U T b . Then 1 , b ′ 2 , . . . , b ′ n m x � Ax − b � 2 = min � Σ y − b ′ � 2 = ( σ j y j − b ′ j ) 2 + ( b ′ j ) 2 . � � min y j =1 j = n +1
If A is of full rank, then all σ j > 0 and y j = b ′ j , 1 ≤ j ≤ n, σ j yields the solution x = V y. If some σ j = 0, then least-squares solution not unique. Often one is interested in the least-squares solution of minimal norm. Arbitrary components y j are set to zero.
Assume that σ 1 ≥ σ 2 . . . ≥ σ ℓ > σ ℓ +1 = . . . = σ n = 0 . Then A is of rank ℓ . Introduce the diagonal matrix Σ † = diag[1 /σ 1 , 1 /σ 2 , . . . , 1 /σ ℓ , 0 , . . . , 0] , n × m. The matrix A † = V Σ † U T is known as the Moore-Penrose pseudoinverse of A .
The solution of the least-squares problem min x � Ax − b � of minimal Eucliden norm can be expressed as x = A † b. Moreover, AA † = P R ( A ) . A † A = I,
Note n A = U Σ V T = σ j u j v T � j . j =1 Define k σ j u j v T � A k := j , 1 ≤ k ≤ ℓ. j =1 Then A k is of rank k ; A k is the sum of k rank-one matrices σ j u j v T j . Moreover, � A − A k � = rank( B ) ≤ k � A − B � = σ k +1 , min i.e., A k is the closest matrix of rank ≤ k to A .
Let b = ˆ b + e , where e denotes an error. Then ℓ u T j b x := A † b � = v j σ j j =1 j ˆ u T u T ℓ ℓ b j e � � = v j + v j σ j σ j j =1 j =1 ℓ u T j e � = x + ˆ v j . σ j j =1 If σ ℓ > 0 tiny, then u T ℓ e σ ℓ might be huge and x a meaningless approximation of ˆ x .
Recall k σ j u j v T � A k = j j =1 best rank- k approximation of A . Pseudoinverse of A k : k A † σ − 1 j v j u T � k := j , σ k > 0 j =1
Approximate ˆ x by A † x k := k b k u T j b � = v j σ j j =1 j ˆ u T u T k k b j e � � = v j + v j . σ j σ j j =1 j =1 for some k ≤ ℓ . How to choose k ?
Example 1 cont’d: Singular values of A 5 10 0 10 −5 10 −10 10 −15 10 −20 10 0 20 40 60 80 100 120 140 160 180 200
Example 1 cont’d: Right-hand side without noise Norm of computed approximate solution 2.8 2.6 2.4 2.2 2 ||x k || 1.8 1.6 1.4 1.2 1 0.8 0 5 10 15 20 25 30 35 40 45 50 k
Example 1 cont’d: Right-hand side without noise Norm of residual vector 0 −2 −4 log 10 ||b−Ax k || −6 −8 −10 −12 −14 0 5 10 15 20 25 30 35 40 45 50 k
Example 1 cont’d: Right-hand side without noise L−curve 0 −2 −4 log 10 ||b−Ax k || −6 −8 −10 −12 −14 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 log 10 ||x k ||
Example 1 cont’d: Right-hand side without noise: Blow-up L−curve −12.4 −12.6 −12.8 log 10 ||b−Ax k || −13 −13.2 −13.4 −13.6 −13.8 −14 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 log 10 ||x k ||
Example 1 cont’d: Exact and computed solutions 1.4 exact x 9 x 10 x 11 1.2 1 0.8 0.6 0.4 0.2 0 0 20 40 60 80 100 120 140 160 180 200
Example 1 cont’d: Right-hand side with relative noise 10 − 3 Norm of computed approximate solution 1.5 1.4 1.3 1.2 ||x k || 1.1 1 0.9 0.8 0 2 4 6 8 10 12 14 16 18 20 k
Example 1 cont’d: Right-hand side with relative noise 10 − 3 Norm of residual 0 10 −1 10 log 10 ||b−Ax k || −2 10 −3 10 −4 10 0 2 4 6 8 10 12 14 16 18 20 k
Example 1 cont’d: Right-hand side with relative noise 10 − 3 L−curve 0 −0.5 −1 log 10 ||b−Ax k || −1.5 −2 −2.5 −3 −3.5 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 log 10 ||x k ||
Example 1 cont’d: Right-hand side with relative noise 10 − 3 1.4 exact x 3 x 4 1.2 x 5 1 0.8 0.6 0.4 0.2 0 0 20 40 60 80 100 120 140 160 180 200
The SVD in Hilbert space
A : X → Y compact linear operator X , Y Hilbert spaces, �· , ·� inner products in X and Y , � · � associated norms. Compute minimal-norm solution ˆ x ∈ X of Ax = ˆ y, y ∈ R ( A ) , ˆ by the SVD of A .
Singular triplets { σ j , u j , v j } ∞ j =1 of A : σ j singular value, σ 1 ≥ σ 2 ≥ σ 3 ≥ . . . > 0 , u j ∈ Y left singular function, v j ∈ X right singular function, satisfy 1 , j = ℓ, � u j , u ℓ � = � v j , v ℓ � = 0 , j � = ℓ, and
A ∗ u j = σ j v j , Av j = σ j u j , j = 1 , 2 , 3 , . . . , ∞ � Ax = σ j � x, v j � u j , ∀ x ∈ X , j =1 ∞ A ∗ y � = σ j � y, u j � v j , ∀ y ∈ Y , j =1 where A ∗ : Y → X the adjoint of A . A compact ⇒ σ j cluster at zero.
Minimal-norm solution ∞ � ˆ y, u j � � x = ˆ v j . σ j j =1 x ∈ X implies Picard condition: ˆ ∞ y, u j �| 2 |� ˆ � < ∞ . σ 2 j j =1 Fourier coefficients c j = � ˆ ˆ y, u j � , j = 1 , 2 , 3 , . . . , have to converge to zero rapidly.
y = ˆ y + e available contaminated rhs Determine approximation of ˆ x by TSVD: k � y, u j � � x k = v j . σ j j =1 ˆ k ≥ 1 smallest index, such that � x ˆ k − ˆ x � = min k ≥ 1 � x k − ˆ x � .
Assume norm of error in rhs known: δ = � y − ˆ y � . z is said to satisfy the discrepancy principle if � Az − y � ≤ τδ, where τ > 1 constant independent of δ . The discrepancy principle selects the smallest index k = k δ , such that � Ax k δ − y � ≤ τδ.
Then • � Ax k − y � decreases monotonically as k increases. • k δ increases monotonically as δ ց 0. • x k δ → ˆ x as δ ց 0.
Tikhonov regularization
Recommend
More recommend
