Non-stationary structure-preserving preconditioning for image - - PowerPoint PPT Presentation

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 g(①) = [Kf ](①) + ν(①) =

• R2 h(① − ②)f (②)d② + ν(①),

① ∈ Ω ⊂ R2, 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 gi =

• j∈Z2

hi−jfj + νi, i ∈ Z2. The entries of the discrete images ❣ = [gi] and ❢ = [fj] 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 ❤ = [hi], 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❢ = ❣, A ∈ Rn2×n2, ❢ , ❣ ∈ Rn2. 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, APk③ = ❣, ③ = P−1

k ❢ .

Thus, the preconditioner Pk is modified in each iteration. By using a suitable sequence of preconditioners Pk, we obtain a preconditioned F-GMRES method that is well suited for image restoration. The preconditioners Pk 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 ∈ Rn2×n2 by Tikhonov filter ˘ λi,j = λi,j |λi,j|2 + α , i, j = 0, 1, . . . , n − 1, 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 ❢ k+1 = ❢ k + Z k

circr k,

Z k

circ = C T(CC T + αkI)−1,

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 ❢ k+1 = ❢ k +Z k

structr k,

Z k

struct = B(C T(CC T +αkI)−1),

r k = ❣ −A❢ k, where the operator B denotes the application of BC to the circulant matrix C T(CC T + αkI)−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 Pk = Z k

struct = B(C T(CC T + αkI)−1),

by F-GMRES. The parameter αk allows the preconditioner Pk to be varied during the iterations.

Pietro Dell’Acqua 9 / 33

Preconditioning

Summarizing, the mask of the Fourier coefficients of the preconditioner Pk 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 Pk is not explicitly formed, but only matrix-vector products with Pk 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, . . . , ✉ℓ ∈ Rn2, the method determines the decomposition AUℓ = Vℓ+1Hℓ+1,ℓ, where Vℓ+1 = [✈ 1, ✈ 2, . . . , ✈ ℓ+1] ∈ Rn2×(ℓ+1) has orthonormal columns with ✈ 1 = ❣/❣, Uℓ = [✉1, ✉2, . . . , ✉ℓ] ∈ Rn2×ℓ and Hℓ+1,ℓ = [hi,j] ∈ ∈ R(ℓ+1)×ℓ is of upper Hessenberg form. Let ❡1 = [1, 0, . . . , 0]T ∈ Rℓ+1 denote the first axis vector. Then min

✇∈range(Uℓ) A✇ − ❣ = 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 = Pk✈ k; ✈ = A✉k 4. for i = 1, 2, . . . , k do 5. hi,k = ✈ T✈ i; ✈ = ✈ − hi,k✈ i 6. end 7. hk+1,k = ✈; ✈ k+1 = ✈/hk+1,k

• 8. end
• 9. define Uℓ = [✉1, ✉2, . . . , ✉ℓ] ∈ Rn2×ℓ
• 10. define Hℓ+1,ℓ = [hi,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 Pk, a well-established choice in the context of iterated Tikhonov methods is the geometric sequence αk = α0qk−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

• btained at each step k by solving with a few steps of Newton’s method

the next non-linear equation r k − CZ k

circr k = qkr k ,

0 < qk < 1, where · denotes the Euclidean norm. Here the parameter qk depends

• n 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 (C − A)③ ≤ ρcirc A③ , ∀ ③ ∈ Rn . 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

structr k = qkr 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

New sequence

We introduce a new technique for computing αk. Assume that r k = ckαp

k,

where ck is a scalar and p > 0 is a fixed exponent. Letting r k = ε yields αk =

p

ε ck . Assuming that ck ≈ ck−1 = r k−1 αp

k−1

, we obtain the sequence α1 = α0, αk =

p

• ε

r k−1 αk−1, k = 2, 3, . . . , where α0 > 0.

Pietro Dell’Acqua 17 / 33

New sequence

The update from αk−1 to αk is based on the ratio between the noise level ε and the norm of the residual relative to the previous iteration. The heuristic motivation behind this strategy is that αk is computed in

• rder to make the norm of the residual r k move towards the desired

value ε, associated with the noise level of the deblurring problem. The key parameter for the generation of the sequence {αk} is p. In the computed examples, we set α0 = 1 and p = 2.

Pietro Dell’Acqua 18 / 33

Test 1: deblurring problem

5 10 15 20 25 5 10 15 20 25

(a) (b) (c)

Figure: (a) True image; (b) PSF; (c) blurred and noisy image.

Pietro Dell’Acqua 19 / 33

Test 1: αk sequence

5 10 15 20 25 30 0.5 1 1.5 iterations α αk Geo αk DH αk New

Figure: Plot of αk related to F-GMRES in case of geometric sequence, Donatelli-Hanke sequence and new sequence.

Pietro Dell’Acqua 20 / 33

Test 1: restorations

(a) (b) (c)

Figure: Restoration associated with discrepancy principle relative to (a) CGLS (RRE 0.0925, IT 27); (b) F-GMRES New (RRE 0.0910, IT 8); (c) Z k

struct-Landweber New (RRE 0.0921, IT 13).

Pietro Dell’Acqua 21 / 33

Test 1: table

Method RRE best res. IT RRE discrepancy IT CGLS 0.0899 36 0.0925 27 GMRES 0.1473 5 — — F-GMRES Geo 0.0901 13 0.0911 12 F-GMRES DH n/a n/a 0.0907 8 F-GMRES New 0.0901 9 0.0910 8 Z k

struct-Landweber Geo

0.0892 23 0.0912 21 Z k

struct-Landweber DH

n/a n/a 0.0950 13 Z k

struct-Landweber New

0.0891 17 0.0921 13

Table: Relative restoration error (RRE) and iteration number (IT) relative to the best restoration and the one associated with discrepancy principle for different methods.

Pietro Dell’Acqua 22 / 33

Test 1: RRE plots

5 10 15 20 25 30 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 iterations RRE CGLS GMRES F−GMRES Geo F−GMRES DH F−GMRES New Zk

struct−Landweber Geo

Zk

struct−Landweber DH

Zk

struct−Landweber New

Pietro Dell’Acqua 23 / 33

Test 2: deblurring problem

5 10 15 20 25 5 10 15 20 25

(a) (b) (c)

Figure: (a) True image; (b) PSF; (c) blurred and noisy image.

Pietro Dell’Acqua 24 / 33

Test 2: αk sequence

5 10 15 20 25 30 0.2 0.4 0.6 0.8 1 iterations α αk Geo αk DH αk New

Figure: Plot of αk related to F-GMRES in case of geometric sequence, Donatelli-Hanke sequence and new sequence.

Pietro Dell’Acqua 25 / 33

Test 2: restorations

(a) (b) (c)

Figure: Restoration associated with discrepancy principle relative to (a) CGLS (RRE 0.1251, IT 27); (b) F-GMRES New (RRE 0.1175, IT 9); (c) Z k

struct-Landweber DH (RRE 0.1145, IT 10).

Pietro Dell’Acqua 26 / 33

Test 2: table

Method RRE best res. IT RRE discrepancy IT CGLS 0.1129 17 0.1251 27 GMRES 0.1287 50 — — F-GMRES Geo 0.1115 9 0.1167 12 F-GMRES DH n/a n/a 0.1142 6 F-GMRES New 0.1114 7 0.1175 9 Z k

struct-Landweber Geo

0.1093 18 — — Z k

struct-Landweber DH

n/a n/a 11.45 10 Z k

struct-Landweber New

0.1091 11 — —

Table: Relative restoration error (RRE) and iteration number (IT) relative to the best restoration and the one associated with discrepancy principle for different methods.

Pietro Dell’Acqua 27 / 33

Test 2: RRE plots

5 10 15 20 25 30 0.105 0.11 0.115 0.12 0.125 0.13 0.135 0.14 0.145 0.15 iterations RRE CGLS GMRES F−GMRES Geo F−GMRES DH F−GMRES New Zk

struct−Landweber Geo

Zk

struct−Landweber DH

Zk

struct−Landweber New

Pietro Dell’Acqua 28 / 33

Analysis of the new sequence

5 10 15 20 25 30 0.2 0.4 0.6 0.8 1 iterations α p=1 p=2 p=3 p=4 p=5 5 10 15 20 25 30 0.2 0.4 0.6 0.8 1 iterations α p=1 p=2 p=3 p=4 p=5

Figure: Plot of αk related to F-GMRES New considering different values of p for Test 1 (on the left) and Test 2 (on the right).

Pietro Dell’Acqua 29 / 33

Analysis of the new sequence

1 2 3 4 5 0.09 0.092 0.094 0.096 0.098 0.1 0.102 p RRE Best restoration Discrepancy principle 1 2 3 4 5 0.1 0.11 0.12 0.13 0.14 0.15 0.16 p RRE Best restoration Discrepancy principle

Figure: Plot of RRE relative to best restoration and discrepancy one computed by F-GMRES New depending on p for Test 1 (on the left) and Test 2 (on the right).

Pietro Dell’Acqua 30 / 33

Analysis of the new sequence

1 2 3 4 5 2 4 6 8 10 12 14 p iterations Best restoration Discrepancy principle 1 2 3 4 5 4 6 8 10 12 14 16 18 p iterations Best restoration Discrepancy principle

Figure: Plot of iteration number relative to best restoration and discrepancy one computed by F-GMRES New depending on p for Test 1 (on the left) and Test 2 (on the right).

Pietro Dell’Acqua 31 / 33

Conclusions

The problem We have considered the image deblurring problem with a non-symmetric point spread function and Anti-Reflective boundary conditions. Krylov methods Usually the solution of this kind of problem by Krylov subspace methods may require a substantial number of iterations. This can make the restoration of large images expensive. Results We have described a non-stationary preconditioner designed to reduce the number of iterations and discuss how it can be applied in conjunction with the F-GMRES iterative method. Numerical examples have shown that this approach is competitive with the state of the art.

Pietro Dell’Acqua 32 / 33

Pietro Dell’Acqua 33 / 33