Parameter estimation in regularization models for Poisson data L. - PowerPoint PPT Presentation

Parameter estimation in regularization models for Poisson data L. Zanni Department of Physics, Computer Science and Mathematics, University of Modena and Reggio Emilia, Italy First French-German Mathematical Image Analysis Conference Paris, 13 - 15 January 2014 Joint work with: M. Bertero , University of Genova, Italy V. Ruggiero , University of Ferrara, Italy L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 1 / 49

Outline Poisson data and the MAP optimization problem 1 Regularization parameter estimation 2 A possible discrepancy principle Solving the discrepancy equation 3 Gradient-based methods for differentiable problems Updating the steplength parameter Updating the scaling matrix Secant-based root finding solver Solving a constrained model 4 ADMM as optimization solver Numerical considerations 5 Conclusions 6 L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 2 / 49

Poisson data ➤ Consider imaging processes where image intensity is measured via the counting of incident particles ➤ Fluctuations in the emission-counting process can be described by modeling the data as realizations of Poisson random variables ➤ The probability of receiving n particles is given by p ( n ) = e − λ λ n , n = 0 , 1 , 2 , . . . n ! where λ is the expected value of the counts ➤ A statistical model appropriate for describing data in different imaging applications ( emission tomography, fluorescence microscopy, optical/infrared astronomy, etc. ) L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 3 / 49

Poisson noisy image restoration: problem setting ➤ x ∗ ∈ R n − → the unknown true image; x ∗ i ≥ 0 ➤ A ∈ R n × n − → the imaging matrix representing the blurring phenomenon ( A = I in denoising) ➤ b ∈ R − → the nonnegative background radiation ➤ ( A x ∗ + b ) − → the image that would be recorded in absence of noise ➤ y ∈ R n − → the observed blurred noisy image; y i ≥ 0 ➩ Given A, b, y , determine an estimate of the true image x ∗ L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 4 / 49

Poisson noisy image restoration: the optimization problem Following the maximum a posteriori (MAP) approach, an estimate x ∗ β of the unknown image can be obtained by solving The MAP constrained optimization problem min x ≥ 0 f 0 ( x ) + βf 1 ( x ) , β > 0 f 0 ( x ) data-fidelity function (generalized Kullback-Leibler divergence) n � y i � � f 0 ( x ) = D KL ( y , x ) = y i log + ( A x + b ) i − y i ( A x + b ) i i =1 f 1 ( x ) = f ( L x ) regularization term ( L linear operator) f 1 ( x ) = 1 2 � x � 2 2 , f 1 ( x ) = � x � 1 , f 1 ( x ) = �|∇ x |� 1 β is the regularization parameter L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 5 / 49

Regularization parameter estimation ➤ Look for estimations of a regularization parameter β suitable for balancing the data fidelity with the regularity of the solution [Engl-Hanke-Neubauer,1996], [Bertero-Boccacci,1998] ➤ Focus on ideas based on the discrepancy principle: a suitable β is such that a measure of the discrepancy D y ( β ) between the corrupted data y and the reconstruction x ∗ β equals some known error τ : D y ( β ) = τ ➤ We consider n � � y i � D y ( β ) = D KL ( y , x ∗ + ( A x ∗ β ) = y i log β + b ) i − y i ( A x ∗ β + b ) i i =1 A i,j ≥ 0 , � � ∀ i, j, i A i,j = 1 , j A i,j > 0 , b > 0 [Bardsley-Goldes,2009], [Bertero et al., 2010], [Carlavan, Blanc-F´ eraud, 2011,2012], [Teuber-Steidl-Chan, 2013] L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 6 / 49

Following the discrepancy principle: what is necessary? n � � y i � + ( A x ∗ D y ( β ) = y i log β + b ) i − y i = τ ( A x ∗ β + b ) i i =1 ➤ Find a suitable value for the constant τ , given the Poisson data assumption and the above discrepancy function. ➤ Exploit effective algorithms for finding β such that D y ( β ) = τ . Alternative approaches: Determine β by solving directly the nonlinear equation [Zanella-Boccacci-Z.-Bertero 2009], [Bertero et al., 2010] Determine β by solving the constrained minimization problem � � min f 1 ( x ) sub. to D y ( β ) ≤ τ x ≥ 0 [Carlavan, Blanc-F´ eraud, 2011,2012], [Teuber-Steidl-Chan, 2013] L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 7 / 49

A possible setting for the constant τ Lemma ( [Zanella et al., Inverse Problems 2009, 2013] ) Let Y λ be a Poisson r.v. with expected value λ and consider � � Y λ � � F ( Y λ ) = 2 Y λ ln + λ − Y λ . λ Then the expected value of F ( Y λ ) satisfies � 1 � E { F ( Y λ ) } = 1 + O , λ → + ∞ . λ The asymptotic estimate of the expected value of D y ( β ) is n 2 . Then τ = n 2 eraud, IEEE T. Image Proc. 2011] τ = m 2 , m = # { y i , y i > 0 } In [Carlavan, Blanc-F´ L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 8 / 49

Find ¯ β such that D y (¯ β ) = τ Consider the nonlinear equation n � � y i � + ( A x ∗ D y ( β ) − τ = y i log β + b ) i − y i − τ = 0 ( A x ∗ β + b ) i i =1 where x ∗ β = argmin f β ( x ) = f 0 ( x ) + βf 1 ( x ) x ≥ 0 and f 0 ( x ) is nonnegative, convex and coercive on R n ≥ 0 . Assume f 1 ( x ) differentiable, nonnegative, convex and such that ∇ 2 f 0 ( x ) ∇ 2 f 1 ( x ) � � � � N ∩ N = { 0 } ➩ f β ( x ) is coercive and strictly convex ⇒ D y ( β ) is well-defined if ¯ D y ( β ) is an increasing function of β ⇒ β exists, it is unique L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 9 / 49

Look at an example: edge preserving regularization n � � f 1 ( x ) = ∆ k k =1 ∆ k = ( x i +1 ,j − x i,j ) 2 + ( x i,j +1 − x i,j ) 2 + δ 2 = � L k x � 2 + δ 2 � � 1 n ] T , L = [ L T 1 , . . . , L T ( L k ) 2 × n , E ( x ) = diag ∆ k I 2 2 k =1 ,...,n � ∆ k A k xx T A T � 1 F ( x ) = diag I 2 − k k =1 ,...,n ∇ f 1 ( x ) = L T E ( x ) L x , ∇ 2 f 1 ( x ) = L T E ( x ) − 1 F ( x ) L , ➩ = { x ∈ R n | x = c 1 n , c ∈ R n } ∇ 2 f 1 ( x ) N � � L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 10 / 49

Look at an example: edge preserving regularization Recall that n � y i � � f 0 ( x ) = y i log + ( A x + b ) i − y i ( A x + b ) i i =1 A i,j ≥ 0 , � � ∀ i, j, x ≥ 0 i A i,j = 1 , j A i,j > 0 , b > 0 , y y ∇ 2 f 0 ( x ) = A T ∇ f 0 ( x ) = 1 n − A T ( A x + b ) , ( A x + b ) 2 A = { x ∈ R n | ( A x ) i = 0 , i ∈ I 1 } , ∇ 2 f 0 ( x ) � � N I 1 = { i | y i > 0 } ➩ ∇ 2 f 0 ( x ) ∇ 2 f 1 ( x ) � � � � N ∩ N = { 0 } ➤ f β ( x ) is strictly convex and D y ( β ) is an increasing function of β L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 11 / 49

Edge preserving regularization: existence of ¯ β such that D y (¯ β ) = τ √ ∆ k f 1 ( x ) = � n Let k =1 β = x ∗ , x ∗ = argmin β → 0 x ∗ (i) lim f 0 ( x ) x ∗ ≥ 0 j =1 ( A T y ) j > b 1 � n (ii) If then n n A i y i � � β →∞ x ∗ A i = lim β = ¯ c 1 n , ¯ c : c + b = n , A i,j A i ¯ i ∈I 1 j =1 ➩ j =1 ( A T y ) j > b , 1 � n f 0 ( x ∗ ) < n c 1 n ) > n If 2 , f 0 (¯ 2 , n β such that D y (¯ ¯ then β ) = τ exists and is unique y = 1 A i = 1 y − b , � (iii) If then c = ¯ ¯ ¯ i ∈I 1 y i n c 1 n ) > n 1 i ∈I 1 y i ln y i > 1 f 0 (¯ ⇔ � 2 + ¯ y ln¯ y 2 n L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 12 / 49

Blurred images corrupted by Poisson noise: two test problems Test environment: Matlab 7.14.0 on a processor Intel Core i7 CPU Q720 1.60 GHz, 4GB RAM Test problems: Cameraman ( 256 × 256 ), Spacecraft ( 256 × 256 ) Original Image Observed Image L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 13 / 49

Edge preserving regularization: cameraman test problem n = 256 2 , ¯ c = 1407 . 84 f 0 ( x ∗ ) = 16403 . 5 < 32768 = n 2 < 14879957 . 3 = f 0 (¯ c 1 n ) 0.5 0.4 0.3 0.2 Behaviour of 0.1 0 F( β ) F ( β ) = 2 −0.1 nD y ( β ) − 1 −0.2 −0.3 −0.4 −0.5 10 −10 10 −8 10 −6 10 −4 10 −2 β L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 14 / 49

Edge preserving regularization: spacecraft test problem n = 256 2 , ¯ c = 154 . 22 f 0 ( x ∗ ) = 32399 . 6 < 32768 = n 2 < 8022676 . 7 = f 0 (¯ c 1 n ) 0.45 0.4 0.35 Behaviour of 0.3 0.25 F( β ) F ( β ) = 2 0.2 nD y ( β ) − 1 0.15 0.1 0.05 0 −0.05 10 −12 10 −10 10 −8 10 −6 10 −4 10 −2 β L. Zanni Parameter estimation in regularization models for Poisson data Paris, 13-15 January 2014 15 / 49