A quasi-Newton proximal splitting method

Introduction

Convex optimization has proved to be extremely useful to all quantitative disciplines of science. A common trend in modern science is the increase in size of datasets, which drives the need for more efficient optimization schemes. For large-scale unconstrained smooth convex problems, two classes of methods have seen the most success: limited memory quasi-Newton methods and non-linear conjugate gradient (CG) methods. Both of these methods generally outperform simpler methods, such as gradient descent.

For problems with non-smooth terms and/or constraints, it is possible to generalize gradient descent with proximal gradient descent (which includes projected gradient descent as a sub-cases), which is just the application of the forward-backward algorithm .

Unlike gradient descent, it is not easy to adapt quasi-Newton and CG methods to problems involving constraints and non-smooth terms. Much work has been written on the topic, and approaches generally follow an active-set methodology. In the limit, as the active-set is correctly identified, the methods behave similar to their unconstrained counterparts. These methods have seen success, but are not as efficient or as elegant as the unconstrained versions. In particular, a sub-problem on the active-set must be solved, and the accuracy of this sub-iteration must be tuned with heuristics in order to obtain competitive results.

Our goal is the generic minimization of functions of the form

The class we consider covers non-smooth convex optimization problems, including those with convex constraints. Here are some examples in regression, machine learning and classification.

One would like to find a linear decision function which minimizes the objective

2 Contributions

This paper introduces a class of scaled norms for which we can compute a proximity operator; these results themselves are significant, for previous results only cover diagonal scaling (the diagonal scaling result is trivial). Then, motivated by the discrepancy between constrained and unconstrained performance, we define a class of limited-memory quasi-Newton methods to solve ( P ) and that extends naturally and elegantly from the unconstrained to the constrained case. Most well-known quasi-Newton methods for constrained problems, such as L-BFGS-B , are only applicable to box constraints $l\leq x\leq u$. The power of our approach is that it applies to a wide-variety of useful non-smooth functionals (see §3.1.4 for a list) and that it does not rely on an active-set strategy. The approach uses the zero-memory SR1 algorithm, and we provide evidence that the non-diagonal term provides significant improvements over diagonal Hessians.

Quasi-Newton forward-backward splitting

In the following, define the quadratic approximation

The standard (non relaxed) version of the forward-backward splitting algorithm (also known as proximal or projected gradient descent) to solve ( P ) updates to a new iterate $x_{k+1}$ according to

Our diagonal+rank 1 quasi-Newton forward-backward splitting algorithm is listed in Algorithm 1 (with details for the quasi-Newton update in Algorithm 2, see §4 for details). These algorithms are listed as simply as possible to emphasize their important components; the actual software used for numerical tests is open-source and available at http://www.greyc.ensicaen.fr/~jfadili/software.html.

2 Relation to prior work

The algorithm in (5) is variously known as proximal descent or iterated shrinkage/thresholding algorithm (IST or ISTA). It has a grounded convergence theory, and also admits over-relaxation factors $\alpha\in(0,1)$ .

The spectral projected gradient (SPG) method was designed as an extension of the Barzilai-Borwein spectral step-length method to constrained problems. In , it was extended to non-smooth problems by allowing general proximity operators; we refer to this as SPG/SpaRSA (N.B. we do not use the SpaRSA implementation since we do not use warm-starts or restarts, in order to be fair to all algorithms). The Barzilai-Borwein method use a specific choice of step-length $t_{k}$ motivated by quasi-Newton methods. Numerical evidence suggests the SPG/SpaRSA method is highly effective, although convergence results are not as strong as for ISTA.

FISTA is a multi-step accelerated version of ISTA inspired by the work of Nesterov. The stepsize $t$ is chosen in a similar way to ISTA; in our implementation, we tweak the original approach by using a Barzilai-Borwein step size, a standard line search, and restart, since this led to improved performance. Nesterov acceleration can be viewed as an over-relaxed version of ISTA with a specific, non-constant over-relaxation parameter $\alpha_{k}$.

The above approaches assume $B_{k}$ is a constant diagonal. The general diagonal case was considered in several papers in the 1980s as a simple quasi-Newton method, but never widely adapted. More recent attempts include a static choice $B_{k}\equiv B$ for a primal-dual method . A convergence rate analysis of forward-backward splitting with static and variable $B_{k}$ where one of the operators is maximal strongly monotone is given in .

By transforming the problem into a standard conic programming problem, the generic problem is amenable to interior-point methods (IPM). IPM requires solving a Newton-step equation, so first-order like “Hessian-free” variants of IPM solve the Newton-step approximately, either by approximately solving the equation or by subsampling the Hessian. The main issues are speed and robust stopping criteria for the approximations.

Yet another approach is to include the non-smooth $h$ term in the quadratic approximation. Yu et al. propose a non-smooth modification of BFGS and L-BFGS, and test on problems where $h$ is typically a hinge-loss or related function.

The projected quasi-Newton (PQN) algorithm is perhaps the most elegant and logical extension of quasi-Newton methods, but it involves solving a sub-iteration. PQN proposes the SPG algorithm for the subproblems, and finds that this is an efficient tradeoff whenever the cost function (which is not involved in the sub-iteration) is relatively much more expensive to evaluate than projecting onto the constraints. Again, the cost of the sub-problem solver (and a suitable stopping criteria for this inner solve) are issues. As discussed in , it is possible to generalize PQN to general non-smooth problems whenever the proximity operator is known (since, as mentioned above, it is possible to extend SPG to this case).

Proximity operators and proximal calculus

We only recall essential definitions. More notions and results from convex analysis can be found in §A.

Let $h\in\Gamma_{0}(\mathcal{H})$. Then, for every $x\in\mathcal{H}$, the function $z\mapsto\frac{1}{2}\left\|x-z\right\|^{2}+h(z)$ achieves its infimum at a unique point denoted by $\prox_{h}x$. The uniquely-valued operator $\prox_{h}:\mathcal{H}\to\mathcal{H}$ thus defined is the proximity operator or proximal mapping of $h$.

Let $h\in\Gamma_{0}(\mathcal{H})$, then for any $x\in\mathcal{H}$

where $v=\alpha D^{-1/2}u$ and $\alpha$ is the unique root of

It is of course straightforward to compute $\prox^{V}_{h^{*}}$ from $\prox^{V}_{h}$ either using Theorem 7, or using this theorem together with Corollary 6 and the Sherman-Morrison inversion lemma.

1.2 Diagonal+rank-1: Separable case

The following corollary is key to our novel optimization algorithm.

where $v=\alpha u$ and $\alpha$ is the unique root of

Let $k=\sum_{i=1}^{N}k_{i}$. Then $\prox^{V}_{h}(x)$ can be obtained exactly by sorting at most the $k$ real values $\left({\tfrac{d_{i}}{u_{i}}(x_{i}-t_{j})}\right)_{(i,j)\in\{1,\ldots,N\}\times\{1,\ldots,k_{i}\}}$.

Recall that (10) has a unique solution. When $\prox_{h_{i}}$ is piecewise affine with $k_{i}$ segments, it is easy to see that $p(\alpha)$ in (12) is also piecewise affine with slopes and intercepts changing at the $k$ transition points $\left({\tfrac{d_{i}}{u_{i}}(x_{i}-t_{j})}\right)_{(i,j)\in\{1,\ldots,N\}\times\{1,\ldots,k_{i}\}}$. To get $\alpha^{\star}$, it is sufficient to isolate the unique segment that intersects the abscissa axis. This can be achieved by sorting the values of the transition points which can cost in average complexity $O(k\log k)$. ∎

Corollary 9 can be extended to the “block” separable (i.e. separable in subsets of coordinates) when $D$ is piecewise constant along the same block indices.

1.3 Semi-smooth Newton method

In many situations (see examples below), the root of $p(\alpha)$ can be found exactly in polynomial complexity. If no closed-form is available, one can appeal to some efficient iterative method to solve (10) (or (12)). As $p$ is Lipschitz-continuous, hence so-called Newton (slantly) differentiable, semi-smooth Newton are good such solvers, with the proviso that one can design a simple slanting function which can be algorithmically exploited.

The semi-smooth Newton method for the solution of (10) can be stated as the iteration

where $g$ is a generalized derivative of $p$.

If $\prox_{h\circ D^{-1/2}}$ is Newton differentiable with generalized derivative $G$, then so is the mapping $p$ with a generalized derivative

This follows from linearity and the chain rule [23, Lemma 3.5]. The second statement follows strict increasing monotonicity of $p$ as established in Theorem 7. ∎

Thus, as $p$ is Newton differentiable with nonsingular generalized derivative whose inverse is also bounded, the general semi-smooth Newton convergence theorem implies that (13) converges super-linearly to the unique root of (10).

1.4 Examples

Many functions can be handled very efficiently using our results above. For instance, Table 1 summarizes a few of them where we can obtain either an exact answer by sorting when possible, or else by minimizing w.r.t. to a scalar variable (i.e. finding the unique root of (10)).

To put Proposition10 on a more concrete footing, we briefly cover the positivity constraint explicitly. Let $V=D+uu^{T}$ and $h(x)=\imath_{\{x:\,x\geqslant 0\}}$. We will calculate

Since we work with $V^{-1}$ and not $V$, we will not use $p(\alpha)$ but rather $\hat{p}(\alpha)$ which will be used in a similar way to $p$.

If $(y,\lambda)$ is a primal-dual solution to (14), the KKT conditions must be satisfied:

Define the scalar $\alpha=u^{T}\lambda$. The key observation is that if $\alpha$ is known, then the problem is solved since it becomes separable and the solution is

where $\left({x_{i}}\right)_{+}:=\max(0,x_{i})$. Let $\lambda_{i}^{(\alpha)}:=\left({-(x_{i}+\alpha u_{i})/d_{i}}\right)_{+}$, so we search for a value of $\alpha$ such that $\alpha=u^{T}\lambda^{(\alpha)}$, or in other words, a root of $\hat{p}(\alpha)=\alpha-u^{T}\lambda^{(\alpha)}$.

Define $\hat{\alpha}_{i}$ to be the sorted values of $(-x_{i}/u_{i})$, so we see that $\hat{p}$ is linear in the regions $[\hat{\alpha}_{i},\hat{\alpha}_{i+1}]$ and so it is trivial to check if $\hat{p}$ has a root in this region. Thus the problem is reduced to finding the correct region $i$, which can be done efficiently by a binary search over $\log_{2}(n)$ values of $i$ since $\hat{p}$ is monotonic. To see that $\hat{p}$ is monotonic, we write it as

where $\chi_{i}(\alpha)$ encodes the positivity constraint in the argument of $\left({\cdot}\right)_{+}$ and is thus either $0$ or $1$, hence the slope is always positive.

A primal rank 1 SR1 algorithm

Following the conventional quasi-Newton notation, we let $B$ denote an approximation to the Hessian of $f$ and $H$ denote an approximation to the inverse Hessian. All quasi-Newton methods update an approximation to the (inverse) Hessian that satisfies the secant condition:

Algorithm 1 follows the SR1 method , which uses a rank-1 update to the inverse Hessian approximation at every step. The SR1 method is perhaps less well-known than BFGS, but it has the crucial property that updates are rank-1, rather than rank-2, and it is described “[SR1] has now taken its place alongside the BFGS method as the pre-eminent updating formula.” .

We propose two important modifications to SR1. The first is to use limited-memory, as is commonly done with BFGS. In particular, we use zero-memory, which means that at every iteration, a new diagonal plus rank-one matrix is formed. The other modification is to extend the SR1 method to the general setting of minimizing $f+h$ where $f$ is smooth but $h$ need not be smooth; this further generalizes the case when $h$ is an indicator function of a convex set. Every step of the algorithm replaces $f$ with a quadratic approximation, and keeps $h$ unchanged. Because $h$ is left unchanged, the subgradient of $h$ is used in an implicit manner, in comparison to methods such as that use an approximation to $h$ as well and therefore take an explicit subgradient step.

In our experience, the choice of $H_{0}$ is best if scaled with a Barzilai-Borwein spectral step length

(we call it $\tau_{\text{BB}2}$ to distinguish it from the other Barzilai-Borwein step size $\tau_{\text{BB}1}=\left\langle s_{k},s_{k}\right\rangle/\left\langle s_{k},y_{k}\right\rangle\geqslant\tau_{\text{BB}2}$).

In SR1 methods, the quantity $\left\langle s_{k}-H_{0}y_{k},y_{k}\right\rangle$ must be positive in order to have a well-defined update for $u_{k}$. The update is:

A value of $\gamma=0.8$ works well in most situations. We have tested picking $\gamma$ adaptively, as well as trying $H_{0}$ to be non-constant on the diagonal, but found no consistent improvements.

Numerical experiments and comparisons

Consider the unconstrained LASSO problem (1). Many codes, such as and L-BFGS-B , handle only non-negativity or box-constraints. Using the standard change of variables by introducing the positive and negative parts of $x$, the LASSO can be recast as

Our second example uses a square operator $A$ with dimensions $n=13^{3}=2197$ chosen as a 3D discrete differential operator. This example stems from a numerical analysis problem to solve a discretized PDE as suggested by . For this example, we set $\lambda=1$. For all the solvers, we use the same parameters as in the previous example. Unlike the previous example, Fig. 11(b) now shows that L-BFGS-B is very slow on this problem. The FPC-AS method, very slow on the earlier test, is now the fastest. However, just as before, our SR1 method is nearly as good as the best algorithm. This robustness is one benefit of our approach, since the method does not rely on active-set identifying parameters and inner iteration tolerances.

Conclusions

In this paper, we proposed a novel variable metric (quasi-Newton) forward-backward splitting algorithm, designed to efficiently solve non-smooth convex problems structured as the sum of a smooth term and a non-smooth one. We introduced a class of weighted norms induced by a diagonal+rank 1 symmetric positive definite matrices, and proposed a whole framework to compute a proximity operator in the weighted norm. The latter result is distinctly new and is of independent interest. We also provided clear evidence that the non-diagonal term provides significant acceleration over diagonal matrices.

The proposed method can be extended in several ways. Although we focused on forward-backward splitting, our approach can be easily extended to the new generalized forward-backward algorithm of . However, if we switch to a primal-dual setting, which is desirable because it can handle more complicated objective functionals, updating $B_{k}$ is non-obvious. Though one can think of non-diagonal pre-conditioning methods.

Another improvement would be to derive efficient calculation for rank-2 proximity terms, thus allowing a 0-memory BFGS method. We are able to extend (result not presented here) Theorem 7 to diagonal+rank $r$ matrices. However, in general, one must solve an $r$-dimensional inner problem using the semismooth Newton method.

A final possible extension is to take $B_{k}$ to be diagonal plus rank-1 on diagonal blocks, since if $h$ is separable, this is still can be solved by our algorithm (see Remark 10). The challenge here is adapting this to a robust quasi-Newton update. For some matrices that are well-approximated by low-rank blocks, such as H-matrices , it may be possible to choose $B_{k}\equiv B$ to be a fixed preconditioner.

SB would like to acknowledge the Fondation Sciences Mathématiques de Paris for his fellowship.

Appendix A Elements from convex analysis

We here collect some results from convex analysis that are key for our proof. Some lemmata are listed without proof and can be either easily proved or found in standard references such as .

Let $\mathcal{C}$ a nonempty subset of $\mathcal{H}$. The indicator function $\imath_{\mathcal{C}}$ of $\mathcal{C}$ is

$\dom(\imath_{\mathcal{C}})=\mathcal{C}$.

$(h\circ A)^{*}=h^{*}\circ\left({A^{-1}}\right)^{*}$ if $A$ is a linear invertible operator.

$(h(x-x_{0}))^{*}(v)=h^{*}(v)+\left\langle v,x_{0}\right\rangle$.

Separability: $\left({\sum_{i=1}^{n}h_{i}(x_{i})}\right)^{*}(v_{1},\cdots,v_{n})=\sum_{i=1}^{n}h_{i}^{*}(v_{i})$, where $(x_{1},\cdots,x_{n})\in\mathcal{H}_{1}\times\cdots\times\mathcal{H}_{n}$.

Conjugate of a sum: assume $h_{1},h_{2}\in\Gamma_{0}(\mathcal{H})$ and the relative interiors of their domains have a nonempty intersection. Then

Let $Q$ be a symmetric positive semi-definite matrix. Let $Q^{+}$ be its Moore-Penrose pseudo-inverse. Then,

The subdifferential of a proper convex function $h\in\Gamma_{0}(\mathcal{H})$ at $x\in\mathcal{H}$ is the set-valued map $\partial h:\mathcal{H}\to 2^{\mathcal{H}}$

An element $v$ of $\partial h$ is called a subgradient.

The subdifferential map $\partial h$ is a maximal monotone operator from $\mathcal{H}\to 2^{\mathcal{H}}$.

If $h$ is (Gâteaux) differentiable at $x$, its only subgradient at $x$ is its gradient $\nabla h(x)$.

The duality formula to be stated shortly will be very useful throughout the rest of the paper.

with the relashionships between $x^{\star}$ and $u^{\star}$, respectively the solutions of the primal and dual problems

A.2 Proximal calculus in ℋ\mathcal{H}

The function $\mathchoice{\hphantom{{}^{{{\rho}}}}h^{{\kern-7.32623pt{\rho}\kern 4.68175pt}}_{{\kern-4.29286pt\kern 4.68175pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-7.32623pt{\rho}\kern 4.68175pt}}_{{\kern-4.29286pt\kern 4.68175pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-4.74384pt{\rho}\kern 2.82318pt}}_{{\kern-2.4343pt\kern 2.82318pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-3.93721pt{\rho}\kern 2.01656pt}}_{{\kern-1.62767pt\kern 2.01656pt}}}(x)=\inf_{z\in\mathcal{H}}\frac{1}{2\rho}\left\|x-z\right\|^{2}+h(z)$ for $0<\rho<+\infty$ is the Moreau envelope of index $\rho$ of $h$.

$\mathchoice{\hphantom{{}^{{{\rho}}}}h^{{\kern-7.32623pt{\rho}\kern 4.68175pt}}_{{\kern-4.29286pt\kern 4.68175pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-7.32623pt{\rho}\kern 4.68175pt}}_{{\kern-4.29286pt\kern 4.68175pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-4.74384pt{\rho}\kern 2.82318pt}}_{{\kern-2.4343pt\kern 2.82318pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-3.93721pt{\rho}\kern 2.01656pt}}_{{\kern-1.62767pt\kern 2.01656pt}}}$ is also the infimal convolution of $h$ with $\frac{1}{2\rho}\left\|\cdot\right\|^{2}$.

Translation: $\prox_{h(\cdot-y)}(x)=y+\prox_{h}(x-y)$.

Scaling: $\forall\rho\in(-\infty,\infty),\prox_{h(\rho\cdot)}(x)=\prox_{\rho^{2}f}(\rho x)/\rho$.

Let $h\in\Gamma_{0}(\mathcal{H})$. Then its Moreau envelope $\mathchoice{\hphantom{{}^{{{\rho}}}}h^{{\kern-7.32623pt{\rho}\kern 4.68175pt}}_{{\kern-4.29286pt\kern 4.68175pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-7.32623pt{\rho}\kern 4.68175pt}}_{{\kern-4.29286pt\kern 4.68175pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-4.74384pt{\rho}\kern 2.82318pt}}_{{\kern-2.4343pt\kern 2.82318pt}}}{\hphantom{{}^{{{\rho}}}}h^{{\kern-3.93721pt{\rho}\kern 2.01656pt}}_{{\kern-1.62767pt\kern 2.01656pt}}}$ is convex and Fréchet-differentiable with $1/\rho$-Lipschitz gradient

Let $h\in\Gamma_{0}(\mathcal{H})$, then for any $x\in\mathcal{H}$

Appendix B Proofs

B.2 Proof of Theorem 7

Let $p=\prox^{V}_{h}(x)$. Then, we have to solve

By virtue of Lemma 25, $\mathchoice{\hphantom{{}^{{{1}}}}{\left({h^{*}\circ D^{1/2}}\right)}^{{\kern-36.96486pt{1}\kern 34.40375pt}}_{{\kern-34.01486pt\kern 34.40375pt}}}{\hphantom{{}^{{{1}}}}{\left({h^{*}\circ D^{1/2}}\right)}^{{\kern-36.96486pt{1}\kern 34.40375pt}}_{{\kern-34.01486pt\kern 34.40375pt}}}{\hphantom{{}^{{{1}}}}{\left({h^{*}\circ D^{1/2}}\right)}^{{\kern-23.13829pt{1}\kern 21.27718pt}}_{{\kern-20.88829pt\kern 21.27718pt}}}{\hphantom{{}^{{{1}}}}{\left({h^{*}\circ D^{1/2}}\right)}^{{\kern-19.34482pt{1}\kern 17.4837pt}}_{{\kern-17.09482pt\kern 17.4837pt}}}$ is continuously differentiable with 1-Lipschitz gradient. Together with Lemma 22(24), 20 and 26, this yields

where $v^{\star}$ is the unique solution to the above dual problem (28). This problem amounts to minimizing a proper convex smooth continuously differentiable objective with a Lipschitz gradient over a linear set. The latter can be parametrized by a real scalar $\alpha$ such that $v=\alpha q=\alpha D^{-1/2}u$, and is then equivalent to solving the scalar strongly convex smooth optimization problem

whose solution $\alpha^{\star}$ is unique. This is equivalent to saying that $\alpha^{\star}$ is the unique root of

where we used again Lemma 25 and 26. Lipschitz continuity of $p(\alpha)$ follows from non-expansiveness of the proximal mapping, and the Lipschitz constant is straightforward from the triangle and Cauchy-Schwartz inequalities.

Let’s turn now to strict increasing monotonicity of $p$. Let $\beta>\alpha$. Denote the operator $P=\prox_{h\circ D^{-1/2}}$. Then,

where the first inequality is a consequence of the fact that the proximal mapping is firmly non-expansive. ∎

References