Non-stationary structure-preserving preconditioning for image restoration Pietro Dell’Acqua joint work with Marco Donatelli, Lothar Reichel Computational Methods for Inverse Problems in Imaging 16-18 July 2018, Como Pietro Dell’Acqua 1 / 33
The blurring model We are concerned with the restoration of blurred and noise-corrupted images in two space-dimensions. The blurring is modeled by a convolution and the image degradation model is of the form � ① ∈ Ω ⊂ R 2 , g ( ① ) = [ Kf ]( ① ) + ν ( ① ) = R 2 h ( ① − ② ) f ( ② ) d ② + ν ( ① ) , where f represents the (desired but unavailable) exact image, h the space invariant point-spread function (PSF) with compact support, ν random noise, and g the (available) blurred and noise-corrupted image. Hence, f and g are real-valued non-negative functions that determine the light intensity of the desired and available images, respectively. Pietro Dell’Acqua 2 / 33
Discretization Discretization of the integral equation at equidistant nodes gives the linear system of algebraic equations � i ∈ Z 2 . g i = h i − j f j + ν i , j ∈ Z 2 The entries of the discrete images ❣ = [ g i ] and ❢ = [ f j ] represent the light intensity at each picture element (pixel) and ν = [ ν i ] models the noise-contamination at these pixels. The pixels with index i ∈ [1 , n ] 2 make up the finite field of view (FOV), which for notational simplicity is assumed to be square. We would like to determine an accurate approximation of the exact image ❢ in the FOV given ❤ = [ h i ], distributional information about ν , and the blurred image ❣ in the FOV. Pietro Dell’Acqua 3 / 33
Boundary conditions A common approach is to impose boundary conditions (BC) on the image to obtain the next linear system A ∈ R n 2 × n 2 , ❢ , ❣ ∈ R n 2 . A ❢ = ❣ , The boundary conditions specify that the values at pixels outside the FOV are linear combinations of values at certain pixels inside the FOV. Popular boundary conditions include: Zero Dirichlet boundary conditions (ZDBC) Periodic boundary conditions (PBC) Reflective boundary conditions (RBC) Anti-Reflective boundary conditions (ARBC) Pietro Dell’Acqua 4 / 33
Boundary conditions Pietro Dell’Acqua 5 / 33
Goals Focus Development of a fast and stable iterative regularization method for the approximate solution of when A is defined by a non-symmetric PSF h with ARBC. GMRES GMRES method is commonly used for the iterative solution of linear systems of equations with a square non-symmetric matrix that is obtained by the discretization of a well-conditioned problem. In our context, preconditioners are employed to accelerate the convergence. Low level of noise We remark that available preconditioning strategies for GMRES can determine accurate restorations, but may require many iterations when the noise level is low. Pietro Dell’Acqua 6 / 33
Flexible GMRES We investigate non-stationary preconditioning with GMRES-type methods with the aim to obtain accurate restorations within few iterations also in presence of a strongly non-symmetric PSF and a small noise level. We use the flexible GMRES (F-GMRES) method to solve, at step k , ③ = P − 1 AP k ③ = ❣ , k ❢ . Thus, the preconditioner P k is modified in each iteration. By using a suitable sequence of preconditioners P k , we obtain a preconditioned F-GMRES method that is well suited for image restoration. The preconditioners P k depend on a thresholding parameter α k that, in turn, depends on the amount of noise in ❣ . Pietro Dell’Acqua 7 / 33
Preconditioning We first discuss left-preconditioned Landweber iteration. ZA ❢ = Z ❣ . Application of Landweber iteration yields the iterates ❢ k +1 = ❢ k + Z ( ❣ − A ❢ k ) . We can determine the preconditioner Z ∈ R n 2 × n 2 by Tikhonov filter λ i , j ˘ λ i , j = , i , j = 0 , 1 , . . . , n − 1 , | λ i , j | 2 + α where α > 0 is a regularization parameter and the bar denotes complex conjugation. The BCCB matrix Z can be defined analogously for other BC. Pietro Dell’Acqua 8 / 33
Preconditioning Circulant non-stationary version of the iterative method circ = C T ( CC T + α k I ) − 1 , ❢ k +1 = ❢ k + Z k Z k circ r k , r k = ❣ − A ❢ k . Here C is the BCCB matrix associated with the PSF that defines the matrix A . Structured non-stationary version of the iterative method struct = B ( C T ( CC T + α k I ) − 1 ) , ❢ k +1 = ❢ k + Z k Z k struct r k , r k = ❣ − A ❢ k , where the operator B denotes the application of BC to the circulant matrix C T ( CC T + α k I ) − 1 ) and, thus, modifying the structure of the latter. The matrix Z k struct may be considered a preconditioner. In particular, we may solve the right-preconditioned linear system struct = B ( C T ( CC T + α k I ) − 1 ) , P k = Z k by F-GMRES. The parameter α k allows the preconditioner P k to be varied during the iterations. Pietro Dell’Acqua 9 / 33
Preconditioning Summarizing, the mask of the Fourier coefficients of the preconditioner P k can be determined by carrying out the following steps: 1 Compute λ i , j by the FFT applied to the mask associated with the PSF. 2 Compute ˘ λ i , j for α = α k . 3 Compute ˘ H by the IFFT applied to ˘ λ i , j . In actual computations, we modify the BCCB matrix C T ( CC T + α I ) − 1 to correspond to ARBC. This yields a structured preconditioner B ( C T ( CC T + α I ) − 1 ), where B is an operator that imposes the ARBC. We remark that P k is not explicitly formed, but only matrix-vector products with P k are evaluated. Pietro Dell’Acqua 10 / 33
Flexible GMRES algorithm F-GMRES is a minimal residual method that is designed for application with a sequence of preconditioners. Given a set of linearly independent vectors ✉ 1 , ✉ 2 , . . . , ✉ ℓ ∈ R n 2 , the method determines the decomposition AU ℓ = V ℓ +1 H ℓ +1 ,ℓ , where V ℓ +1 = [ ✈ 1 , ✈ 2 , . . . , ✈ ℓ +1 ] ∈ R n 2 × ( ℓ +1) has orthonormal columns with ✈ 1 = ❣ / � ❣ � , U ℓ = [ ✉ 1 , ✉ 2 , . . . , ✉ ℓ ] ∈ R n 2 × ℓ and H ℓ +1 ,ℓ = [ h i , j ] ∈ ∈ R ( ℓ +1) × ℓ is of upper Hessenberg form. Let ❡ 1 = [1 , 0 , . . . , 0] T ∈ R ℓ +1 denote the first axis vector. Then ✇ ∈ range ( U ℓ ) � A ✇ − ❣ � = min min ② ∈ R ℓ � AU ℓ ② − ❣ � = min ② ∈ R ℓ � H ℓ +1 ,ℓ ② − ❡ 1 � ❣ � � . Pietro Dell’Acqua 11 / 33
Flexible GMRES algorithm 1. ✈ 1 = ❣ / � ❣ � 2. for k = 1 , 2 , . . . , ℓ do 3. ✉ k = P k ✈ k ; ✈ = A ✉ k 4. for i = 1 , 2 , . . . , k do h i , k = ✈ T ✈ i ; ✈ = ✈ − h i , k ✈ i 5. 6. end 7. h k +1 , k = � ✈ � ; ✈ k +1 = ✈ / h k +1 , k 8. end 9. define U ℓ = [ ✉ 1 , ✉ 2 , . . . , ✉ ℓ ] ∈ R n 2 × ℓ 10. define H ℓ +1 ,ℓ = [ h i , j ] ∈ R ( ℓ +1) × ℓ upper Hessenberg 11. compute ② ℓ := arg min ② ∈ R ℓ � H ℓ +1 ,ℓ ② − ❡ 1 � ❣ � � and ❢ ℓ = U ℓ ② ℓ Pietro Dell’Acqua 12 / 33
Discrepancy principle We assume that a fairly accurate bound ε for the norm of the error ν in the vector ❣ is available. Thus, � ν � ≤ ε. Then we can apply the discrepancy principle to decide when to stop the iterations. Let ❢ 1 , ❢ 2 , ❢ 3 , . . . be a sequence of approximate solutions determined by the algorithm and define the associated residual vectors r ℓ = ❣ − A ❢ ℓ , ℓ = 1 , 2 , 3 , . . . . The discrepancy principle prescribes that the iterations be terminated as soon as a residual vector r ℓ that satisfies � r ℓ � ≤ ε has been determined. Pietro Dell’Acqua 13 / 33
Geometric sequence Regarding the determination of the parameters α k of the preconditioners P k , a well-established choice in the context of iterated Tikhonov methods is the geometric sequence α k = α 0 q k − 1 , k = 1 , 2 , . . . , where α 0 > 0 and 0 < q < 1. The value of the initial regularization parameter α 0 is not critical as long as it is not too small. In our numerical experiments, we set α 0 = 1 and q = 0 . 8. Pietro Dell’Acqua 14 / 33
DH sequence In a paper by Donatelli and Hanke, the regularization parameter α k is obtained at each step k by solving with a few steps of Newton’s method the next non-linear equation � r k − CZ k circ r k � = q k � r k � , 0 < q k < 1 , where � · � denotes the Euclidean norm. Here the parameter q k depends on the noise level and it is related to a value 0 < ρ circ < 1 / 2 which should be as small as possible and satisfies the assumption ∀ ③ ∈ R n . � ( C − A ) ③ � ≤ ρ circ � A ③ � , The parameter ρ circ needs to be set in the algorithm. It is associated with the relative distance between A and C , so if A is well approximated by its BCCB counterpart C , then ρ circ can be chosen as a very small value. Pietro Dell’Acqua 15 / 33
DH sequence Recently the same approach has been applied to the structured case and so the goal is to estimate α k by solving � r k − AZ k struct r k � = q k � r k � , which is not computationally practicable. Thus, the regularization parameter α k is estimated by using equation of circulant case, that can be seen as a computable approximation of the desired one. The iterations are stopped under a special choice of the discrepancy principle, that is � r ℓ � ≤ 1 + 2 ρ struct ε. 1 − 2 ρ struct In the examples presented here, we set ρ struct = 10 − 2 for Test 1 and ρ struct = 10 − 1 for Test 2. Pietro Dell’Acqua 16 / 33
Recommend
More recommend
