Rahul Mourya

Loïc Denis

Jean-Marie Becker

Eric Thiébaut

AUGMENTED LAGRANGIAN WITHOUT ALTERNATING DIRECTIONS: PRACTICAL ALGORITHMS FOR INVERSE PROBLEMS IN IMAGING

Keywords: Constrained convex optimization, nonsmooth optimization, hierarchical optimization, ADMM, proximal operator, deblurring, total variation

Several problems in signal processing and machine learning can be casted as optimization problems. In many cases, they are of large-scale, nonlinear, have constraints, and may be nonsmooth in the unknown parameters. There exists plethora of fast algorithms for smooth convex optimization, but these algorithms are not readily applicable to nonsmooth problems, which has led to a considerable amount of research in this direction. In this paper, we propose a general algorithm for nonsmooth bound-constrained convex optimization problems. Our algorithm is instance of the so-called augmented Lagrangian, for which theoretical convergence is well established for convex problems. The proposed algorithm is a blend of superlinearly convergent limited memory quasi-Newton method, and proximal projection operator. The initial promising numerical results for total-variation based image deblurring show that they are as fast as the best existing algorithms in the same class, but with fewer and less sensitive tuning parameters, which makes a huge difference in practice.

INTRODUCTION

Many problems in signal processing and machine learning can be stated follows:

x * := arg min

x F (x) subject to x ∈ Ω, (1)
where F (x) = {f (x) + r(x)}; f : R n → R is twice continuously differentiable, and convex; r : R n → R is continuous and convex but not necessarily differentiable everywhere, and Ω is a bound constraint set. Often, f corresponds to a loss, and r to a regularizer. This class of problems, in general, do not have close-form solutions, and rely on iterative method. There exists numerous well established algorithms [START_REF] Zhu | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF] (LBFGS-B), [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] (VMLM-B), [START_REF] Steven | A Limited Memory Variable Metric Method in Subspaces and Bound Constrained Optimization Problems[END_REF] (BLMVM), [START_REF] Hager | A New Active Set Algorithm for Box Constrained Optimization[END_REF] (ASA-CG), [START_REF] Mark W Schmidt | Optimizing costly functions with simple constraints: A limited-memory projected quasi-newton algorithm[END_REF] (minConf-TMP), based on limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] with superlinear convergence rate for the problem (1) when ∀x, r(x) = 0, but these techniques do not work off-the-shelf when ∀x, r(x) = 0. Some efforts have been made in an adhoc manner to apply LBFGS methods directly to nonsmooth convex problems that are differentiable almost everywhere, and convergence to the optimum has been noted [START_REF] Lemaréchal | Numerical Experiments in Nonsmooth Optimization[END_REF] in the cases where no nonsmooth point is encountered, otherwise [START_REF] Luksan | Globally Convergent Variable Metric Method for Convex Nonsmooth Unconstrained Minimization[END_REF][START_REF] Haarala | New limited memory bundle method for large-scale nonsmooth optimization[END_REF][START_REF] Lewis | Nonsmooth optimization via quasi-Newton methods[END_REF] report catastrophic failures (convergence to a nonoptimum) of the direct methods. The traditional algorithms for nonsmooth optimization are based on a stabilization steepest descent by exploiting gradient or subgradient information evaluated at multiple points, which is the essential idea behind subgradient methods [START_REF] Nedic | Convergence Rate of Incremental Subgradient Algorithms[END_REF][START_REF] Yu | A Quasi-Newton Approach to Nonsmooth Convex Optimization Problems in Machine Learning[END_REF], the bundle methods [START_REF] Haarala | New limited memory bundle method for large-scale nonsmooth optimization[END_REF][START_REF] Hui Teo | Bundle Methods for Regularized Risk Minimization[END_REF], and the gradient sampling algorithm [START_REF] Burke | A Robust Gradient Sampling Algorithm for Nonsmooth, Nonconvex Optimization[END_REF][START_REF] Krzysztof | Convergence of the Gradient Sampling Algorithm for Nonsmooth Nonconvex Optimization[END_REF]. These algorithms are successful for nonsmooth problems with sublinear rate of convergence.

Another class of algorithms for solving the problem (1) use the general idea of variable splitting: introduction of auxiliary variables in order to decouple minimizations of f and of r. The proximal forward-backward iterative scheme introduced in [START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in Hilbert space[END_REF] and [START_REF] Bruck | On the Weak Convergence of an Ergodic Iteration for the Solution of Variational Inequalities for Monontone Operators in Hilbert Space[END_REF] is a notable representative algorithm of this class; see the recent survey [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] and the references therein for very general convergence results of proximal forward-backward algorithms under various conditions and settings. relevant to problem [START_REF] Zhu | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF]. Unlike the algorithms using second-order information (quasi-Newton methods), the proximal forward-backward algorithms use only first order information, and are slower having only sublinear convergence rate, but enjoy global convergence guarantees for nonsmooth convex problems. The convergence rates of these algorithms have been further improved in [START_REF] José | A New TwIST: Two-Step Iterative Shrinkage/Thresholding Algorithms for Image Restoration[END_REF][START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF] by using the idea of Nesterov [START_REF] Nesterov | A method for unconstrained convex minimization problem with the rate of convergence O (1/k2)[END_REF] developed in 1983, which suggest to use the information from the previous iterations in a smart way at each new iteration.

The two instances of the general proximal forwardbackward algorithm are Alternating Minimization Algorithm (AMA) [START_REF] Tseng | Applications of a Splitting Algorithm to Decomposition in Convex Programming and Variational Inequalities[END_REF], and closely related Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF] whose developments backs to 1970s. These algorithms use augmented Lagrangian (AL) terms to handle the constraints, and are closely related to algorithms such as dual decomposition, the method of multipliers, Douglas-Rachford splittings, Dykstra's alternating projections, Bregman iterative algorithm, and others. ADMM is a blend of decomposability of dual ascent and superior convergence properties of the method of multipliers. Because of its flexibility to in handling different types of objective functions/constraints, and its simplicity in implementation for distributed optimization, ADMM has gained popularity in both signal processing and machine learning communities since the last two decades. However, only sublinear rate of convergence have been achievable by the ADMM, and convergence speed is highly dependent upon the chosen penalty parameters. Optimal tuning of those parameters remains largely an open question. In this paper we propose an algorithm for solving nonsmooth bound-constrained convex optimization problem, which uses variable splitting and AL to handle the nonsmooth part, and the robust quasi-Netwon method to handle the bound-constraint and smooth part. The proposed algorithm is as fast as the other algorithms in the same class with fewer number of tuning parameters, and is immune to large range variation in the tuning parameter.

PROPOSED ALGORITHM

The problem (1) without the bound constraint can be written equivalently after variable splitting as: arg min

x, z

{ f (x) + r(z) } s.t. x -z = 0. (2)
The Augmented Lagrangian of problem (2) can be written:

L ρ (x, z, u) = f (x) + r(z) + ρ 2 x -z + u 2 2 , (3)
where u is dual scaled variable, and ρ > 0 is a augmented penalty parameter.

The method of multipliers solves the problem (2) with bound constraint, (x ∈ Ω), by the following iterations (k being the iteration counter):

{x (k+1) , z (k+1) } := arg min z,x∈Ω { f (x) + r(z) + ρ 2 x -z + u (k) 2 2 } (4) u (k+1) := u (k) + x (k+1) -z (k+1)
(5) Rather than jointly minimizing with respect to x and z variables, ADMM solves the same problem by following iterations:

x

(k+1) := arg min x∈Ω { f (x) + ρ 2 x -z (k) + u (k) 2 2 } (6) z (k+1) := arg min z { r(z) + ρ 2 x (k+1) -z + u (k) 2 2 } (7) u (k+1) := u (k) + x (k+1) -z (k+1) (8)
The above two methods for solving the problem (2) are very similar. ADMM is viewed as a version of the method of multipliers with a single Gauss-Seidel iteration over x and z. The real advantage of these methods is their convergence to the solution under very general conditions [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF]: (i) the extended-real-valued functions f and r in problem (2) are closed, proper, and convex; and (ii) the unaugmented Lagrangian L 0 has a saddle point. Moreover, the two methods converge even when x-, z-updates are not carried out exactly, provided that the errors between the approximate and exact solution at each iteration are summable [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF].

Motivation

In general the global convergence rate of ADMM is limited to sublinear rate, however, recently numerous algorithms [START_REF] Afonos | Fast Image Recovery Using Variable Splitting and Constrained Optimization[END_REF][START_REF] Manya V Afonso | An Augmented Lagrangian Approach to the Constrained Optimization Formulation of Imaging Inverse Problems[END_REF][START_REF] Matakos | Accelerated Edge-Preserving Image Restoration Without Boundary Artifacts[END_REF] based on ADMM for solving the instances of the general problem [START_REF] Zhu | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF] in applications related to image reconstruction and other linear inverse problems have proven to be more efficient, and converges faster in comparison to many stateof-the-art algorithms [START_REF] Wright | Sparse Reconstruction by Separable Approximation[END_REF][START_REF] José | A New TwIST: Two-Step Iterative Shrinkage/Thresholding Algorithms for Image Restoration[END_REF][START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF]. The convergence speed of ADMM is highly dependent upon the penalty parameters, and only optimally tunned parameters result in high convergence speed. The optimal tunning of the penalty parameters is still an open problem. It has also been observed that the convergence speed is higher when separate penalty parameters are used for each augmented term introduced in the ADMM. The experimental results in [START_REF] Matakos | Accelerated Edge-Preserving Image Restoration Without Boundary Artifacts[END_REF] clearly demonstrate that the ADMM with i). sufficient number of variable splittings (so that each variable updates can be carried out in closed-form), and ii). optimally tunned penalty parameters, converges faster than the other variants of ADMM. Following this suggestion, a nonsmooth convex problem such as image deblurring problem with total variation(TV) regularization, the one considered in section (3), requires multiple variable splittings, which requires to tune multiple parameters hindering the convergence speed. In this paper, we consider going back to the method of multipliers and avoiding as many splitting as possible so as to reduce the number of penalty parameters. We show that by using a hierarchical optimization approach, an efficient algorithm is obtained.

Derivation of the algorithm

Going back to the method of multipliers (4)-( 5), the joint minimization problem can be formulated in a hierarchical way:

x * := arg min

x∈Ω f (x) + r(z * (x)) + ρ 2 x -z * (x) + u 2 2 with z * (x) = arg min z r(z) + ρ 2 x -z + u 2 2 (9)
The inner optimization problem is often known in closed form. We use a quasi-Newton method to compute the solution of the outer minimization problem.

Proposition 1. The gradient of the partially optimized augmented Lagrangian with respect to x is given by:

∇ x L ρ (x, z * (x), u) = ∇f (x) + ρ(x -z * (x) + u) . ( 10
)
Proof. When r is a smooth function, this follows directly from the chain rule formula since the gradient of augmented Lagrangian with respect to z is zero at the optimal value z * . In the non-smooth case, the augmented Lagrangian is still differentiable by virtue of the augmented term that smooths r [START_REF] Parikh | Proximal Algorithms[END_REF].

In order to prove our proposition in the non-smooth case, we use a classical characterization of the gradient of a convex function, namely that the function lies above its gradient. As z * (x) is the solution of the inner minimization:

∀z, ∀x, r(z) + ρ 2 x -z + u 2 2 ≥ r(z * (x)) + ρ 2 x -z * (x) + u 2 2 .
(11) For a fixed z, the gradient of the augmented Lagrangian L ρ (x, z, u) with respect to x is:

∇ x L ρ (x, z, u) = ∇f (x) + ρ(x -z + u) . (12
) By convexity of f , we have: ∀x, ∀z, L ρ (x, z, u) ≥ L ρ (x, z, u) 13) Equation ( 11) applied to z * (x) and x and equation ( 13) applied to z * (x) lead to: ∀x, L ρ (x, z * (x), u) ≥ L ρ (x, z * (x), u)

+ [∇f (x) + ρ(x -z + u)] t (x -x) . (
+ [∇f (x) + ρ(x -z * (x) + u)] t (x -x) . ( 14
)
This inequality is close to the gradient inequality sought, except that z * (x) appears instead of z * (x) in the second line. We thus need to prove that (z * (x)z * (x)) t (xx) ≥ 0 for all x. By expanding the squared norm in equation ( 11), we can show that ∀z, ∀x, r(z)

+ ρ 2 z 2 2 -ρz t (x + u) ≥ r(z * (x)) + ρ 2 z * (x) 2 2 -ρz * (x) t (x + u) . ( 15
)
Combining equation [START_REF] Krzysztof | Convergence of the Gradient Sampling Algorithm for Nonsmooth Nonconvex Optimization[END_REF] applied to (z := z * (x), x := x) and to (z := z * (x), x) gives the desired result. In conclusion, the following inequality holds: ∀x, L ρ (x, z * (x), u) ≥ L ρ (x, z * (x), u)

+ [∇f (x) + ρ(x -z * (x) + u)] t (x -x) , (16)
which proves the proposition in Eq.( 10).

Proposition 2. The gradient difference is independent of the value of dual variable u: 1) , z * (x (1) ), u = ∇f (x (2) )-∇f (x (1) )+ρ x (2) -x (1) +z * (x (1) )-z * (x (2) )

∇ x L ρ x (2) , z * (x (2) ), u -∇ x L ρ x (
Using proposition 1, any smooth optimization method based solely on cost function and gradient evaluations can be used to perform joint minimization over x and z. If the method uses gradient differences to collect second order information (e.g., quasi-Newton methods), the memorized previous steps can be used even after the dual variables have been updated, since gradient differences are not affected by dual updates, see proposition 2. The applicability of smooth optimization methods is thus extended to non-smooth problems typically encountered in image reconstruction.

ILLUSTRATION ON IMAGE DEBLURRING

We evaluate the efficiency of our algorithm on an image deblurring problem. We consider maximum a posteriori deblurring with TV regularization and positivity constraints:

x * := arg min x≥0 1 2 y -Hx 2 W + λ TV(x) , (17)
where y ∈ R n is the observed (blurry and noisy) image, H ∈ R n×n is the blurring operator (discrete convolution matrix), W ∈ R n×n is the inverse of noise covariance matrix (diagonal in our case), λ > 0 is tuned to reach good compromise between smoothness of the solution and data fitting. TV(x) is the isotropic-TV on the unknown image, x, defined as:

TV(x) = n i=1 (∇x) i 2 , where ∇ = [∇ T x ∇ T y ]
T is first-order finite difference operator in two dimensions (see e.g., [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]). The v 2 W is defined as v t Wv, ∀v. As done in [START_REF] Matakos | Accelerated Edge-Preserving Image Restoration Without Boundary Artifacts[END_REF], to correctly handle borders, we reconstruct x that is larger than the available blurred image. Pixel values outside the field-of-view are given zero weight by matrix W (i.e., if pixel k is not seen, then W k,k = 0). We apply the deblurring model [START_REF] Bruck | On the Weak Convergence of an Ergodic Iteration for the Solution of Variational Inequalities for Monontone Operators in Hilbert Space[END_REF] on a portion of blurry and noisy Lena image shown in Figure 3. Note the zero-padding on the border of the observed image to handle it correctly.

Previous research [START_REF] Matakos | Accelerated Edge-Preserving Image Restoration Without Boundary Artifacts[END_REF] has shown the superiority of ADMM with multiple splittings over other state-of-theart methods. We compare the convergence speed of (i) ADMM (hereafter method ADMM 3x) with three splittings: y -Hx = ξ; ∇x = ω; z = x, and the resulting AL:

1 2 ξ 2 W + ρ 2 y -Hx -ξ + u 1 2 2 + g(z) + ν 2 x -z + u 2 2 2 + λ n i=1 ω i 2 + γ 2 n i=1 (∇x -ω + u 3 ) i 2 2 (18)
that makes possible to solve each sub-problems in closedform. g is an indicator function for the constraint z ≥ 0.

The three penalty parameters are tuned after many experimental trails; (ii) ADMM (hereafter ADMM 1x) with a single splitting: ∇x = ω with similar AL as in [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. The x-update is solved approximately by a quasi-Newton method (VMLM-B) and ω-update in closed form by the proximal operator (soft-thresholding); (iii) our proposed algorithm (hereafter proposed) with or without restarting the quasi-Newton method (VMLM-B) with the same variable splitting and AL as in ADMM 1x. All these methods are implemented in MATLAB, except the core of the VMLM-B (LBFGS update and line-search) is implemented in C language. In all of methods the most expensive operation is the circular convolution, which is done using Fast Fourier Transform (FFT).

Other operations have linear cost and are moreover same in all the methods, thus we neglect them in the convergence speed comparison. The Figure 1 compares convergence speed of the four methods against the number of FFTs consumed, and clearly the proposed method without restarting VMLM-B has speed comparable to the fastest ADMM 3x.

In ADMM 1x, and proposed methods, the number of maximum successful iterations in VMLM-B was fixed to 150 and 60 respectively. Table 1 reports the relative distance between the cost function after 1500 FFTs evaluations and the optimal value (estimated by the value reached after 10000 FFTs by the best performing algorithm). When the sub-minimizations are performed using an iterative method (VMLM-B), the number of Fig. 3:

On left, the observed image (241x241, BSNR=37.25dB); on top right corner, the blur kernel (15x15); on right, the estimated image (255x255, ISNR=6.12dB) after 1000 FFT evaluations by the proposed method. λ is tunned by hand to have possible best image quality. iterations influence the convergence speed. We therefore report different values corresponding to 30, 60 or 100 inner iterations. As observed in [START_REF] Matakos | Accelerated Edge-Preserving Image Restoration Without Boundary Artifacts[END_REF], doing several splittings and using well tuned penalty parameters is preferable than performing an approximate minimization of the more complex subproblem posed with a single splitting. Our proposed method reaches a better solution after the considered number of iterations if previous gradient steps are re-used after dual variable update (last row of the table). Loosing this second-order information noticeably impacts the performance.

Figure 2 illustrates the influence of the penalty parameters on the convergence speed of the four optimization methods, and the proposed methods are immune to large range of variation in the penalty parameter. ADMM 3x can converge very quickly for well tunned parameters, the convergence degrades strongly for a bad tuning. Given the number of parameters to be jointly tuned, this is a clear practical drawback. Reducing the number of splitting by resorting to iterative minimization of the sub-problems, as done with ADMM 1x, simplifies parameter tuning at the cost of a degraded convergence. Our proposed method offers a remarkable robustness to changes of the single penalty parameter while showing a convergence speed comparable to ADMM 3x with optimal tuning.

CONCLUSION

We proposed a general approach to handle non-smooth largescale optimization problems, which follows the augmented Lagrangian approach. We have shown that hierarchical optimization is preferable to performing an alternate minimization, as usually done with ADMM. Quasi-Newton methods can then be applied efficiently thanks to a simple expression of the gradient, and the possibility to accumulate gradient difference knowledge independently of dual variable updates. Our results on image deblurring problem are promising: competitive convergence speed is reached with a method that is much simpler to tune.

1 Fig. 1 : 1 Fig. 2 :

 1112 Fig. 1: Convergence speed comparison of the optimization methods with possible best penalty parameters (ρ, γ) for each.

Table 1 :

 1 Relative difference between the cost function after 1500 FFT evaluations and at the solution (in %)

	ADMM 3x	.35
	ADMM 1x	2.6 (30 it) 2.7 (60 it	2.2 (100 it)
	proposed(rst) 2.2 (30 it) .55 (60 it) .20 (100 it)
	proposed	.06 (30 it) .04 (60 it) .038 (100 it)

method

[F (x (1500) ) -F (x * )]/F (x * )

The author is supported by a PhD grant funded by the ARC6, Région Rhône-Alpes.