Inverse problems and large scale optimization Original image - - PowerPoint PPT Presentation

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 1/26

A BLOCK PARALLEL MAJORIZE-MINIMIZE MEMORY GRADIENT ALGORITHM

Emilie Chouzenoux, LIGM, UPEM

(joint work with Sara Cadoni, Jean-Christophe Pesquet and Caroline Chaux)

S´ eminaire Parisien des Math´ ematiques Appliqu´ ees ` a l’Imagerie Institut Henri Poincar´ e

3 November 2016

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 2/26

Inverse problems and large scale optimization

Original image Degraded image

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 2/26

Inverse problems and large scale optimization

Original image Degraded image x ∈ RN y = D(Hx) ∈ RM

◮ H ∈ RM×N: matrix associated with the degradation

• perator.

◮ D: RM → RM: noise degradation.

How to find a good estimate of x from the observations y and the model H in the context of large scale processing?

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 3/26

Inverse problems and large scale optimization

Variational approach:

An image estimate ˆ x ∈ RN is generated by minimizing (∀x ∈ RN) F(x) =

S

• s=1

fs(Lsx) with fs : RPs → R, Ls ∈ RPs×N, Ps > 0.

In the context of maximum a posteriori estimation :

◮ L1: Degradation operator, i.e. H; ◮ f1: Data fidelity (e.g. least squares); ◮ (fs)2sS: Regularization functions on some linear transforms (Ls)2sS of the sought solution.

→ Often no closed form expression or solution expensive to compute (especially in large scale context). ◮ Need for an efficient iterative minimization strategy !

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 4/26

Outline

∗ MAJORIZE-MINIMIZE MEMORY GRADIENT ALGORITHM

◮ Majorize-Minimize principle ◮ Subspace acceleration ◮ Convergence theorem

∗ BLOCK PARALLEL 3MG ALGORITHM

◮ Block alternating 3MG ◮ Block separable majorant ◮ Practical implementation ◮ Convergence theorem

∗ APPLICATION TO 3D DECONVOLUTION

◮ Variational approach ◮ Parallel implementation ◮ Numerical results

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 5/26

Majorize-Minimize Memory Gradient algorithm

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 6/26

Majorize-Minimize principle

• 1. Find a tractable surrogate for F Majorization step

F(·) Q(·, xk)

xk xk+1

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 6/26

Majorize-Minimize principle

• 1. Find a tractable surrogate for F Majorization step

Quadratic tangent majorant of F at xk (∀x ∈ RN) Q(x, xk) = F(xk) + ∇F(xk)⊤(x − xk) + 1 2(x − xk)⊤A(xk)(x − xk) where, for every x ∈ RN, A(x) ∈ RN×N is a symmetric definite positive matrix such that (∀x ∈ RN) Q(x, xk) F(x). ∗ Several methods available to construct matrix A(x) in the context of inverse problems in image processing.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 7/26

Subspace acceleration

• 2. Minimize in a subspace

Minimization step (∀k ∈ N∗) xk+1 ∈ Argmin

x∈ran Dk

Q(x, xk), with Dk ∈ RN×Mk.

◮ ran Dk = RN ⇒ half-quadratic algorithm. ◮ Mk small ⇒ low-complexity per iteration.

Memory-Gradient subspace: Dk =

• [−∇F(xk), xk − xk−1]

if k 1 −∇F(x0) if k = 0 3MG algorithm

(similar ideas in NLCG, L-BFGS, TWIST, FISTA, ...)

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 8/26

3MG algorithm

Initialize x0 ∈ RN For k = 0, 1, 2, . . .                Compute ∇F(xk) If k = 0

• Dk = −∇F(x0)

Else

• Dk = [−∇F(xk), xk − xk−1]

Sk = D⊤

k A(xk)Dk

uk = S†

kD⊤ k ∇F(xk)

xk+1 = xk + Dkuk

Low computational cost since Sk is of dimension Mk × Mk, with Mk ∈ {1, 2}. Complexity reductions possible by taking into account the structures of F and Dk.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 9/26

Convergence theorem

Let assume that:

• 1. F : RN → R is a coercive, differentiable function.
• 2. There exists (ν, ν) ∈]0, +∞[2 such that (∀k ∈ N)

ν Id A(xk) ν Id, Then, the following hold:

• ∇F(xk) → 0 and F(xk) ց F(

x) where x is a critical point of F.

• If F is convex, any sequential cluster point of (xk)k∈N is a

minimizer of F.

• If F is strongly convex, then (xk)k∈N converges to the

unique (global) minimizer x of F

• If F satisfies the Kurdyka-Łojasiewicz inequality, then the

sequence (xk)k∈N converges to a critical point of F.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 10/26

3MG in practical situations

3MG algorithm outperforms state-of-the arts optimization algorithms in many image processing applications. Problem: Computational issues with very large-size problems. Main reasons:

◮ High computational time for calculating the gradient

direction ∇F(xk) and the matrix Sk = D⊤

k A(xk)Dk; ◮ High storage cost for ∇F(xk), Dk and xk.

↓ Block parallel approach

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 11/26

Block parallel 3MG algorithm

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 12/26

Block parallel strategy

The vector of unknowns x is partitioned into block subsets. At each iteration, some blocks are updated in parallel. Advantages:

◮ Control of the memory thanks to the block alternating strategy; ◮ Reduction of the computational time thanks to the parallel procedure.

x = x(1) x(j) x(J) x(S) = (xp)p∈S =

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 13/26

Block alternating 3MG

• 1. Select a block subset: Choose a non empty Sk ⊂ {1, . . . , J}.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 13/26

Block alternating 3MG

• 1. Select a block subset: Choose a non empty Sk ⊂ {1, . . . , J}.
• 2. Find a tractable surrogate in this subset:

Set A(Sk)(xk) = ([A(xk)]p,p)p∈Sk. The restriction of F to Sk is majorized at xk by (∀v ∈ R|Sk|) Q(Sk)(v, xk) = F(xk) + ∇F (Sk)(xk)⊤(v − x(Sk)

k

) + 1 2(v − x(Sk)

k

)⊤A(Sk)(xk)(v − x(Sk)

k

).

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 13/26

Block alternating 3MG

• 1. Select a block subset: Choose a non empty Sk ⊂ {1, . . . , J}.
• 2. Find a tractable surrogate in this subset:

Set A(Sk)(xk) = ([A(xk)]p,p)p∈Sk. The restriction of F to Sk is majorized at xk by (∀v ∈ R|Sk|) Q(Sk)(v, xk) = F(xk) + ∇F (Sk)(xk)⊤(v − x(Sk)

k

) + 1 2(v − x(Sk)

k

)⊤A(Sk)(xk)(v − x(Sk)

k

).

• 3. Minimize within the memory gradient subspace

x(Sk)

k+1 =

Argmin

v∈ran D(Sk)

k

Q(Sk)(v, xk) where (∀j ∈ Sk) D(j)

k

=

• −∇F (j)(xk)

if j / ∈ k−1

ℓ=0 Sℓ,

• − ∇F (j)(xk)
• x(j)

k

− x(j)

k−1]

• therwise.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 13/26

Block alternating 3MG

• 1. Select a block subset: Choose a non empty Sk ⊂ {1, . . . , J}.
• 2. Find a tractable surrogate in this subset:

Set A(Sk)(xk) = ([A(xk)]p,p)p∈Sk. The restriction of F to Sk is majorized at xk by (∀v ∈ R|Sk|) Q(Sk)(v, xk) = F(xk) + ∇F (Sk)(xk)⊤(v − x(Sk)

k

) + 1 2(v − x(Sk)

k

)⊤A(Sk)(xk)(v − x(Sk)

k

).

• 3. Minimize within the memory gradient subspace

x(Sk)

k+1 =

Argmin

v∈ran D(Sk)

k

Q(Sk)(v, xk) where (∀j ∈ Sk) D(j)

k

=

• −∇F (j)(xk)

if j / ∈ k−1

ℓ=0 Sℓ,

• − ∇F (j)(xk)
• x(j)

k

− x(j)

k−1]

• therwise.

Problem: Matrices A(S) do not have any block diagonal structure = ⇒ Difficult to perform Step 3 in parallel !

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 14/26

Block separable majorant matrix

Let us assume that: (∀x ∈ RN) A(x) =

S

• s=1

L⊤

s Diag {ωs(Lsx)} Ls,

Let S ⊂ {1, . . . , J} non empty. Then, (∀x ∈ RN) A(S)(x)B(S)(x) = BDiag

• B(j)(x)
• j∈S
• ,

where, for every j ∈ S, matrix B(j)(x) ∈ RNj×Nj is given by: B(j)(x) =

S

• s=1
• (L(j)

s )⊤Diag {bs(Lsx)} L(j) s

• ,

with, for every s ∈ {1, . . . , S} and p ∈ {1, . . . , Ps}, [bs(Lsx)]p = [ωs(Lsx)]p[|L(S)

s

|1|S|]p/[|L(j)

s |1Nj ]p.

Proof: Rely on Jensen’s inequality.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 15/26

BP3MG Algorithm

Initialize x0 ∈ RN For k = 0, 1, 2, . . .                     Select Sk ⊂ {1, . . . , J} s.t. |Sk| = C Parfor j ∈ Sk            Compute ∇F (j)(xk) Compute B(j)

k (xk)

Construct D(j)

k

u(j)

k

= −

• (D(j)

k )⊤B(j) k (xk)D(j) k

† (D(j)

k )⊤∇F (j)(xk)

x(j)

k+1 = x(j) k

+ D(j)

k u(j) k

Set , for every j ∈ {1, . . . , J} \ Sk, x(j)

k+1 = x(j) k .

Share

• x(j)

k+1

• j∈Sk

between all cores.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 16/26

Practical implementation

In practice, it is usually not necessary to send the full vector xk+1 to all the cores, at each iteration k. ◮ The update of the j-th block only require the knowledge of the current iterate at indices Nj =

S

• s=1
• n ∈ {1, . . . , N}|(∃p ∈ Ps,j) [Ls]p,n = 0
• ,

where Ps,j = {p ∈ {1, . . . , Ps} |(∃i ∈ Jj) [Ls]p,i = 0}. ∗ The cardinality of Nj is usually very small with respect to N. Example: S = 1 and L1 is a discrete gradient operator with

• ne pixel neighborhood ⇒ |Nj| = 3.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 17/26

Convergence theorem (ongoing work)

Let assume that:

• 1. F : RN → R is a coercive, differentiable function.
• 2. There exists a constant K J such that, for every k ∈ N,

{1, . . . , J} ⊂ k+K−1

ℓ=k

Sℓ.

• 3. There exists (ν, ν) ∈]0, +∞[2 such that (∀k ∈ N)

ν Id B(Sk)(xk) ν Id, Then, the following hold:

• ∇F(xk) → 0 and F(xk) ց F(

x) where x is a critical point of F.

• If F is convex, any sequential cluster point of (xk)k∈N is a

minimizer of F.

• If F is strongly convex, then (xk)k∈N converges to the

unique (global) minimizer x of F

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 18/26

Application to 3D image deconvolution

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 19/26

Problem statement

➜ ➜ Original 3D image Degradations Measured 3D image x ∈ RN H ∈ RN×N, b ∈ RN y = Hx + b ◮ H: 3D convolution operator representing depth-variant 3D Gaussian blur (kernel size 5 × 5 × 11). For each depth z ∈ {1, . . . , NZ}, different variance and rotation parameters. ◮ b: additive Gaussian i.i.d. zero-mean noise.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 20/26

Variational approach

(∀x ∈ RN) F(x) = 1 2Hx − y2 + R(x) OBJECTIVE FUNCTION Hybrid penalization term R = R1 + R2 + R3:

◮ R1(x) = η N n=1 d2 [xmin,xmax](xn) ◮ R2(x) = λ N n=1

• ([V Xx]n)2 + ([V Yx]n)2 + δ2

◮ R3(x) = κ N n=1([V Zx])2

• (η, λ, δ, κ) ∈ (0, +∞)4: regularization parameters;
• [xmin, xmax]: range of pixel intensity values; dC: distance to set C;
• V X,V Y,V Z ∈ RN×N: discrete gradients along X,Y and Z.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 21/26

Parallel implementation

◮ Blocks: NZ slices of the 3D

volume.

◮ Message Passing Interface

command SPMD of MATLAB R

• ◮ Master-Slave implementation:

➜ 1 master core: Main loop of the algorithm. ➜ C slave cores: Perform their tasks simultaneously.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 21/26

Parallel implementation

◮ Blocks: NZ slices of the 3D

volume.

◮ Message Passing Interface

command SPMD of MATLAB R

• ◮ Master-Slave implementation:

➜ 1 master core: Main loop of the algorithm. ➜ C slave cores: Perform their tasks simultaneously.

Iteration k x(1) x(J) Slave 1 Slave 2 Slave 3 Slave 4

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 21/26

Parallel implementation

◮ Blocks: NZ slices of the 3D

volume.

◮ Message Passing Interface

command SPMD of MATLAB R

• ◮ Master-Slave implementation:

➜ 1 master core: Main loop of the algorithm. ➜ C slave cores: Perform their tasks simultaneously.

Iteration k + 1 x(1) x(J) Slave 1 Slave 2 Slave 3 Slave 4

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 22/26

Restoration results: FlyBrain

(a) Original (b) Degraded (c) Restored Images corresponding to slice z = 18 of the 3D volume FlyBrain (256 × 256 × 48). Initial SNR 13.42 dB. Final SNR 16.98 dB.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 23/26

Restoration results: Tube

(a) Original (b) Degraded (c) Restored Images corresponding to slice z = 31 of the 3D volume Tube (284 × 280 × 48). Initial SNR 11.53 dB. Final SNR 14.47 dB.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 24/26

Acceleration

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 1 2 3 4 5 6 7 8 9 10 11

Core number C Acceleration

(a) FlyBrain

Core number C

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

Acceleration

1 2 3 4 5 6 7 8 9 10

(b) Tube

Ratio between the computation time for one core and the computation time for C cores (+) with linear fitting (· · · ).

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 25/26

Conclusion

The Block Parallel Majorize-Minimize Memory Gradient (BP3MG) Algorithm handles smooth optimization problems of very large dimension. Reduced complexity / memory requirement. High efficiency in the context of 3D image restoration. Great potential for parallelization. Future work will involve implementation in other languages.

Introduction 3MG Algorithm Block Parallel 3MG Algorithm Experimental results Conclusion Imaging in Paris - IHP 26/26

Some references

• S. Cadoni, E. Chouzenoux, J.-C. Pesquet, C. Chaux.

A Block Parallel Majorize-Minimize Memory Gradient Algorithm. In Proc. IEEE Conf. Image Process. (ICIP 2016), pp. 3194–3198, Phoenix, AZ, 25-28 Sep 2016. Best paper award finalist.

• E. Chouzenoux, J. Idier, S. Moussaoui

A Majorize-Minimize strategy for subspace optimization applied to image restoration in IEEE Trans. Image Process., vol.20, no.18, pp. 1517-1528, 2011.

• E. Chouzenoux,A. Jezierska, J.-C. Pesquet, H. Talbot

A Majorize-Minimize Subspace Approach for ℓ2-ℓ0 Image Regularization in SIAM J. Imag. Sci., vol.6, no.1, pp. 563-591,2013.

• E. Chouzenoux, L. Lamass´

e, S. Anthoine, C. Chaux, A. Jaouen, I. Vanzetta, F. Debarbieux Approche variationnelle pour la d´ econvolution rapide de donn´ ees 3D en microscopie biphotonique in Actes du 25e colloque GRETSI, Lyon, France, 8-11 septembre 2015.

• A. Florescu, E. Chouzenoux, J.-C. Pesquet, P

. Ciuciu and S. Ciochina. A Majorize-Minimize Memory Gradient Method for Complex-Valued Inverse Problems in Signal Processing, Vol. 103, pp. 285-295, 2014.