[PPT] - Primal-dual interior-point optimization for a regularized PowerPoint Presentation

SLIDE 1

Primal-dual interior-point optimization for a regularized reconstruction of NMR relaxation time distributions

Emilie Chouzenoux1, Sa¨ ıd Moussaoui2, J´ erˆ

me Idier2 and Fran¸

cois Mariette3

1 Universit´

e Paris-Est

2Ecole Centrale de Nantes 3IRSTEA

IGM, UMR CNRS 8049 IRCCyN, CNRS UMR 6597 UR TERE, IRM Food Marne-La Vall´ ee, France Nantes, France Rennes, France

contact: said.moussaoui@irccyn.ec-nantes.fr 31st May 2013 • Vancouver, Canada Special session – Signal Processing for Chemical Sensing

SLIDE 2

Outline

1. NMR relaxation
Principle
Problem statement
2. Reconstruction method
Optimization framework
Main algorithm steps
3. Illustration
4. Conclusion

2/15

SLIDE 3

1. NMR Relaxation

How to identify the molecular structure of a material by observing its dynamics? 1.1 Principle

M0 M(τ) τ xy z B1 B0 M1 Φ

Static field B0 ⇒ nuclear spin alignment

(z axis)

Short magnetic pulse B1 ⇒ flip angle Φ
Relaxation: return to the equilibrium state
1. Longitudinal dynamics (z axis)

⇒ T1 relaxation: x1(τ1) = Mz(τ1)

2. Transverse dynamics (xy plane)

⇒ T2 relaxation: x2(τ1) = Mxz(τ1)

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

3/15

SLIDE 4

1) One-dimensionnal analysis Find T1 or T2 relaxation time constants distribution xi(τi) =

ki(τi, Ti) S(Ti) dTi −

→ y = Ks + e with k1(τ1) = 1−(1−cos Φ)e−τ1/T1 in T1 relaxation and k2(τ2) = e−τ2/T2 for T2 relaxation

0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 Repetition time τ1 [s] x1 Longitudinal relaxation signal 0.5 1 1.5 2 20 40 60 80 T1 distribution Relaxation time T1 [s] S1

− → Numerical inversion of a Laplace transfrom

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

4/15

SLIDE 5

2) Two-dimensional analysis

[English 1991]

Apply two successive magnetic pulses with a predefined time spacing τ1 x(τ1, τ2) = k1(τ1, T1)S(T1, T2)k2(τ2, T2) dT1 dT2 Y = K1SK⊺

2 + E ⇐

⇒ y = (K1 ⊗ K2)s + e Find the joint distribution S(T1, T2) of the relaxation time constants

1 2 3 2 4 1 2 τ1 [s] T1−T2 data τ2 [s]

1 2 3 1 2 3 25 50 75 100 T2 (s) T1 (s)

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

5/15

SLIDE 6

1.2 Relaxation time estimation ... problem statement

Ill-conditioned matrices K1, K2. The singular values of K1 and K2 decay exponentially

∢ Direct inversion yields unstable results

Large-size problem in the case of T1-T2 analysis

⊛ Typical setup

1. m1 = 50 repetition time values τ1
2. m2 = 5000 echo time instants τ2
3. N1 = N2 = 300 values of T1 and T2

∢ Matrix K := K1 ⊗ K2 of size m1m2 × N1N2 contains over 1010 elements !!

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

6/15

SLIDE 7

1.3 Relaxation time estimation ... regularization framework

The relaxation time distribution is a solution of

min s∈RN+

F(s) = 1

2Ks − y2

2 + βR(s)

where R(s) is a convex and differentiable regularization criterion

− → Solve a non-negativity constrained optimization problem − → Avoid the storing of matrix K in the 2D case

Previous works
1. Data compression and Tikhonov regularization [Venkataramanan, 2002]
2. Maximum entropy and truncated Newton algorithm [Chouzenoux, 2010]
Our proposal

Adopt and adapt an inexact primal-dual interior-point method

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

7/15

SLIDE 8

2. Primal-dual interior point optimization

2.1 Problem formulation min s F(s) s.t. s 0

Primal problem

and max λ g(λ) s.t. λ 0

Dual problem

where g(λ) is the Lagrange dual function: g(λ) = inf s0 (L(s) := F(s) − λ⊺s) 1) Optimality conditions (Karush Kuhn Tucker) (C1) ∇F(s) − λ = 0, (C2) Λs = 0, (C3) s 0, (C4) λ 0 But in practice, take: (C2) Λs = µk with µk > 0 a perturbation parameter such that lim

k→∞ µk = 0. ICASSP 2013 • 31st May 2013 • Vancouver, Canada

8/15

SLIDE 9

2) Interior point-algorithm ... four steps [Armand, 2000] ① Calculate the primal and dual directions (ds

k , dλ k)

A Newton step on (C1) and (C2) gives:

∇2F(sk)

−I ΛkI Diag(sk) ds

k

dλ

k

=
λk − ∇F(sk)

µk − Λksk

∢ System of large size ... infeasible in 2D NMR !

② Find a step-size αk by a backtracking linesearch and Armijo’s condition on: Fµ(s, λ) = F(s) + λ⊺s − µk

N

n=1

log(λns2

n)

③ Update primal and dual variables:

sk+1 = sk + αkds

k, λk+1 = λk + αkdλ k

④

Decrease the perturbation parameter value: µk+1 = θλ⊺

k+1sk+1

N , with θ ∈]0, 1).

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

9/15

SLIDE 10

3) Primal and dual direction calculation A variable substitution gives: dλ

k = Diag(sk)−1 [µk − Λksk − Λkds k] , and

∇2F(sk) + Diag(sk)−1Λk
ds

k = −∇F(sk) + Diag(sk)−1µk

where ∇2F(sk) = (K⊺

1K1) ⊗ (K⊺ 2K2) + ∇2R(sk)

∢ Still remains a huge system in 2D relaxation ∢ Approximate resolution using a preconditioned conjugate gradient algorithm

1. Perform TSVDs of K1 and K2 to construct an efficient preconditioner,
2. Calculate with a low complexity Hessian-vector and Preconditioner-vector products,
3. See the paper for the stopping criteria.

∢ The convergence proof is established when ds

k is obtained by an approximate

resolution.

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

10/15

SLIDE 11

3. Illustration

3.1 Synthetic data

Mixture of two Gaussian distributions,
(m1 = 50, m2=5000) values of (τ1, τ2),
A flip angle Φ = 90◦ in the T1-T2 model,
Additive Gaussian noise with SNR=20 dB.

T2 [s] T1 [s] T1−T2 distribution 0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 5 10 5 10 15 0.2 0.4 0.6 0.8 τ1 [s] T1−T2 data τ2 [s]

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

11/15

SLIDE 12

⊛ Reconstruction

1. (N1 = 300, N2 = 300) with uniform spacing
2. Use a Tikhonov regularization criterion R(s) = s2

2

3. Set the regularization parameter β = 100 (unsupervised tuning [Chouzenoux, 2010])

T2 [s] T1 [s] Estimated S(T1,T2) 0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 1 2 T1 [s] s1(T1): Est.(−) vs. Ref. (−−) 0.5 1 1.5 2 2.5 3 2 4 T2 [s] s2(T2): Est.(−) vs. Ref. (−−)

4. Computation time: (20 it., 1s) in 1D and (36 it., 55 s) for 2D reconstruction.

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

12/15

SLIDE 13

3.2 Real data: Analysis of an organic matter (apple)

Measurements: (m1 = 50, m2 = 10000)
Reconstruction for N1 = N2 = 300
The flip angle is set to Φ = 85o
T1-T2 Computation time: 260 s for 61

iterations.

T1 (resp. T2) computation time: 0.3 s for

21 iterations (resp. 9 s for 35 iterations).

0.5 1 1.5 2 2.5 3 10 20

T1 [s] Est.T1 and T2 distributions

0.5 1 1.5 2 2.5 3 10 20

T2 [s]

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

13/15

SLIDE 14

Summary

⊛ Main contributions

Address the inverse problem of 2D NMR relaxation times estimation
Propose an efficient optimization algorithm for a differentiable convex regularization
Exploit de forward model structure to reduce the computational complexity.

⊛ Future investigations

1. Estimate the flip angle
2. Gaussian noise assumption is not valid on the nuclear spin module
3. Compare with alternative constrained optimization methods.

ICASSP 2013 • 31st May 2013 • Vancouver, Canada

14/15

SLIDE 15

References

[1] A. E. English, K. P. Whittall, M. L. G. Joy, and R. M. Henkelman, Quantitative two-dimensional time correlation relaxometry, Magnetic Resonance in Medecine, vol. 22, pp. 425–434, 1991. [2] L. Venkataramanan, Y. Q. Song, and M. D. H¨ urlimann, Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions, IEEE Transactions on Signal Processing, vol. 50, no. 5,

pp. 1017-1026, 2002.

[3] E. Chouzenoux, S. Moussaoui, J. Idier, and F. Mariette, Efficient maximum entropy reconstruction of nuclear magnetic resonance T1-T2 spectra, IEEE Transactions on Signal Processing, vol. 58, no. 2, pp. 6040-605, December 2010. [4] P. Armand, J. C. Gilbert, and S. Jan-J´ egou, A feasible BFGS interior point algorithm for solving strongly convex minimization problems, SIAM Journal on Optimization, vol. 11, pp. 199-222, 2000.