Sparsity Constrained Image Restoration: An Approach Using the Newton Projection Method

. Image restoration under sparsity constraints has received increased attention in recent years. This problem can be formulated as a nondiﬀerentiable convex optimization problem whose solution is challenging. In this work, the non-diﬀerentiability of the objective is addressed by reformulating the image restoration problem as a nonnega-tively constrained quadratic program which is then solved by a specialized Newton projection method where the search direction computation only requires matrix-vector operations. A comparative study with state-of-the-art methods is performed in order to illustrate the eﬃciency and eﬀectiveness of the proposed approach.


Introduction
This work is concerned with the general problem of image restoration under sparsity constraints formulated as where A is a m × n real matrix (usually m ≤ n), x ∈ R n , b ∈ R m and λ is a positive parameter.(Throughout the paper, • will denote the Euclidean norm).In some image restoration applications, A = KW where K ∈ R m×n is a discretized linear operator and W ∈ R n×n is a transformation matrix from a domain where the image is a priori known to have a sparse representation.The variable x contains the coefficients of the unknown image and the data b is the measurements vector which is assumed to be affected by Gaussian white noise intrinsic to the detection process.The formulation (1) is usually referred to as synthesis formulation since it is based on the synthesis equation of the unknown image from its coefficients x.
The penalization of the 1 -norm of the coefficients vector x in (1) simultaneously favors sparsity and avoids overfitting.For this reason, sparsity constrained image restoration has received considerable attention in the recent literature and has been successfully used in various areas.The efficient solution of problem (1) is a critical issue since the nondifferentiability of the 1 -norm makes standard unconstrained optimization methods unusable.Among the current state-of-theart methods there are gradient descent-type methods as TwIST [5], SparSA [14], FISTA [2] and NESTA [3].GPSR [6] is a gradient-projection algorithm for the equivalent convex quadratic program obtained by splitting the variable x in its positive and negative parts.Fixed-point continuation methods [9], as well as methods based on Bregman iterations [7] and variable splitting, as SALSA [1], have also been recently proposed.In [12], the classic Newton projection method is used to solve the bound-constrained quadratic program formulation of (1) obtained by splitting x.A Modified Newton projection (MNP) method has been recently proposed in [11] for the analysis formulation of the 1 -regularized least squares problem where W is the identity matrix and x represents the image itself.The MNP method uses a fair regularized approximation to the Hessian matrix so that products of its inverse and vectors can be computed at low computational cost.As a result, the only operations required for the search direction computation are matrix-vector products.
The main contribution of this work is to extend the MNP method of [11], developed for the case W = I n , to the synthesis formulation of problem (1) where W = I n .In the proposed approach, problem (1) is firstly formulated as a nonnegatively constrained quadratic programming problem by splitting the variable x into the positive and negative parts.Then, the quadratic program is solved by a special purpose MNP method where a fair regularized approximation to the Hessian matrix is proposed so that products of its inverse and vectors can be computed at low computational cost.As a result, the search direction can be efficiently obtained.The convergence of the proposed MNP method is analyzed.Even if the size of the problem is doubled, the low computational cost per iteration and less iterative steps make MNP quite efficient.The performance of MNP is evaluated on several image restoration problems and is compared with that of some state-of-the-art methods.The results of the comparative study show that MNP is competitive and in some cases is outperforming some state-of-theart methods in terms of computational complexity and achieved accuracy.
The rest of the paper can be outlined as follows.In section 2, the quadratic program formulation of (1) is derived.The MNP method is presented and its convergence is analyzed in section 3.In this section, the efficient computation of the search direction is also discussed.In section 4, the numerical results are presented.Conclusions are given in section 5.

Nonnegatively constrained quadratic program formulation
The proposed approach firstly needs to reformulate (1) as a nonnegatively constrained quadratic program (NCQP).The NCQP formulation is obtained by splitting the variable x into its positive and negative parts [6], i.e Problem (1) can be written as the following NCQP: where 1 denotes the n-dimensional column vector of ones.The gradient g and Hessian H of F(u, v) are respectively defined by We remark that the computation of the objective function and its gradient values requires only one multiplication by A and one by A H , nevertheless the double of the problem size.Since H is positive semidefinite, we propose to approximate it with the positive definite matrix H τ : where τ is a positive parameter and I is the identity matrix of size n.
Proposition 21 Let σ 1 , σ 2 , . . ., σ n be the nonnegative eigenvalues of A in nonincreasing order: Then, H τ is a positive definite matrix whose eigenvalues are The proof is immediate since the spectrum of H τ is the union of the spectra of The following proposition shows that an explicit formula for the inverse of H τ can be derived.

Proposition 22
The inverse of the matrix H τ is the matrix M τ defined as where Proof.We have Similarly, we can prove that We now show that M 1 M 2 = M 2 M 1 .We have After some simple algebra, it can be proved that and thus From ( 8), ( 9) and ( 10), it follows 3 The Modified Newton projection method The algorithm.The Newton projection method [4] for problem (2) can be written as where [•] + denotes the projection on the positive orthant, g u and g v respectively indicate the partial derivatives of F with respect to u and v at the current iterate.The scaling matrix S (k) is a partially diagonal matrix with respect to the index set A (k) defined as and ε is a small positive parameter.The step-length α (k) is computed with the Armijo rule along the projection arc [4].Let E (k) and F (k) be the diagonal matrices [13] such that In MNP, we propose to define the scaling matrix S (k) as Therefore, the complexity of computation of the search direction is mainly due to one multiplication of the inverse Hessian approximation M τ with a vector since matrix-vector products involving E (k) and F (k) only extracts some components of the vector and do not need not to be explicitly performed.
We remark that, in the Newton projection method proposed in [4], the scaling matrix is which requires to extract a submatrix of H τ and then to invert it.
Convergence analysis.As proved in [4], the convergence of Newton-type projection methods can be proved under the general following assumptions which basically require the scaling matrices S (k) to be positive definite matrices with uniformly bounded eigenvalues.
A1 The gradient g is Lipschitz continuous on each bounded set of R 2n .A2 There exist positive scalars c 1 and c 2 such that The key convergence result is provided in Proposition 2 of [4] which is restated here for the shake of completeness.
Proposition 31 [4, Proposition 2] Let {[u (k) , v (k) ]} be a sequence generated by iteration (11) where S (k) is a positive definite symmetric matrix which is diagonal with respect to A (k) and α k is computed by the Armijo rule along the projection arc.Under assumptions A1 and A2 above, every limit point of a sequence {[u (k) , v (k) ]} is a critical point with respect to problem (2).
Since the objective F of ( 2) is twice continuously differentiable, it satisfies assumption A1.From propositions 21 and 22, it follows that M τ is a symmetric positive definite matrix and hence, the scaling matrix S (k) defined by ( 12) is a positive definite symmetric matrix which is diagonal with respect to A (k) .The global convergence of the MNP method is therefore guaranteed provided S (k) verifies assumption A2.
Proposition 32 Let S (k) be the scaling matrix defined as k) .Then, there exist two positive scalars c 1 and c 2 such that Proof.Proposition 21 implies that the largest and smallest eigenvalue of M τ are respectively 1/τ and 1/(2σ 1 + τ ), therefore We have y H S (k) y = E (k) y H M τ E (k) y + y H F (k) y.From (13) it follows that The thesis follows by setting c 1 = min{ 1 2σ1+τ , 1} and c 2 = max{ 1 τ , 1}.
Computing the search direction.We suppose that K is the matrix representation of a spatially invariant convolution operator with periodic boundary conditions so that K is a block circulant with circulant blocks (BCCB) matrix and matrix-vector products can be efficiently performed via the FFT.Moreover, we assume that the columns of W form an orthogonal basis for which fast sparsifying algorithms exist, such as a wavelet basis, for example.Under these assumptions, A is a full and dense matrix but the computational cost of matrixvector operations with A and A H is relatively cheap.As shown by (12), the computation of the search direction p (k) = S (k) g (k) requires the multiplication of a vector by M τ .Let z, w ∈ R 2n be a given vector, then it immediately follows that Formula ( 14) needs the inversion of 2A H A + τ I. Our experimental results indicate that the search direction can be efficiently and effectively computed as follows.Using the Sherman-Morrison-Woodbury formula, we obtain Substituting ( 15) and ( 16) in ( 14), we obtain Since K is BCCB, it is diagonalized by the Discrete Fourier Transform (DFT), i.e.K = U H DU where U denotes the unitary matrix representing the DFT and D is a diagonal matrix.Thus, we have: Substituting ( 19) and (20) in ( 17) and (18), we obtain Equations ( 14), ( 21) and ( 22) show that, at each iteration, the computation of the search direction p (k) requires two products by W, two products by W H , two products by U and two products by U H .The last products can be performed efficiently by using the FFT algorithm.

Numerical results
In this section, we present the numerical results of some image restoration test problems.The numerical experiments aim at illustrating the performance of MNP compared with some state-of-the-art methods as SALSA [1], CGIST [8], the nonmonotonic version of GPSR [6], and the Split Bregman method [7].Even if SALSA has been shown to outperform GPSR [1], we consider GPSR in our comparative study since it solves, as MNP, the quadratic program (2).The Matlab source code of the considered methods, made publicly available by the authors, has been used in the numerical experiments.The numerical experiments have been executed in Matlab R2012a on a personal computer with an Intel Core i7-2600, 3.40 GHz processor.
The numerical experiments are based on the well-known Barbara image (figure 1), whose size is 512 × 512 and whose pixels have been scaled into the range between 0 and 1.In our experiments, the matrix W represents an orthogonal Haar wavelet transform with four levels.For all the considered methods, the initial iterate x (0) has been chosen as x (0) = W H b; the regularization parameter λ has been heuristically chosen.In MNP, the parameter τ of the Hessian approximation has been fixed at τ = 100λ.This value has been fixed after a wide experimentation and has been used in all the presented numerical experiments.
The methods iteration is terminated when the relative distance between two successive objective values becomes less than a tolerance tol φ .A maximum number of 100 iterations has been allowed for each method.
In the first experiment, the Barbara image has been convolved with a Gaussian PSF with variance equal to 2, obtained with the code psfGauss from [10], and then, the blurred image has been corrupted by Gaussian noise with noise level equal to 7.5 • 10 −3 and 1.5 • 10 −2 .(The noise level NL is defined as NL := η Ax original where x original is the original image and η is the noise vector.)In the second experiment, the Barbara image has been corrupted by out-of-focus blur, obtained with the code psfDefocus from [10] and by Gaussian noise with noise level equal to 2.5 • 10 −3 and 7.5 • 10 −3 .The degraded images and the MNP restorations are shown in figure 1 for NL = 7.5 • 10 −3 .Table 1 reports the Mean Squared Error (MSE) values, the objective function values and the CPU times in seconds obtained by using the stopping tolerance values tol φ = 10 −1 , 10 −2 , 10 −3 .In figure 2, the MSE behavior and the decreasing of the objective function versus time (in seconds) are illustrated.
The reported numerical results indicate that MNP is competitive with the considered state-of-the-art and, in terms of MSE reduction, MNP reaches the minimum MSE value very early.

Conclusions
In this work, the MNP method has been proposed for sparsity constrained image restoration.In order to gain low computational complexity, the MNP method uses a fair approximation of the Hessian matrix so that the search direction can be computed efficiently by only using FFTs and fast sparsifying algorithms.The results of numerical experiments show that MNP may be competitive with some state-of-the-art methods both in terms of computational efficiency and accuracy.