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

Primal-dual interior-point optimization for a regularized reconstruction of NMR relaxation time distributions ome Idier 2 and Fran¸ Emilie Chouzenoux 1 , Sa¨ ıd Moussaoui 2 , J´ cois Mariette 3 erˆ 1 Universit´ 2 Ecole Centrale de Nantes 3 IRSTEA e Paris-Est 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

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

1. NMR Relaxation How to identify the molecular structure of a material by observing its dynamics? 1.1 Principle z • Static field B 0 ⇒ nuclear spin alignment M 0 B 1 ( z axis) M( τ ) • Short magnetic pulse B 1 ⇒ flip angle Φ Φ B 0 τ • Relaxation: return to the equilibrium state M 1 1. Longitudinal dynamics ( z axis) xy ⇒ T 1 relaxation: x 1 ( τ 1 ) = M z ( τ 1 ) 2. Transverse dynamics ( xy plane) ⇒ T 2 relaxation: x 2 ( τ 1 ) = M xz ( τ 1 ) 3/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

1 ) One-dimensionnal analysis � Find T1 or T2 relaxation time constants distribution � x i ( τ i ) = k i ( τ i , T i ) S ( T i ) dT i − → y = Ks + e with k 1 ( τ 1 ) = 1 − (1 − cos Φ) e − τ 1 /T 1 in T 1 relaxation and k 2 ( τ 2 ) = e − τ 2 /T 2 for T 2 relaxation Longitudinal relaxation signal T1 distribution 3 80 2.5 60 2 S 1 x 1 40 1.5 1 20 0.5 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 Repetition time τ 1 [s] Relaxation time T 1 [s] − → Numerical inversion of a Laplace transfrom 4/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

2 ) Two-dimensional analysis [English 1991] � Apply two successive magnetic pulses with a predefined time spacing τ 1 � � x ( τ 1 , τ 2 ) = k 1 ( τ 1 , T 1 ) S ( T 1 , T 2 ) k 2 ( τ 2 , T 2 ) dT 1 dT 2 Y = K 1 SK ⊺ 2 + E ⇐ ⇒ y = ( K 1 ⊗ K 2 ) s + e � Find the joint distribution S ( T 1 , T 2 ) of the relaxation time constants T1−T2 data 100 75 2 50 25 1 0 3 0 3 2 4 1 2 2 1 2 1 3 τ 1 [s] 0 0 τ 2 [s] T 1 (s) T 2 (s) 5/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

1.2 Relaxation time estimation ... problem statement • Ill-conditioned matrices K 1 , K 2 . The singular values of K 1 and K 2 decay exponentially ∢ Direct inversion yields unstable results • Large-size problem in the case of T1-T2 analysis ⊛ Typical setup 1. m 1 = 50 repetition time values τ 1 2. m 2 = 5000 echo time instants τ 2 3. N 1 = N 2 = 300 values of T 1 and T 2 ∢ Matrix K := K 1 ⊗ K 2 of size m 1 m 2 × N 1 N 2 contains over 10 10 elements !! 6/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

1.3 Relaxation time estimation ... regularization framework • The relaxation time distribution is a solution of � � F ( s ) = 1 2 � Ks − y � 2 min 2 + βR ( s ) s ∈ R N + 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 7/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

2. Primal-dual interior point optimization 2.1 Problem formulation min s F ( s ) s.t. s � 0 and max g ( λ ) s.t. λ � 0 λ � �� � � �� � Primal problem Dual problem ( L ( s ) := F ( s ) − λ ⊺ s ) where g ( λ ) is the Lagrange dual function: g ( λ ) = inf s � 0 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 . 8/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

2 ) Interior point-algorithm ... four steps [Armand, 2000] Calculate the primal and dual directions ( d s k , d λ ① k ) A Newton step on (C1) and (C2) gives: � � � � � � d s ∇ 2 F ( s k ) − I λ k − ∇ F ( s k ) k = d λ Diag( s k ) µ k − Λ k s k Λ k I k ∢ System of large size ... infeasible in 2D NMR ! ② Find a step-size α k by a backtracking linesearch and Armijo’s condition on: N � F µ ( s , λ ) = F ( s ) + λ ⊺ s − µ k log( λ n s 2 n ) n =1 � � s k +1 = s k + α k d s k , λ k +1 = λ k + α k d λ ③ Update primal and dual variables: k Decrease the perturbation parameter value: µ k +1 = θ λ ⊺ k +1 s k +1 ④ , with θ ∈ ]0 , 1) . N 9/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

3 ) Primal and dual direction calculation k = Diag( s k ) − 1 [ µ k − Λ k s k − Λ k d s A variable substitution gives: d λ k ] , and � � d s ∇ 2 F ( s k ) + Diag( s k ) − 1 Λ k k = −∇ F ( s k ) + Diag( s k ) − 1 µ k where ∇ 2 F ( s k ) = ( K ⊺ 1 K 1 ) ⊗ ( K ⊺ 2 K 2 ) + ∇ 2 R ( s k ) ∢ Still remains a huge system in 2D relaxation ∢ Approximate resolution using a preconditioned conjugate gradient algorithm 1. Perform TSVDs of K 1 and K 2 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 d s k is obtained by an approximate resolution. 10/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

3. Illustration 3.1 Synthetic data • Mixture of two Gaussian distributions, T1−T2 distribution T1−T2 data 3 2.5 • ( m 1 = 50 , m 2 =5000) values of ( τ 1 , τ 2 ), 0.8 2 0.6 T 1 [s] 1.5 0.4 • A flip angle Φ = 90 ◦ in the T1-T2 model, 1 0.2 0 0.5 0 15 5 10 10 5 0 τ 1 [s] 0 0.5 1 1.5 2 2.5 3 τ 2 [s] T 2 [s] • Additive Gaussian noise with SNR= 20 dB. 11/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

⊛ Reconstruction 1. ( N 1 = 300 , N 2 = 300) with uniform spacing 2. Use a Tikhonov regularization criterion R ( s ) = � s � 2 2 3. Set the regularization parameter β = 100 (unsupervised tuning [Chouzenoux, 2010] ) Estimated S(T 1 ,T 2 ) s 1 (T 1 ): Est.(−) vs. Ref. (−−) 3 2 1 2.5 0 2 0 0.5 1 1.5 2 2.5 3 T 1 [s] T 1 [s] 1.5 s 2 (T 2 ): Est.(−) vs. Ref. (−−) 1 4 0.5 2 0 0 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 T 2 [s] T 2 [s] 4. Computation time: (20 it., 1 s ) in 1D and (36 it., 55 s) for 2D reconstruction. 12/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

3.2 Real data: Analysis of an organic matter (apple) • Measurements: ( m 1 = 50 , m 2 = 10000) • Reconstruction for N 1 = N 2 = 300 • The flip angle is set to Φ = 85 o • T1-T2 Computation time: 260 s for 61 iterations. Est.T1 and T2 distributions • T1 (resp. T2) computation time: 0 . 3 s for 20 10 21 iterations (resp. 9 s for 35 iterations). 0 0 0.5 1 1.5 2 2.5 3 T 1 [s] 20 10 0 0 0.5 1 1.5 2 2.5 3 T 2 [s] 13/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

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. 14/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada

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. 15/15 ICASSP 2013 • 31st May 2013 • Vancouver, Canada