New convergence results for the scaled gradient projection method

Introduction

Several inverse problems in applied sciences can be addressed by means of a constrained optimization problem

where $\boldsymbol{y}^{(k)}$ is the scaled Euclidean projection of $\boldsymbol{x}^{(k)}-\alpha_{k}D_{k}\nabla f(\boldsymbol{x}^{(k)})$ onto $\Omega$, i.e.

We also recall the definitions of stationary point and descent direction for problem (1) (see for example ).

A point ${\boldsymbol{x}}\in\Omega$ is a stationary point for problem (1) if

Let $\boldsymbol{x}$ be any point of the set $\Omega$.

Finally, we report the definitions of convex, globally and locally Lipschitz and level bounded function.

locally Lipschitz, if for every compact set $K\subseteq\Omega$ there exists $L_{K}>0$ such that

General results about Armijo based gradient projection methods

In this section, we recall the basic properties of the most popular linesearch procedure, the Armijo linesearch, given in Algorithm 1. These results allow to prove a general convergence result which applies to any method where the objective function over two successive iterates decreases at least as it would decrease by applying the Armijo linesearch procedure along a suitable descent direction.

For the linesearch procedure based on the Armijo rule we recall the following basic theorem, which can be derived from known results .

where $\lambda^{(k)}$ is computed with Algorithm 1, then we have

Proof. Inequality (8) can be rewritten as

Summing the previous inequality for $k=0,...,j$ gives

Let $\boldsymbol{x}^{(k)}\in\Omega$ and $\boldsymbol{d}^{(k)}$ be defined as in (2)–(3). Then we have

Moreover, $\boldsymbol{d}^{(k)}=\boldsymbol{0}$ if and only if $\boldsymbol{x}^{(k)}$ is stationary for problem (1).

We are now ready to give the more general convergence result based on the above mentioned properties of the descent direction and with the Armijo rule establishing the sufficient decrease of the objective function. Its proof is omitted since it can be easily derived by the analogous results in .

It is worth stressing that the previous result applies also to nonconvex problems and the gradient of the objective function is not required to be Lipschitz continuous. Moreover, the only limitations to the algorithms parameters choice are that $\alpha_{k}$ and the eigenvalues of $D_{k}$ have to be bounded above and below away from zero. When $\nabla f$ satisfies some Lipschitz property, the next proposition states that the Armijo steplengths are bounded away from zero. This also means that there exists a finite upper bound for the number of backtracking reductions at any iteration (a similar result can be found in [27, Theorem 3.2]). This result will be useful in the convergence rate analysis of the next section.

Assume that $\nabla f$ satisfies one of the following conditions:

$\nabla f$ is globally Lipschitz on $\Omega$;

$\nabla f$ is locally Lipschitz and $f$ is level bounded on $\Omega$.

Proof. If $\nabla f$ is Lipschitz continuous on $\Omega$ with Lipschitz constant $L$, then from the descent lemma [25, p.667] we have

where $\gamma={\mu\alpha_{\max}}$. The previous inequality ensures that the Armijo condition

Convergence analysis of the scaled gradient projection algorithm

Before to give the main convergence result, we prove the following lemma.

Thus, if the series on the right hand side of (19) converges, the quantities $\theta_{k}$ are bounded for all $k$. We observe that, since $\mu_{j}^{2}=1+\zeta_{j}$, by the known limit $\lim_{\zeta_{j}\rightarrow 0}\log(1+\zeta_{j})/\zeta_{j}=1$, the series $\sum_{j=0}^{\infty}\log(\mu_{j}^{2})$ and $\sum_{j=0}^{\infty}\zeta_{j}$ have the same behaviour. Thus, since by hypothesis the latter one is convergent, the theorem follows. $\square$ The next theorem states that, when $f$ is convex and admits finite minimum, if the scaling matrices $D_{k}$ asymptotically reduce to the identity matrix at a certain rate, then the sequence generated by SGP converges to a solution of (1). The line of the proof is similar to that of [28, Theorem 1], which can be considered as a special case of it. After giving the proof of our result, we discuss the relations of our approach with the related work already present in the literature.

Proof. We recall first the basic norm equality

which holds true for any positive definite matrix $E$. Moreover, it is easy to see that if $D_{k}\in\mathcal{M}_{\mu_{k}}$, then $D_{k}^{-1}\in\mathcal{M}_{\mu_{k}}$. Let $\hat{\boldsymbol{x}}\in X^{*}$. By definition of $\boldsymbol{y}^{(k)}$ we have

which, for ${\boldsymbol{x}}=\hat{\boldsymbol{x}}$ gives

where the inequality follows from the convexity of $f$ and the last equality by definition of $\boldsymbol{x}^{(k+1)}$. By equality (20) with ${\boldsymbol{x}}=\boldsymbol{x}^{(k+1)}$, ${\boldsymbol{y}}=\boldsymbol{x}^{(k)}$, ${\boldsymbol{z}}=\hat{\boldsymbol{x}}$, $E={D_{k}^{-1}}$ we obtain

which, since $\lambda^{(k)}\leq 1$, results in

(since $f(\boldsymbol{x}^{(k)})-f(\hat{\boldsymbol{x}})\geq 0$). From the last inequality and in view of (4), it follows that

Recalling that the scalar product at the right-hand-side is nonpositive, since $\mu_{k}\geq 1$ and $\alpha_{k}\leq\alpha_{\max}$ this results in

By repeatedly applying the previous inequality we obtain

where $\theta^{k}_{j}=\prod_{i=j}^{k}\mu_{j}^{2}$. Since $\mu_{j}^{2}\geq 1$, we have $\theta^{k}_{j}\leq\theta_{0}^{k}$, and by Lemma 3.1 we obtain

and can be described by the following iteration

Clearly, when $g$ is the indicator function of the convex set $\Omega$, problem (1) is equivalent to (24) and the SGP iteration can be expressed in the same form of (25). In , the convergence of the iterates (25) is proved for objective functions with Lipschitz continuous gradients, under the condition

where $\zeta_{k}$ is a summable sequence. Variable metrics were considered also in [31, Chapter 5] in the context of subgradient methods for nonsmooth, convex, unconstrained minimization. In this case, setting $D_{k}=B_{k}B_{k}^{T}$, the scaling matrices are assumed to satisfy $\prod_{k=0}^{\infty}\|B_{k+1}^{-1}B_{k}\|^{2}<\infty$ and

We remark that our condition, $D_{k}\in{\mathcal{M}}_{\mu_{k}}$, is quite different from both (26) and (27) since it does not impose a strict connection between the scaling matrices at two successive iterates. This freedom of choosing the metric at each iteration allows for example to adopt a suitable adaptation of a well performing scaling technique, based on a gradient splitting , which may lead to significant improvements of the convergence behaviour, as we will show in section 4. In the following we give a complexity result about SGP, showing that it has a ${\mathcal{O}}(1/k)$ convergence rate on the objective function value. Similar results can be found in for forward–backward methods with linesearch along the projection arc (i.e. of the form (25) with $D_{k}=I$, $\lambda_{k}=1$ for all $k$ and with $\alpha_{k}$ determined by a backtracking procedure).

Assume that the hypotheses of Theorem 3.1 hold and, in addition, that assumption a) or b) of Proposition 2.2 is satisfied. Let $f^{*}$ be the optimal function value for problem (1). Then, we have

Proof. Setting $a=2\lambda_{\min}\alpha_{\min}$, where $\lambda_{\min}$ is defined in Proposition 2.2, from (3) we have

where the second inequality follows from the fact that $\nabla f(\boldsymbol{x}^{(k)})^{T}(\boldsymbol{y}^{(k)}-\boldsymbol{x}^{(k)})$ and $f(\hat{\boldsymbol{x}})-f(\boldsymbol{x}^{(k)})$ are negative quantities. Thanks to inequality (4), we can write

By multiplying the last inequality by $\mu_{k}$ we obtain

where the last inequality follows from the fact that $\mu_{k}\geq 1$. By repeatedly applying the last inequality we obtain

where, as in the proof of Theorem 3.1, we set $\theta^{k}_{j}=\prod_{i=j}^{k}\mu_{j}^{2}$ and $M$ is the upper bound of all $\theta^{k}_{j}$. Thanks to inequality (13), we have

where we also added the positive quantity $a(f({\boldsymbol{x}}^{(0)})-f(\hat{\boldsymbol{x}}))$ to the right hand side of (28). Moreover, exploiting the inequality

establishing the result. $\square$ In the recent literature, several authors developed the so-called intertial methods, which are first order methods including an extrapolation step which allows to prove a ${\mathcal{O}}(1/{k^{2}})$ convergence rate on the objective function values (see for example ). However, as we will show in section 4, the practical performances of SGP can be comparable with those of ${\mathcal{O}}(1/{k^{2}})$ methods, even if the theoretical convergence rate estimate is only ${\mathcal{O}}(1/k)$.

Numerical illustration

In this section we consider some relevant applications and we show that they can be effectively solved by algorithms which can be framed in the analysis of the previous sections. We give also some hints on how to choose the parameters $\alpha_{k}$ and $D_{k}$ at each iteration, even if a specific treatment of this issue is far beyond the scope of this paper. Both sets of numerical tests concern the image deconvolution problem in the presence of Poisson noise. In particular, in the next subsection we will consider a fit-to-data + regularization model with an arbitrarily fixed regularization parameter, while in the following tests we will investigate the same problem combined with an automatic procedure for the choice of this parameter recently proposed by Zanni et al. .

where $KL({\boldsymbol{x}})$ is the generalized Kullback–Leibler divergence

$R({\boldsymbol{x}})$ is some regularization functional, chosen according to the a priori information on the desired solution, and $\nu>0$ is the regularization parameter balancing the relative weight of the two terms. In order to preserve the edges in the restored image, a good choice for the regularization term is the following hypersurface (HS) functional

where $\rho>0$ and $\boldsymbol{D}_{i}^{h}$, $\boldsymbol{D}_{i}^{v}$ are finite difference approximations of the horizontal and vertical image gradient, respectively. If $\rho$ is small, it can be considered as an approximation of the total variation functional, but it has been shown that better reconstructions can be obtained for large values of the smoothing parameter . Thus, we consider the following convex optimization problem

whose main features have been studied in . As for the SGP method, borrowing the ideas in , at each iteration $k$ we adopt the following diagonal scaling matrix

where $V^{R}_{i}(\boldsymbol{x}^{(k)})$ is defined as in [18, formula (25)], while $\mu_{k}=\sqrt{1+10^{10}/k^{2}}$ so that Theorem 3.1 applies. The steplength parameter $\alpha_{k}$ is then computed in two different ways:

the adaptive alternation of the scaled Barzilai–Borwein (BB) rules as proposed in ;

the Ritz-like values proposed by Fletcher for a steepest descent method in the case of unconstrained optimization and recently extended to the SGP algorithm applied to a general constrained problem (1) .

Besides SGP, we consider also for comparison the “plain” gradient projection (GP) method with Euclidean projection and variable steplength (chosen with the same two rules exploited in the scaled case), the PidSplit+ algorithm , which is an alternating direction method of multipliers specific for the minimization of the Kullback–Leibler plus the discrete total variation functional ($\rho=0$), adapted to the smoothed case with $\rho>0$, and the accelerated proximal-gradient method with inertial/extrapolation with backtracking (FISTA-b) . As test problems, we consider:

the Shepp-Logan (SL) phantom of size $256\times 256$, multiplied by a factor of 500, corrupted with Gaussian blur of variance 9 and with Poisson noise simulated using the imnoise Matlab function on the blurred image including the additive background. The background constant is $b=10$;

the confocal microscopy (CM) phantom of size $128\times 128$ described in [44, section V.C], with values in the range $$and with a constant background$b=1$.

We assume periodic boundary conditions, so that the matrix $A$ is block circulant with circulant blocks (BCCB) and the matrix-vector products involving $A$ can be performed with a $\mathcal{O}(n\log(n))$ complexity by means of the fast Fourier transform . In figure 1 we report the original objects, the corrupted images and the solutions $\boldsymbol{x}^{*}$ of problem (30) for both test problems. The parameters $(\nu,\rho)$ in (33) have been empirically tuned to obtain a visually satisfactory solution and have been set equal to $(0.0415,1)$ for SL and $(0.06,1)$ for CM. Moreover, the ‘$\gamma$’ parameter of PidSplit+ has been set equal to $50/\nu$ (SL) and $1/\nu$ (CM) and the initial steplength parameter for FISTA-b is 100 in both cases.

In order to illustrate the convergence behaviour of the methods, we first compute a ground truth solution ${\boldsymbol{x}}_{\nu}$ (see figure 1, right panel) by running 1500 iterations of SGP. Then, we evaluate the progress towards this solution by computing at each iterate the relative difference of the objective function value with respect to the estimated minimum $f({\boldsymbol{x}}_{\nu})$ (see figure 2). We include in our comparison also the version of SGP with fixed bounds on the scaling matrix $\mu_{k}=\mu=10^{5}$, which is denoted by SGP∗.

From figure 2 we can observe what follows:

the choice of a suitable projection operator can have a significant impact on the practical performances of the gradient projection methods, since GP is outperformed by SGP with both choices for the steplength parameters;

SGP with variable bounds on the scaling matrix gives the best performances: in particular, condition (18), which is employed in Theorem 3.1 to prove the convergence of the method on convex problems, seems also to significantly improve its practical performances, especially when the iterates are close to the solution;

in spite of the theoretical convergence rate given in Theorem 3.2, the practical behaviour of SGP is comparable with the ${\mathcal{O}}(1/{k^{2}})$ method FISTA-b.

2 Automatic parameter estimation

The choice of the regularization parameter in Poisson data inversion is an active field and several different strategies have been proposed in the last years . Here we consider that proposed by Bertero et al. , which consists of selecting the value of $\nu$ in (30) such that

where $\eta$ is a given number close to 1 (here we will assume $\eta=1$). In particular, in the authors introduced an effective secant-type solver for the discrepancy equation (35), called modified Dai-Fletcher (MDF) method, able to reduce the number of required solutions of problems (30). At each step of the secant method, an approximation of the solution of problem (30) for a given value of $\nu$ is provided by running an optimization method until the stopping criterium

where $\varepsilon=5\times 10^{-8}$, is satisfied or when a maximum number of iterations equal to 5000 is reached. In this section we consider again the KL + HS model (33) and we investigate the impact of (some of) the strategies used for the previous tests within this automatic scheme for the choice of $\nu$. In particular, we restrict our analysis to the GP, SGP∗ and SGP methods equipped with the Ritz-like steplengths, and the PidSplit+ algorithm with the adaptive choice of its parameter $\gamma$ described in [35, equation (24)], which resulted to be less dependent on the parameter settings than the standard approach. The test problems we considered are based on the Satellite dataset already used in several papers and available at www.mathcs.emory.edu/$\sim$nagy/RestoreTools/index.html. The original image is sized $256\times 256$ and assumes values in the range $$. The blurred image has been obtained by convolving the object with a point spread function simulating a ground-based telescope response, and a constant background$b=10$has been added to the resulting image before introducing Poisson noise. Two further datasets have been obtained by multiplying object and background by factors of 10 and 100 before the blurring step. The three test sets will be denoted by S2550, S25500 and S255000 and the corrupted images are shown in figure 3 together with the original one. As concerns the parameter$\rho$defining the HS regularization term, we followed the suggestion in and set$\rho=10^{-4}\max(\boldsymbol{g})$.

The results obtained by the algorithms are shown in table 1, where we reported the number of steps $k$ of the secant-based method required to satisfy either the relation

being $\varepsilon_{1}=5\times 10^{-4}$ and $\varepsilon_{2}=5\times 10^{-3}$, the total number of iterations $k_{\rm{tot}}$ performed by each method in the $k$ steps, the final regularization parameter $\nu_{k}$, the relative reconstruction error between ${\boldsymbol{x}}_{\nu_{k}}$ and $\boldsymbol{x}^{*}$ and the execution time in seconds. These numerical experiments has been carried out on a Dual CPU Intel(R) Xeon(R) X5690 at 3.47GHz with 188 GB RAM (see also fermi.unife.it) in a Matlab2013a environment.

The performances summarized in table 1 confirm what already observed in the previous section, since SGP equipped with the scaling matrices with variable bounds succeeds in reducing the overall number of iterations required to provide the regularization parameter and the corresponding reconstruction if compared with SGP with fixed bounds for the scaling matrices or GP (which, in two of the three tests, often fails in satisfying the stopping criterium (36) within the maximum number of iterations allowed). As concerns the comparison with PidSplit+, we can observe that the number of iterations performed by this latter strategy is lower than that of SGP, but the higher cost per iteration which characterizes PidSplit+ makes the procedure more expensive in terms of total CPU time with respect to the SGP method.

Conclusions

In this paper we revisited the SGP method, originally published in 2009 and exploited in the successive years in several inverse problems as image denoising/deblurring, Fourier-based image reconstruction, blind deconvolution, system identification and non-negative matrix factorization, with several applications in astronomy, microscopy and engineering. Despite all the good numerical results provided in solving these problems, the only theoretical convergence result proved so far is the stationarity of any limit point of the sequence generated by SGP. In this paper we showed that stronger results can be proved in the convex case, if the sequence of scaling matrices characterizing the SGP iterations is chosen as convergent to the identity matrix at a certain rate. Moreover, in the same setting we provided also a convergence rate estimate on the objective function values, as provided in the literature for several other optimization methods. Some numerical tests showed that the specific rule introduced on the scaling matrices to prove the theoretical convergence results helps also to improve the performances of the method, making SGP competitive also with methods for which the ${\mathcal{O}}(1/{k^{2}})$ convergence rate has been demonstrated.

Acknowledgments

This work has been partially supported by MIUR (Italian Ministry for University and Research), under the projects FIRB - Futuro in Ricerca 2012, contract RBFR12M3AC, and PRIN 2012, contract 2012MTE38N. The Italian GNCS - INdAM (Gruppo Nazionale per il Calcolo Scientifico - Istituto Nazionale di Alta Matematica) is also acknowledged.

References

References