The aim of this paper is to consider a modification of a block coordinate gradient projection method with Armijo linesearch along the descent direction in which the projection on the feasible set is performed according to a variable non Euclidean metric. The stationarity of the limit points of the resulting scheme has recently been proved under some general assumptions on the generalized gradient projections employed. Here we tested some examples of methods belonging to this class on a blind deconvolution problem from data affected by Poisson noise, and we illustrate the impact of the projection operator choice on the practical performances of the corresponding algorithm

Bonettini, Silvia

Prato, Marco

Rebegoldi, Simone

REBEGOLDI, SIMONE

BONETTINI, Silvia

PRATO, Marco

Archivio istituzionale della ricerca - Università di Modena e Reggio Emilia

18 September 2017intestazione repositorydell’ateneoApplication of cyclic block generalized gradient projection methods to Poisson blind deconvolution / Rebegoldi, Simone;Bonettini, Silvia; Prato, Marco. - ELETTRONICO. - (2015), pp. 225-229. ((Intervento presentato al convegno 2015 23rdEuropean Signal Processing Conference (EUSIPCO) tenutosi a Nice nel Aug. 31-Sept. 4 2015.OriginalApplication of cyclic block generalized gradient projection methods to Poisson blind deconvolutionPublisher:PublishedDOI:10.1109/EUSIPCO.2015.7362378Terms of use:openAccessPublisher copyright(Article begins on next page)Testo definito dall’ateneo relativo alle clausole di concessione d’usoAvailability:This version is available at: 11380/1073401 since: 2017-07-27T09:48:23ZIEEE Institute of Electrical and Electronics Engineers Inc.This is a pre print version of the following article:APPLICATION OF CYCLIC BLOCK GENERALIZED GRADIENTPROJECTION METHODS TO POISSON BLIND DECONVOLUTIONS. Rebegoldi∗, S. Bonettini†, M. Prato∗∗ Universita` di Modena e Reggio Emilia, Dipartimento di Scienze Fisiche, Informatiche e Matematiche,Via Campi 213/b, 41125 Modena, Italy† Universita` di Ferrara, Dipartimento di Matematica e Informatica,Via Saragat 1, 44122 Ferrara, ItalyABSTRACTThe aim of this paper is to consider a modification of a blockcoordinate gradient projection method with Armijo linesearchalong the descent direction in which the projection on the fea-sible set is performed according to a variable non Euclideanmetric. The stationarity of the limit points of the resultingscheme has recently been proved under some general assump-tions on the generalized gradient projections employed. Herewe tested some examples of methods belonging to this classon a blind deconvolution problem from data affected by Pois-son noise, and we illustrate the impact of the projection opera-tor choice on the practical performances of the correspondingalgorithm.Index Terms— Constrained optimization, gradient pro-jection methods, alternating algorithms, nonconvex optimiza-tion.1. INTRODUCTIONBlock descent methods are useful when one addresses a gen-eral optimization problemminx∈Ωf(x), (1)in which the constraint set, as in several relevant applications,has a separable structure, i.e. Ω = Ω1 × ... × Ωm, withΩi ⊆ Rni,∑mi=1 ni = n, so that any x ∈ Ω can be blockpartitioned as x = (x1, ...,xm).Such methods are based on the idea of performing successiveminimizations over each block, as in the classical nonlinearGauss-Seidel method [1]:x(k+1)i ∈ arg minx∈Ωif(x(k+1)1 , ...,x(k+1)i−1 ,x,x(k)i+1, ...,x(k)m ).(2)This work has been partially supported by MIUR under the two projectsFIRB - Futuro in Ricerca 2012, contract RBFR12M3AC and PRIN 2012,contract 2012MTE38N, and by the Italian Spinner2013 PhD Project “High-complexity inverse problems in biomedical applications and social systems”.The Italian GNCS - INdAM is also acknowledged.However, the convergence of this approach is not ensuredwithout quite restrictive convexity assumptions (see [1, 2])and, in addition, computing an exact minimum of f , even ifrestricted to a single block, can be impractical.On the other side, inspired by the idea of (2), effective meth-ods able to handle general nonconvex problems and withglobal convergence properties can be designed [3, 4]. Inparticular, in [3] the author proposed a cyclic block gradientprojection method, based on inexact solutions of subproblems(2) obtained by exploiting the sufficient decrease property ofthe Armijo linesearch [1]. In a recent paper [5], this approachhas been further developed, allowing generalized projectionswith variable parameters choice. In this work we considerthree instances of projections belonging to this class and weapply the resulting schemes to a large-scale problem in as-tronomy, namely the blind deconvolution of stellar imagesacquired by a ground-based telescope. Future work will in-clude the convergence analysis of the proposed scheme forKurdyka-Łojasiewicz functions [6–8].2. CYCLIC BLOCK GENERALIZED GRADIENTPROJECTION METHODWe begin this section by defining the set D(Ω) of the func-tions d : Ω × Ω → R≥0 which are convex, continuouslydifferentiable and such thatd(x,y) = 0⇔ x = y (3)for all x,y ∈ Ω. Given an array of parameters σ ∈ S ⊆ Rqand a function dσ ∈ D(Ω), we denote by H(f,Ω, S) the setof the metric functions defined ashσ(x,y) = ∇f(y)T (x− y) + dσ(x,y) (4)and, for any hσ ∈ H(f,Ω, S), we define the associated gen-eralized gradient projection operator p( · ;hσ) : Ω→ Ω asp(x;hσ) = argminz∈Ωhσ(z,x) ∀x ∈ Ω. (5)Examples of metric functions which can be rewritten in theform (4) are:a) the standard Euclidean projection p(x;hσ) = PΩ(x −σ∇f(x)), obtained by choosingdσ(x,y) =12σ‖x− y‖2, σ > 0; (6)b) the scaled Euclidean projection, considered for examplein [9, 10], corresponding to the choiced(α,D)(x,y) =12α(x− y)TD−1(x− y). (7)In this case the array of parameters σ is given by the pair(α,D), where α ∈ R>0 and D ∈ Rn×n is a symmetricpositive definite matrix;c) the Bregman distance associated to a strictly convex func-tion b : Ω→ R, which is defined asdσ(x,y) =1σ(b(x)− b(y)−∇b(y)T (x−y)), σ > 0.(8)In order to formally introduce our method for problem (1),we choose the metric function hσ ∈ H(f,Ω, S), where S =S1 × ... × Sm, Si ⊂ Rqi, such that the parameter σ can bepartitioned as σ = (σ1, ...,σm). Moreover, we define hσ sothat it is separable over the m blocks with respect to its firstvariable, i.e.hσ(x,y) =m∑i=1hiσi(xi,y) (9)=m∑i=1∇if(y)T (xi − yi) + diσi(xi,yi),where the functions diσibelong to D(Ωi) (i = 1, ...,m) andare chosen e.g. as in (6)–(8). It is easy to see that the metricfunction hσ defined in (9) belongs toH(f,Ω, S) and the asso-ciated generalized gradient projection can be also partitionedby blocks asp(x;hσ) =p1(x;h1σ1)...pm(x;hmσm) ,wherepi(x;hiσi) = arg minzi∈Ωihiσi(zi,x).The resulting cyclic block generalized gradient projectionmethod is outlined in Algorithm 1.The convergence result for this scheme has been provedin [5]. For sake of completeness, we report the statement ofthe main theorem.Theorem 1 Let {x(k)}k∈N be the sequence generated by Al-gorithm 1 and assume that x¯ is a limit point of {x(k)}k∈N.Then x¯ is a stationary point for problem (1).Algorithm 1 Cyclic block generalized gradient projectionmethodDefine a compact set S and a metric function hσ ∈H(f,Ω, S) as in (9). Choose x(0) ∈ Ω, the upper boundsfor the inner iterations numbers L1, ..., Lm and β, δ ∈ (0, 1).FOR k = 0, 1, 2, ...FOR i = 1, ...,mSet x(k,0)i = x(k)iChoose the inner iterations number L(k)i ≤ LiFOR ℓ = 0, ..., L(k)i − 1Choose the parameter σ(k,ℓ)i ∈ SiSet x˜(k,ℓ)=(x(k+1)1 , ...,x(k+1)i−1 ,x(k,ℓ)i , ...,x(k)m )Compute the descent directiond(k,ℓ)i = pi(x˜(k,ℓ); hiσ(k,ℓ)i)− x(k,ℓ)iSet d˜(k,ℓ)=(0, ...,0,d(k,ℓ)i , ...,0)Compute the linesearch parameter λ(k,ℓ)i = δmk,l ,where mk,l is the smallest non-negative integersuch thatf(x˜(k,ℓ) + λ(k,ℓ)i d˜(k,ℓ)) ≤ f(x˜(k,ℓ))+ βλ(k,ℓ)i ∇f(x˜(k,ℓ))T d˜(k,ℓ)Set x(k,ℓ+1)i = x(k,ℓ)i + λ(k,ℓ)i d(k,ℓ)iENDSet x(k+1)i = x(k,L(k)i)iENDEND3. POISSON BLIND DECONVOLUTIONIn the Poisson image blind deconvolution problem, the basicassumption is that the available data g ∈ Rn is a realizationof a Poisson random variable whose mean is ω∗ ⊗ x∗ + be,where ω∗ ∈ Rn is an unknown point spread function (PSF),⊗ denotes the convolution operator (periodic boundary condi-tions are assumed here), b is a positive parameter representingthe background radiation, e ∈ Rn is the vector of all ones andx∗ ∈ Rn is the image we would like to recover. In the fol-lowing, we will assume that the PSF ω∗ is normalized to one.Following a maximum likelihood approach, one can considerthe optimization problemminx∈Ωx,ω∈ΩωKL(x,ω), (10)where KL(x,ω) is the Kullback–Leibler divergenceKL(x,ω) =n∑i=1gi log(gi(ω ⊗ x)i + b)(11)+ (ω ⊗ x)i + b− gi,and introduce regularization by solving it approximatelythrough the early stopping of an iterative procedure.As concerns the feasible sets, we consider non-negativity andflux conservation for the image x, while for the PSFω we im-pose non-negativity, normalization to 1 and an upper bound swhich can be estimated in the case of adaptive optics devicesfrom the knowledge of the so-called Strehl ratio (SR), i.e.the ratio of peak diffraction intensity of an aberrated versusperfect waveform (see e.g. [11]). These constraints in a blinddeconvolution framework have been used e.g. in [12–14] andallow good reconstructions even in presence of a large scalenonconvex problem as (10). The resulting sets are then givenbyΩx = {x ∈ Rn | x ≥ 0,n∑i=1xi =n∑i=1gi − nb},Ωω = {ω ∈ Rn | 0 ≤ ω ≤ s,n∑i=1ωi = 1}.We consider a realistic simulation in the astronomical fieldby a) generating a 512 × 512 image x∗ of a cluster of 100stars with different magnitudes (brightest value ≈ 3.2 · 107,dimmest value ≈ 4.2 · 106); b) convolving it with a PSFω∗ (SR = 0.81) mimicking the response of a single mirrorof the large binocular telescope (LBT); c) adding a realisticbackground radiation (b ≈ 2.6 · 104); and d) corrupting theresulting blurred image with Poisson noise. In Figure 1 wereport both the PSF used in this experiment and the simulatedmeasured image. The reconstruction algorithms are instancesof Algorithm 1 where the projection is defined by means ofthe three choices detailed in section 2, whose main featuresare detailed below.Scaled gradient projection (SGP). The scaled Euclideanprojection (7) has been proposed in [10] within the so-calledscaled gradient projection method and allowed a notableacceleration of the same method employing the standardEuclidean projection (6), as remarked in several recent pa-pers [13, 15–21]. Here we adopt, at each iteration k, thefollowing diagonal scaling matrix[Dk]ii = max{1µ,min{µ,x(k)}}, (12)where µ is a prefixed threshold. The steplength parameterαk is then computed by the adaptive alternation of the scaledBarzilai–Borwein (BB) rules as proposed in [10]. If the scal-ing matrix is set equal to the identity we recover the usualgradient projection (GP) method.Gradient method with Bregman projection (GBP). A fur-ther instance of projections belonging to the family describedin section 2 corresponds to the choice of the metric (4) withrelated distance-like function dσ defined in (8), where b(x) =∑ni=1(xi + γ) log(xi + γ), with γ > 0, is the “regularizedentropy” [22]. The resulting projection operator (5) is givenby[p(x, hσ)]i = max{(xi + γ)e−σ∇if(x) − γ, 0}.The steplength parameter σ is adaptively computed at eachiteration in the following way. First, we observe that, by theTaylor expansion of the exponential function, we have(xi + γ)e−σ∇if(x) = (xi + γ)− qi(σ)(xi + γ)∇if(x),where qi(σ) =∑∞j=0(−1)j σj+1(j+1)!∇if(x)j. The termqi(σ) can be explicitly expressed also as qi(σ) = (1 −e−σ∇if(x))/∇if(x) when ∇if(x) 6= 0, qi(σ) = σ when∇if(x) = 0. Then, the GBP method can be considered alsoas an approximated scaled gradient method employing thefollowing scaling matrix[Dk]ii(σ) = (x(k)i + γ)qi(σ). (13)Thus, it is reasonable to determine the steplength parameteraccording to the quasi-Newton approachminσ∈[σmin,σmax]‖Dk(σ)−1s(k) −w(k)‖2,minσ∈[σmin,σmax]‖s(k) −Dk(σ)w(k)‖2,where s(k) = x(k) − x(k−1) and w(k) = ∇f(x(k)) −∇f(x(k−1)). The previous one-dimensional minimum prob-lems can be easily solved (for example by means of thefminbnd Matlab function), giving two possible valuesfor the steplength σk. In our experiments, we adopt thesame adaptive alternation strategy applied in [10] for the twoBarzilai-Borwein rules.Following the suggestion in [13], we used L(k)1 = L1 = 50(k ∈ N) inner iterations for the image step, L(k)2 = L2 = 1(k ∈ N) iteration for the PSF step, a constant image as x(0)and the autocorrelation of the ideal PSF of LBT as ω(0). Theouter iterations have been arbitrarily stopped at 3000.In Figure 2 we show the reconstruction of the PSF providedby the three approaches together with the horizontal and ver-tical central cuts of the pictures compared with those of thetarget PSF. Moreover, in Figure 3 we plotted the reconstruc-tion errors and the decrease of the objective function versusthe number of iterations, where for the PSF we used the stan-dard root mean square error (RMSE) ‖ω(k) − ω∗‖/‖ω∗‖while for the image we computed the RMSE for each star|x(k)i − x∗i |/|x∗i | (i = 1, ..., 100) and then calculated themean of the 100 resulting values. From the results obtainedFig. 1. Star cluster test problem: original PSF (left) and blurred and noisy image (right). Both images are in log scale.244 250 256 262 26800.0050.010.0150.020.0250.03PSF − CENTRAL ROW  TRUEGPSGPGBP244 250 256 262 26800.0050.010.0150.020.025PSF − CENTRAL COLUMN  TRUEGPSGPGBPFig. 2. PSF for the star cluster test problem. First row: reconstructions with GP (left), SGP (middle) and GBP (right) in logscale. Second row: horizontal and vertical central cuts of original and restored PSFs (red and black plots coincide).in this test problem, we can conclude that the choice of theprojection strongly affects the behaviour of the minimiza-tion algorithm, and the consequences are emphasized in thecase of nonconvex objective function with multiple stationarypoints. The SGP and GBP choices seem to be attracted bythe same limit point, even if going through different paths.With these approaches the reconstructions are very satisfac-tory, since both the image and the PSF are restored with anerror below 1%. On the contrary, the standard projection inEuclidean norm leads to a significantly different pair (x,ω),with a higher precision in recovering the correct magnitudeof the stars coupled with a worse reconstruction of the PSF(RMSE > 20%), as clearly attested also by the plots shownin the second row of Figure 2.REFERENCES[1] D. Bertsekas, Nonlinear programming, Athena Scien-tific, Belmont, 1999.[2] L. Grippo and M. Sciandrone, “On the convergence ofthe block nonlinear Gauss-Seidel method under convexconstraints,” Oper. Res. Lett., vol. 26, no. 3, pp. 127–136, 2000.[3] S. Bonettini, “Inexact block coordinate descent meth-ods with application to the nonnegative matrix factor-ization,” IMA J. Numer. Anal., vol. 31, no. 4, pp. 1431–1452, 2011.[4] A. Cassioli, D. Di Lorenzo, and M. Sciandrone, “On theconvergence of inexact block coordinate descent meth-ods for constrained optimization,” Eur. J. Oper. Res.,100 101 102 10310−310−210−1RMSE ON THE IMAGE  GPSGPGBP100 101 102 10310−310−210−1100RMSE ON THE PSF  GPSGPGBP100 101 102 103104105106107108OBJECTIVE FUNCTION DECREASE  GPSGPGBPFig. 3. Blind deconvolution: RMSE on image (left), RMSE on PSF (middle) and objective function decrease (right) in logarith-mic scale versus the iterations number (best RMSE on image in unblind case ≈ 1.6 · 10−4).vol. 231, no. 2, pp. 274–281, 2013.[5] S. Bonettini, M. Prato, and S. Rebegoldi, “A cyclicblock coordinate descent method with generalized gra-dient projections,” ArXiv e-prints, p. 1502.06737, 2015.[6] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alter-nating linearized minimization for nonconvex and non-smooth problems,” Math. Program., vol. 146, no. 1–2,pp. 459–494, 2014.[7] E. Chouzenoux, J.-C. Pesquet, and A. Repetti, “Ablock coordinate variable metric forward-backward al-gorithm,” <hal-00945918>, Feb. 2014.[8] P. Frankel, G. Garrigos, and J. Peypouquet, “Splittingmethods with variable metric for Kurdyka–Łojasiewiczfunctions and general convergence rates,” J. Optim.Theory Appl., vol. 165, no. 3, pp. 874–900, 2015.[9] E. G. Birgin, J. M. Martinez, and M. Raydan, “Inex-act spectral projected gradient methods on convex sets,”IMA J. Numer. Anal., vol. 23, no. 4, pp. 539–559, 2003.[10] S. Bonettini, R. Zanella, and L. Zanni, “A scaled gra-dient projection method for constrained image deblur-ring,” Inverse Probl., vol. 25, no. 1, pp. 015002, 2009.[11] V. N. Mahajan, “Strehl ratio for primary aberrations interms of their aberration variance,” J. Opt. Soc. Am.,vol. 73, no. 6, pp. 860–861, 1983.[12] G. Desidera` and M. Carbillet, “Strehl-constrained itera-tive blind deconvolution for post-adaptive-optics data,”Astron. Astrophys., vol. 507, no. 3, pp. 1759–1762,2009.[13] M. Prato, A. La Camera, S. Bonettini, and M. Bertero,“A convergent blind deconvolution method for post-adaptive-optics astronomical imaging,” Inverse Probl.,vol. 29, no. 6, pp. 065017, 2013.[14] M. Prato, A. La Camera, S. Bonettini, S. Rebegoldi,M. Bertero, and P. Boccacci, “A blind deconvolutionmethod for ground-based telescopes and Fizeau inter-ferometers,” manuscript submitted for publication.[15] S. Bonettini and M. Prato, “Nonnegative image recon-struction from sparse Fourier data: a new deconvolutionalgorithm,” Inverse Probl., vol. 26, no. 9, pp. 095001,2010.[16] S. Bonettini, A. Cornelio, and M. Prato, “A new semi-blind deconvolution approach for Fourier-based imagerestoration: an application in astronomy,” SIAM J.Imaging Sci., vol. 6, no. 3, pp. 1736–1757, 2013.[17] S. Bonettini and M. Prato, “Accelerated gradient meth-ods for the X-ray imaging of solar flares,” InverseProbl., vol. 30, no. 5, pp. 055004, 2014.[18] I. Loris, M. Bertero, C. De Mol, R. Zanella, andL. Zanni, “Accelerating gradient projection methodsfor ℓ1-constrained signal recovery by steplength selec-tion rules,” Appl. Computat. Harmon. A., vol. 27, no. 2,pp. 247–254, 2009.[19] M. Prato, R. Cavicchioli, L. Zanni, P. Boccacci, andM. Bertero, “Efficient deconvolution methods for as-tronomical imaging: algorithms and IDL-GPU codes,”Astron. Astrophys., vol. 539, pp. A133, 2012.[20] R. Zanella, P. Boccacci, L. Zanni, and M. Bertero, “Ef-ficient gradient projection methods for edge-preservingremoval of Poisson noise,” Inverse Probl., vol. 25, no.4, pp. 045010, 2009.[21] R. Zanella, G. Zanghirati, R. Cavicchioli, L. Zanni,P. Boccacci, M. Bertero, and G. Vicidomini, “Towardsreal-time image deconvolution: application to confocaland sted microscopy,” Sci. Rep., vol. 3, pp. 2523, 2013.[22] A. Nemirovski, “Prox-method with rate of conver-gence o(1/t) for variational inequalities with Lipschitzcontinuous monotone operators and smooth convex–concave saddle point problems,” SIAM J. Optim., vol.15, no. 1, pp. 229–251, 2004.

Application of cyclic block generalized gradient projection methods to Poisson blind deconvolution

Archivio istituzionale della ricerca - Università di Ferrara

Application of cyclic block generalized gradient projection methods to poisson blind deconvolution

Application of cyclic block generalized gradient projection methods to Poisson blind deconvolution

Abstract

Similar works

Full text

Available Versions

Archivio istituzionale della ricerca - Università di Modena e Reggio Emilia

Archivio istituzionale della ricerca - Università di Ferrara