Generalized Row-Action Methods for Tomographic Imaging Sparse Tomo - PowerPoint PPT Presentation

Generalized Row-Action Methods for Tomographic Imaging Sparse Tomo Days, Technical University of Denmark, March 27, 2014 Martin S. Andersen joint work with Per Christian Hansen Scientific Computing Section

Outline Algebraic methods for X-ray computed tomography 1 Relaxed incremental proximal methods 2 Numerical experiments 3 1 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

X-ray Computed Tomography Measurement model � � � I 1 I 1 = I 0 exp − µ ( u 1 , u 2 ) ds l log( I 0 /I 1 ) ≈ a T i x I 0 b = Ax + e Parallel beam measurement geometry Detector X-ray source 2 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Algebraic reconstruction technique Projection on hyperplane H i = { x | a T i x = b i } � x − x k � = x k − a i ( a T i x k − b i ) P H i ( x k ) = argmin � a i � 2 x ∈H i Kaczmarz’s method / ART x k +1 = ρP H ik ( x k ) + (1 − ρ ) x k , i k ∈ { 1 , . . . , m } Relaxation parameter ρ ∈ (0 , 2) 3 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — consistent system 17 equations, ρ = 1 . 0 35 equations, ρ = 1 . 0 4 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — consistent system 17 equations, ρ = 0 . 5 35 equations, ρ = 0 . 5 4 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — inconsistent system 17 equations, ρ = 1 . 0 35 equations, ρ = 1 . 0 5 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — inconsistent system 17 equations, ρ = 0 . 5 35 equations, ρ = 0 . 5 5 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — inconsistent system 17 equations, ρ = 0 . 2 35 equations, ρ = 0 . 2 5 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — inconsistent system, randomization 17 equations, ρ = 1 . 0 35 equations, ρ = 1 . 0 6 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — inconsistent system, randomization 17 equations, ρ = 0 . 5 35 equations, ρ = 0 . 5 6 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — inconsistent system, randomization 17 equations, ρ = 0 . 2 35 equations, ρ = 0 . 2 6 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — underdetermined system � 5 0 z 5 � 2 � 2 � 1 � 1 0 0 x 1 1 y 2 2 7 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — consistent system ρ = 1 . 0 8 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — consistent system ρ = 0 . 5 8 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Example — consistent system ρ = 1 . 5 8 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Tomographic image reconstruction Ill-posed inverse problem — we need regularization! Incorporate a priori knowledge in reconstruction problem • spatial information (smoothness, piecewise constant/affine, ...) • bounds (nonnegativity, box constraints, ...) • sparsity • ... Large-scale optimization problem • gradient computation is expensive • may not be differentiable 9 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Superiorization Kaczmarz’s method is pertubation resilient (Herman et al., 2009) � � x k +1 = P x k + t k v k , P = P H m ◦ · · · ◦ P H 2 ◦ P H 1 Converges if Ax = b is consistent and t k → 0 for k → ∞ Perturbed iteration yields a “superior” solution in some sense 10 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Incremental methods (I) m � minimize f i ( x ) i =1 subject to x ∈ C Incremental (sub)gradient iteration: � � x k − t k � � x k +1 = P C ∇ f i k ( x k ) , ∇ f i k ( x k ) ∈ ∂f i k ( x k ) • sublinear rate of convergence — initial convergence often very fast • diminishing stepsize or “oscillation” that depends on stepsize • goes back to Kibardin (1980), Litvakov (1966) • Bertsekas (1996): incremental least-squares and extended Kalman filter 11 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Incremental methods (II) Incremental proximal iteration: � � x − x k � 2 � f i k ( x ) + 1 x k +1 = argmin 2 t k x ∈C equivalently x k +1 = x k − t k g k +1 , g k +1 ∈ ∂f i k ( x k +1 ) + N C ( x k +1 ) Linearized proximal iteration: let g k ∈ ∂f i k ( x k ) � � � x − x k � 2 �� k x + 1 f i k ( x k ) + g T x k +1 = P C argmin 2 t k x ∈ R n � � = P C x k − t k g k 12 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Incremental proximal gradient methods (I) � m � � minimize g i ( x ) + h i ( x ) i =1 x ∈ C subject to Algorithm 1 (Bertsekas, 2011) � � u − x k � 2 � g i k ( u ) + 1 z k = argmin 2 t k u ∈ R n � � z k − t k � x k +1 = P C ∇ h i k ( z k ) Interpretation � � x k − t k � ∇ g i k ( z k ) − t k � x k +1 = P C ∇ h i k ( z k ) 13 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Incremental proximal gradient methods (II) Algorithm 2 (Bertsekas, 2011) z k = x k − t k � ∇ h i k ( x k ) � � u − z k � 2 � g i k ( u ) + 1 x k +1 = argmin 2 t k u ∈C Interpretation � � u − x k � 2 � ∇ h i k ( x k ) T u + 1 g i k ( u ) + � x k +1 = argmin 2 t k u ∈C � � x k − t k � ∇ g i k ( x k +1 ) − t k � = P C ∇ h i k ( x k ) 14 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Relaxed incremental proximal gradient methods Algorithm 1 � � u − x k � 2 � g i k ( u ) + 1 w k = argmin 2 t k u ∈ R n z k = w k − t k � ∇ h i k ( w k ) � � x k +1 = P C ρ z k + (1 − ρ ) x k Algorithm 2 w k = w k − t k � ∇ h i k ( w k ) � � u − w k � 2 � g i k ( u ) + 1 z k = argmin 2 t k u ∈ R n � � x k +1 = P C ρ z k + (1 − ρ ) x k 15 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Convergence results Cyclic order or random order? • cyclic order yields worst-case performance bounds • random order yields expected performance bounds Constant stepsize or diminishing stepize? • constant stepsize yields convergence within an error bound • diminishing stepsize yields exact convergence Randomized cyclic order works well in practice 16 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

R-IPG vs. ART (1 / 2) � m i x − b i ) 2 i =1 ( a T minimize subject to x ∈ C i x − b i ) 2 and h i ( x ) = 0 : Let g i ( x ) = (1 / 2)( a T � � x k − ρa i k ( a T i k x k − b i k ) x k +1 = P C t − 1 + � a i k � 2 k Interpretation: ART with damped step size 17 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

R-IPG vs. ART — numerical example 2D tomography problem • 256 × 256 Shepp–Logan phantom • n = 256 2 = 65536 variables • p = 120 uniformly spaced angles • r = 362 parallel rays • m = pr = 43440 measurements • Gaussian noise: e i ∼ N (0 , 0 . 02 · � b � ∞ ) Algorithm: R-IPG with constant step size 18 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

R-IPG vs. ART — numerical example 1 1 Norm of relative error t k = 0 . 001 t k = 0 . 01 0 . 8 0 . 8 0 . 6 0 . 6 0 . 4 0 . 4 0 . 2 0 . 2 0 2 4 6 8 10 0 2 4 6 8 10 1 1 Norm of relative error t k = 0 . 1 t k = 1 . 0 0 . 8 0 . 8 0 . 6 0 . 6 0 . 4 0 . 4 0 . 2 0 . 2 0 2 4 6 8 10 0 2 4 6 8 10 Cycles Cycles ρ = 0 . 1 ρ = 0 . 5 ρ = 1 . 0 ρ = 1 . 5 ρ = 1 . 9 19 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Regularized data fitting f ( x ) ≡ (1 / 2) � Ax − b � 2 + λ h ( x ) minimize subject to x ∈ C Example : express f ( x ) as m m � � λ (1 / 2)( a T i x − b i ) 2 f ( x ) = + mh ( x ) � �� � � �� � i =1 i =1 g i ( x ) h i ( x ) 20 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014

Total-variation regularization 2 + λ � n (1 / 2) � Ax − b � 2 minimize i =1 � D i x � 2 x ∈ C subject to 2 and h i ( x ) = ( λ/p ) � n Let g i ( x ) = � A i x − b i � 2 i =1 � D i x � 2 � � u − x k � 2 � g i k ( u ) + 1 z k = argmin 2 t k u ∈ R n i k + t − 1 k I ) − 1 ( A i k x k − b i k ) = x k − A T i k ( A i k A T � � ρ ( z k − t k � x k +1 = P C ∇ h i k ( z k )) + (1 − ρ ) x k 21 / 28 DTU Compute, Technical University of Denmark Generalized Row-Action Methods March 27, 2014