Scaled gradient projection methods in image deblurring and denoising - PowerPoint PPT Presentation

Scaled gradient projection methods in image deblurring and denoising Mario Bertero 1 Patrizia Boccacci 1 Silvia Bonettini 2 Riccardo Zanella 3 Luca Zanni 3 1 Dipartmento di Matematica, Università di Genova 2 Dipartimento di Matematica, Università di Ferrara 3 Dipartimento di Matematica, Università di Modena e Reggio Emilia Conference on Applied Inverse Problems, Vienna July 20–24 2009 Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 1 / 26

Outline Examples of Imaging problems 1 Optimization problem 2 Gradient methods and step-length selections 3 Scaled Gradient Projection (SGP) Method 4 Test results 5 Conclusions and Future Works 6 Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 2 / 26

Image Deblurring example Image acquisition model: y = H x + b + n , where: y ∈ R n observed image, H ∈ R n × n blurring operator, b ∈ R n background radiation, n ∈ R n unknown noise. Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 3 / 26

Image Deblurring example Image acquisition model: y = H x + b + n , where: y ∈ R n observed image, H ∈ R n × n blurring operator, b ∈ R n background radiation, n ∈ R n unknown noise. Goal: Find an “approximation” of the true image x ∈ R n Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 3 / 26

Image Deblurring example Image acquisition model: y = H x + b + n , where: y ∈ R n observed image, H ∈ R n × n blurring operator, b ∈ R n background radiation, n ∈ R n unknown noise. Goal: Find an “approximation” of the true image x ∈ R n Maximum Likelihood Approach (and early stopping) min L y ( x ) sub. to x ∈ Ω Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 3 / 26

Image Denoising example Image acquisition model: y = x + n , where: y ∈ R n observed image, n ∈ R n unknown noise. Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 4 / 26

Image Denoising example Image acquisition model: y = x + n , where: y ∈ R n observed image, n ∈ R n unknown noise. Goal: Remove noise from y ∈ R n , while preserving some features Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 4 / 26

Image Denoising example Image acquisition model: y = x + n , where: y ∈ R n observed image, n ∈ R n unknown noise. Goal: Remove noise from y ∈ R n , while preserving some features Regularized Approach min J ( 0 ) y ( x ) + µ J R ( x ) sub. to x ∈ Ω where J R ( x ) is (for example): || x || 2 Tikhonov, 2 || x || 1 sparsity inducing, � Ω |∇ x | Total Variation. Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 4 / 26

Problem setting Both examples lead to: Constrained optimization problem min f ( x ) sub. to x ∈ Ω Ω is a convex and closed set f ( x ) is countinuously differentiable in Ω Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 5 / 26

Why gradient type methods? Gradient methods are first order optimization methods. Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 6 / 26

Why gradient type methods? Gradient methods are first order optimization methods. pros Simplicity of implementation first order iterative method Low memory requirements suitable to face high dimensional problems Ability to provide medium-accurate solutions Semiconvergence from numerical practice Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 6 / 26

Why gradient type methods? Gradient methods are first order optimization methods. pros Simplicity of implementation first order iterative method Low memory requirements suitable to face high dimensional problems Ability to provide medium-accurate solutions Semiconvergence from numerical practice cons Low convergence rate hundreds or thousands of iterations Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 6 / 26

The Barzilai-Borwein (BB) step-length selection rules Consider the gradient method: x ( k + 1 ) = x ( k ) − α k g ( k ) k = 0 , 1 , . . . , with g ( x ) = ∇ f ( x ) . Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 7 / 26

The Barzilai-Borwein (BB) step-length selection rules Consider the gradient method: x ( k + 1 ) = x ( k ) − α k g ( k ) k = 0 , 1 , . . . , with g ( x ) = ∇ f ( x ) . Problem: How the step-length α k > 0 can be chosen to improve the convergence rate? Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 7 / 26

The Barzilai-Borwein (BB) step-length selection rules Consider the gradient method: x ( k + 1 ) = x ( k ) − α k g ( k ) k = 0 , 1 , . . . , with g ( x ) = ∇ f ( x ) . Solution: Regard the matrix B ( α k ) = ( α k I ) − 1 as an approximation of the Hessian ∇ 2 f ( x ( k ) ) Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 7 / 26

The Barzilai-Borwein (BB) step-length selection rules Consider the gradient method: x ( k + 1 ) = x ( k ) − α k g ( k ) k = 0 , 1 , . . . , with g ( x ) = ∇ f ( x ) . Solution: Regard the matrix B ( α k ) = ( α k I ) − 1 as an approximation of the Hessian ∇ 2 f ( x ( k ) ) Determine α k by forcing a quasi-Newton property on B ( α k ) : BB1 = argmin α ∈ R � B ( α ) s ( k − 1 ) − z ( k − 1 ) � α k or BB2 = argmin α ∈ R � s ( k − 1 ) − B ( α ) − 1 z ( k − 1 ) � , α k where s ( k − 1 ) = x ( k ) − x ( k − 1 ) � and z ( k − 1 ) = ( g ( k ) − g ( k − 1 ) ) . � Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 7 / 26

The BB step-length selection rules (cont.) It follows that: BB1 = s ( k − 1 ) T s ( k − 1 ) BB2 = s ( k − 1 ) T z ( k − 1 ) α k or α k s ( k − 1 ) T z ( k − 1 ) z ( k − 1 ) T z ( k − 1 ) where s ( k − 1 ) = x ( k ) − x ( k − 1 ) � and z ( k − 1 ) = ( g ( k ) − g ( k − 1 ) ) . � Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 8 / 26

The BB step-length selection rules (cont.) It follows that: BB1 = s ( k − 1 ) T s ( k − 1 ) BB2 = s ( k − 1 ) T z ( k − 1 ) α k or α k s ( k − 1 ) T z ( k − 1 ) z ( k − 1 ) T z ( k − 1 ) where s ( k − 1 ) = x ( k ) − x ( k − 1 ) � and z ( k − 1 ) = ( g ( k ) − g ( k − 1 ) ) . � Remarkable improvements in comparison with the steepest descent method are observed: [Barzilai-Borwein, IMA J. Num. Anal. 1988] [Raydan, IMA J. Num. Anal. 1993] [Friedlander et al., SIAM J. Num. Anal. 1999] [Raydan, SIAM J. Optim. 1997] [Fletcher, Tech. Rep. 207, 2001] [Dai-Liao, IMA J. Num. Anal. 2002] Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 8 / 26

Effective use of the BB rules Further improvements are obtained by using adaptive alternations of the two BB rules; for example:  α k = α BB2 if α BB2 /α BB1 < τ , k k k   α k = α BB1  otherwise ,  k Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 9 / 26

Effective use of the BB rules Further improvements are obtained by using adaptive alternations of the two BB rules; for example:  α k = α BB2 if α BB2 /α BB1 < τ , k k k   α k = α BB1  otherwise ,  k Many suggestions for the alternation are available: [Dai, Optim., 2003] [Dai-Fletcher, Math. Prog. 2005] [Serafini et al., Opt. Meth. Soft. 2005] [Dai et al., IMA J. Num. Anal. 2006] [Zhuo et al., Comput. Opt. Appl., 2006 ] [Frassoldati et al., J. Ind. Manag. Opt. 2008] Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 9 / 26

The BB step-lengths and Scaled Gradient Methods Consider the scaled gradient method: x ( k + 1 ) = x ( k ) − α k D k g ( k ) k = 0 , 1 , . . . , where D k is the symmetric positive definite scaling matrix. Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 10 / 26

The BB step-lengths and Scaled Gradient Methods Consider the scaled gradient method: x ( k + 1 ) = x ( k ) − α k D k g ( k ) k = 0 , 1 , . . . , where D k is the symmetric positive definite scaling matrix. By forcing the quasi-Newton properties on B ( α k ) = ( α k D k ) − 1 we have BB1 = s ( k − 1 ) T D − 1 k D − 1 k s ( k − 1 ) α k s ( k − 1 ) T D − 1 k z ( k − 1 ) and s ( k − 1 ) T D k z ( k − 1 ) BB2 = α k z ( k − 1 ) T D k D k z ( k − 1 ) , where s ( k − 1 ) = x ( k ) − x ( k − 1 ) � and z ( k − 1 ) = ( g ( k ) − g ( k − 1 ) ) . � Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 10 / 26

Scaled Gradient Projection (SGP) method: basic notations [Bonettini et al., Inv. Prob. 2009] Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 11 / 26

Scaled Gradient Projection (SGP) method: basic notations [Bonettini et al., Inv. Prob. 2009] Scaling matrix: D k ∈ D L = { D s.p.d. ∈ R n × n | � D � ≤ L , � D − 1 � ≤ L } , L > 1 , if D k is diagonal, the requirement leads to: L − 1 ≤ ( D k ) ii ≤ L . Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 11 / 26

Scaled Gradient Projection (SGP) method: basic notations [Bonettini et al., Inv. Prob. 2009] Scaling matrix: D k ∈ D L = { D s.p.d. ∈ R n × n | � D � ≤ L , � D − 1 � ≤ L } , L > 1 , if D k is diagonal, the requirement leads to: L − 1 ≤ ( D k ) ii ≤ L . Projection operator: √ P Ω , D ( x ) ≡ argmin y ∈ Ω � x − y � D , where � x � D = x T D x . Zanella (UniMoRe) Gradient projection methods in imaging AIP 2009 11 / 26