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

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.


INTRODUCTION
Block descent methods are useful when one addresses a general optimization problem in which the constraint set, as in several relevant applications, has a separable structure, i.e.Ω = Ω 1 × ... × Ω m , with i=1 n i = n, so that any x ∈ Ω can be block partitioned as x = (x 1 , ..., x m ).Such methods are based on the idea of performing successive minimizations over each block, as in the classical nonlinear Gauss-Seidel method [1]: However, the convergence of this approach is not ensured without quite restrictive convexity assumptions (see [1,2]) and, in addition, computing an exact minimum of f , even if restricted to a single block, can be impractical.
On the other side, inspired by the idea of (2), effective methods able to handle general nonconvex problems and with global convergence properties can be designed [3,4].In particular, in [3] the author proposed a cyclic block gradient projection method, based on inexact solutions of subproblems (2) obtained by exploiting the sufficient decrease property of the Armijo linesearch [1].In a recent paper [5], this approach has been further developed, allowing generalized projections with variable parameters choice.In this work we consider three instances of projections belonging to this class and we apply the resulting schemes to a large-scale problem in astronomy, namely the blind deconvolution of stellar images acquired by a ground-based telescope.Future work will include the convergence analysis of the proposed scheme for Kurdyka-Łojasiewicz functions [6][7][8].

CYCLIC BLOCK GENERALIZED GRADIENT PROJECTION METHOD
We begin this section by defining the set D(Ω) of the functions d : Ω × Ω → R ≥0 which are convex, continuously differentiable and such that for all x, y ∈ Ω.Given an array of parameters σ ∈ S ⊆ R q and a function d σ ∈ D(Ω), we denote by H(f, Ω, S) the set of the metric functions defined as and, for any h σ ∈ H(f, Ω, S), we define the associated generalized gradient projection operator p( • ; h σ ) : Ω → Ω as Examples of metric functions which can be rewritten in the form (4) are: a) the standard Euclidean projection p(x; h σ ) = P Ω (x − σ∇f (x)), obtained by choosing b) the scaled Euclidean projection, considered for example in [9,10], corresponding to the choice In this case the array of parameters σ is given by the pair (α, D), where α ∈ R >0 and D ∈ R n×n is a symmetric positive definite matrix; c) the Bregman distance associated to a strictly convex function b : Ω → R, which is defined as (8) In order to formally introduce our method for problem (1), we choose the metric function h σ ∈ H(f, Ω, S), where S = S 1 × ... × S m , S i ⊂ R qi , such that the parameter σ can be partitioned as σ = (σ 1 , ..., σ m ).Moreover, we define h σ so that it is separable over the m blocks with respect to its first variable, i.e.
where the functions d i σi belong to D(Ω i ) (i = 1, ..., m) and are chosen e.g. as in ( 6)- (8).It is easy to see that the metric function h σ defined in (9) belongs to H(f, Ω, S) and the associated generalized gradient projection can be also partitioned by blocks as The resulting cyclic block generalized gradient projection method is outlined in Algorithm 1.
The convergence result for this scheme has been proved in [5].For sake of completeness, we report the statement of the main theorem.
Theorem 1 Let {x (k) } k∈N be the sequence generated by Algorithm 1 and assume that x is a limit point of {x (k) } k∈N .Then x is a stationary point for problem (1).
Choose the inner iterations number Choose the parameter σ Compute the linesearch parameter λ , where m k,l is the smallest non-negative integer such that i ) i

POISSON BLIND DECONVOLUTION
In the Poisson image blind deconvolution problem, the basic assumption is that the available data g ∈ R n is a realization of a Poisson random variable whose mean is ω * ⊗ x * + be, where ω * ∈ R n is an unknown point spread function (PSF), ⊗ denotes the convolution operator (periodic boundary conditions are assumed here), b is a positive parameter representing the background radiation, e ∈ R n is the vector of all ones and x * ∈ R n is the image we would like to recover.In the following, we will assume that the PSF ω * is normalized to one.Following a maximum likelihood approach, one can consider the optimization problem where KL(x, ω) is the Kullback-Leibler divergence and introduce regularization by solving it approximately through the early stopping of an iterative procedure.
As concerns the feasible sets, we consider non-negativity and flux conservation for the image x, while for the PSF ω we impose non-negativity, normalization to 1 and an upper bound s which can be estimated in the case of adaptive optics devices from the knowledge of the so-called Strehl ratio (SR), i.e. the ratio of peak diffraction intensity of an aberrated versus perfect waveform (see e.g.[11]).These constraints in a blind deconvolution framework have been used e.g. in [12][13][14] and allow good reconstructions even in presence of a large scale nonconvex problem as (10).The resulting sets are then given by We consider a realistic simulation in the astronomical field by a) generating a 512 × 512 image x * of a cluster of 100 stars with different magnitudes (brightest value ≈ 3.In Figure 1 we report both the PSF used in this experiment and the simulated measured image.The reconstruction algorithms are instances of Algorithm 1 where the projection is defined by means of the three choices detailed in section 2, whose main features are detailed below.

Scaled gradient projection (SGP).
The scaled Euclidean projection (7) has been proposed in [10] within the so-called scaled gradient projection method and allowed a notable acceleration of the same method employing the standard Euclidean projection (6), as remarked in several recent papers [13,[15][16][17][18][19][20][21].Here we adopt, at each iteration k, the following diagonal scaling matrix where µ is a prefixed threshold.The steplength parameter α k is then computed by the adaptive alternation of the scaled Barzilai-Borwein (BB) rules as proposed in [10].If the scaling matrix is set equal to the identity we recover the usual gradient projection (GP) method.

Gradient method with Bregman projection (GBP).
A further instance of projections belonging to the family described in section 2 corresponds to the choice of the metric (4) with related distance-like function d σ defined in (8), where b(x) = n i=1 (x i + γ) log(x i + γ), with γ > 0, is the "regularized entropy" [22].The resulting projection operator ( 5) is given by The steplength parameter σ is adaptively computed at each iteration in the following way.First, we observe that, by the Taylor expansion of the exponential function, we have The term q i (σ) can be explicitly expressed also as Then, the GBP method can be considered also as an approximated scaled gradient method employing the following scaling matrix Thus, it is reasonable to determine the steplength parameter according to the quasi-Newton approach where s (k) = x (k) − x (k−1) and w ).The previous one-dimensional minimum problems can be easily solved (for example by means of the fminbnd Matlab function), giving two possible values for the steplength σ k .In our experiments, we adopt the same adaptive alternation strategy applied in [10] for the two Barzilai-Borwein rules.Following the suggestion in [13], we used L (k) 1 = L 1 = 50 (k ∈ N) inner iterations for the image step, L (k) 2 = L 2 = 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) .The outer iterations have been arbitrarily stopped at 3000.In Figure 2 we show the reconstruction of the PSF provided by the three approaches together with the horizontal and vertical central cuts of the pictures compared with those of the target PSF.Moreover, in Figure 3 we plotted the reconstruction errors and the decrease of the objective function versus the number of iterations, where for the PSF we used the standard root mean square error (RMSE) ω (k) − ω * / ω * while for the image we computed the RMSE for each star |x  in this test problem, we can conclude that the choice of the projection strongly affects the behaviour of the minimization algorithm, and the consequences are emphasized in the case of nonconvex objective function with multiple stationary points.The SGP and GBP choices seem to be attracted by the same limit point, even if going through different paths.With these approaches the reconstructions are very satisfactory, since both the image and the PSF are restored with an error below 1%.On the contrary, the standard projection in Euclidean norm leads to a significantly different pair (x, ω), with a higher precision in recovering the correct magnitude of the stars coupled with a worse reconstruction of the PSF (RMSE > 20%), as clearly attested also by the plots shown in the second row of Figure 2.
− x * i |/|x * i | (i = 1, ..., 100) and then calculated the mean of the 100 resulting values.From the results obtained