An Inexact Successive Quadratic Approximation Method for Convex L-1 Regularized Optimization

In this paper, we study an inexact Newton-like method for solving optimization problems of the form

where $f$ is a smooth convex function and $\mu>0$ is a (fixed) regularization parameter. The method constructs, at every iteration, a piecewise quadratic model of $\phi$ and minimizes this model inexactly to obtain a new estimate of the solution.

The piecewise quadratic model is defined, at an iterate $x_{k}$, as

where $g(x_{k})\stackrel{{\scriptstyle\rm def}}{{=}}\nabla f(x_{k})$ and $H_{k}$ denotes the Hessian $\nabla^{2}f(x_{k})$ or a quasi-Newton approximation to it. After computing an approximate solution $\hat{x}$ of this model, the algorithm performs a backtracking line search along the direction $d_{k}=\hat{x}-x_{k}$ to ensure decrease in the objective $\phi$.

We refer to this method as the successive quadratic approximation method in analogy to the successive quadratic programming method for nonlinear programming. This method is also known in the literature as a “proximal Newton method” [20, 24], but we prefer not to use the term “proximal” in this context since the quadratic term in $(\ref{quadm})$ is better interpreted as a second-order model rather than as a term that simply restricts the size of the step. The paper covers both the cases when the quadratic model $q_{k}$ is constructed with an exact Hessian or a quasi-Newton approximation.

The two crucial ingredients in the inexact successive quadratic approximation method are the algorithm used for the minimization of the model $q_{k}$, and the criterion that controls the degree of inexactness in this minimization. In the first part of the paper, we propose an inexactness criterion for the minimization of $q_{k}$ and prove that it guarantees global convergence of the iterates, and that it can be used to control the local rate of convergence. This criterion is based on the optimality conditions for the minimization of $(\ref{quadm})$, expressed in the form of a semi-smooth function that is derived from the soft-thresholding operator.

The second part of the paper is devoted to the practical implementation of the method. Here, the choice of algorithm for the inner minimization of the model $q_{k}$ is vital, and we consider two options: fista [4], which is a first-order method, and an orthant-based method [1, 7, 8]. The latter is a second-order method where each iteration consists of an orthant-face identification phase, followed by the minimization of a smooth model restricted to that orthant. The subspace minimization can be performed by computing a quasi-Newton step or a Newton-CG step (we explore both options). A projected bactracking line search is then applied; see section 5.3.

Some recent work on successive quadratic approximation methods for problem (1.1) include: Hsie et al. [12], where (1.2) is solved using a coordinate descent method, and which focuses on the inverse covariance selection method; [25] which also employs coordinate descent but uses a different working set identification than [12], and makes use of a quasi-Newton model; and Olsen et al. [18], where the inner solver is fista. None of these papers address convergence for inexact solutions of the subproblem. Recently Lee, Sun and Saunders [13] presented an inexact proximal Newton method that, at first glance, appears to be very close to the method presented here. Their inexactness criterion is, however, different from ours and suffers from a number of drawbacks, as discussed in section 2.

Inexact methods for solving generalized equations have been studied by Patricksson [21], and more recently by Dontchev and Rockafellar [10]. Special cases of the general methods described in those papers result in inexact sequential quadratic approximation algorithms. Patricksson [21] presents convergence analyses based on two conditions for controlling inexactness. The first is based on running the subproblem solver for a limited number of steps. The second rule requires that the residual norm be sufficiently small, but it does not cover the inexactness conditions presented in this paper (since the residual is computed differently and their inexactness measure is is different from ours). The rule suggested in Dontchev and Rockafellar [10] is very general, but it too does not cover the condition presented in this paper. Our rule, and those presented in [10, 13], is inspired by the classical inexactness condition proposed by Dembo et al. [9], and reduces to it for the smooth unconstrained minimization case (i.e. when $\mu=0$).

Another line of research that is relevant to this paper is the global and rate of convergence analysis for inexact proximal-gradient algorithms, which can be seen as special cases of sequential quadratic approximation without acceleration [15, 23, 26]. The inexactness conditions applied in those papers require that the subproblem objective function value be $\epsilon$-close to the optimal subproblem objective [15, 26], or that the approximate solution be exact with respect to an $\epsilon$-perturbed subdifferential [23], for a decreasing sequence $\{\epsilon\}$.

Our interest in the successive quadratic approximation method is motivated by the fact that it has not received sufficient attention from a practical perspective, where inexact solutions to the inner problem (1.2) are imperative. Although a number of studies have been devoted to the formulation and analysis of proximal Newton methods for convex composite optimization problems, as mentioned above, the viability of the approach in practice has not been fully explored.

This paper is organized in 5 sections. In section 2 we outline the algorithm, including the inexactness criteria that govern the solution of the subproblem $(\ref{quadm})$. In sections 3 and 4, we analyze the global and local convergence properties of the algorithm. Numerical experiments are reported in section 5. The paper concludes in section 6 with a summary of our findings, and a list of questions to explore.

Given an iterate $x_{k}$, an iteration of the algorithm begins by forming the model $(\ref{quadm})$, where $\mu>0$ is a given scalar and $H_{k}\succ 0$ is an approximation to the Hessian $\nabla^{2}f(x_{k})$. Next, the algorithm computes an approximate minimizer $\hat{x}$ of the subproblem

The point $\hat{x}$ defines the search direction $d_{k}=\hat{x}-x_{k}$. The algorithm then performs a backtracking line search along the direction $d_{k}$ that ensures sufficient decrease in the objective $\phi$. The minimization of (2.1) should be performed by a method that exploits the structure of this problem.

In order to compute an adequate approximate solution to $(\ref{quadm})$, we need some measure of closeness to optimality. In the case of smooth unconstrained optimization, (i.e. $(\ref{l1prob})$ with $\mu=0$), the norm of the gradient is a standard measure of optimality, and it is common [9] to require the approximate solution $\hat{x}$ to satisfy the condition

The term on the left side of $(\ref{smoothcond})$ is a measure of optimality for the model $q_{k}(x)$, in the unconstrained case.

For problem (1.1), the length of the iterative soft-thresholding (ISTA) step is a natural measure of optimality. The ISTA iteration is given by

where $\tau>0$ is a fixed parameter. It is easy to verify that $\|x_{\rm ista}-x_{k}\|$ is zero if and only if $x_{k}$ is a solution of problem (1.1). We need to express $\|x_{\rm ista}-x_{k}\|$ in a way that is convenient for our analysis, and for this purpose we note [16] that some algebraic manipulations show that $\|x_{\rm ista}-x_{k}\|=\tau\|F(x_{k})\|$, where

Here $P(x)_{[-\mu,\mu]}$ denotes the component-wise projection of $x$ onto the interval $[-\mu,\mu]$, and $\tau$ is a positive scalar.

One can directly verify that (2.4) is a valid optimality measure by noting that $F(x)=0$ is equivalent to the standard necessary optimality condition for $(\ref{l1prob})$:

For the objective $q_{k}$ of (2.1), this function takes the form

Using the measures $(\ref{Fdef})$ and $(\ref{Fqdef})$ in a manner similar to $(\ref{smoothcond})$, leads to the condition $\|F_{q}(x_{k};\hat{x})\|\leq\eta_{k}\|F(x_{k})\|$. However, depending on the method used to approximately solve $(\ref{lassop})$, this does not guarantee that $\hat{x}-x_{k}$ is a descent direction for $\phi$. To achieve this, we impose the additional condition that the quadratic model is decreased at $\hat{x}$.

The method is summarized in Algorithm 2.1.

For now, we simply assume that the sequence $\{\eta_{k}\}$ in (2.6) satisfies $\eta_{k}\in[0,1)$, but in section 4 we show that by choosing $\{\eta_{k}\}$ and the parameter $\tau$ appropriately, the algorithm achieves a fast rate of convergence. One may wonder whether the backtracking line search of Step 3 might hinder sparsity of the iterates. Our numerical experience indicates that this is not the case because, in our tests, Algorithm 2 almost always accepts the unit steplength ($\alpha=1$).

It is worth pointing out that Lee et al. [13] recently proposed and analyzed an inexactness criterion that is similar to the first inequality of $(\ref{inx})$. The main difference is that they use the subgradient of $q_{k}$ on the left side of the inequality, and both norms are scaled by $H_{k}^{-1}$. They claim similar convergence results to ours, but a worrying consequence of the lack of continuity of the subgradient of $q_{k}$ is that their inexactness condition can fail for vectors $x$ arbitrarily close to the exact minimizer of $q_{k}$. As a result, their criterion is not an appropriate termination test for the inner iteration. (In addition, their use of the scaling $H_{k}^{-1}$ precludes setting $H_{k}=\nabla^{2}f(x_{k})$, except for small or highly structured problems.)

In this section, we show that Algorithm 2.1 is globally convergent under certain assumptions on the function $f$ and the (approximate) Hessians $H_{k}$. Specifically, we assume that $f$ is a differentiable function with Lipschitz continuous gradient, i.e., there is a constant $M>0$ such that

for all $x,y$. We denote by $\lambda_{\min}(H_{k})$ and $\lambda_{\max}(H_{k})$ the smallest and largest eigenvalues of $H_{k}$, respectively.

Proof. We first show that if $\hat{x}$ is an approximate solution of $(\ref{quadm})$ that satisfies the inexactness conditions $(\ref{inx})$, then there is a constant $\gamma>0$ (independent of $k$) such that for all $k\in\{0,1,\cdots\}$

where $\ell_{k}$ and $F$ are defined in $(\ref{ell})$ and $(\ref{Fdef})$. To see this, note that by (2.6)

and therefore

Next, since $F(x_{k})=F_{q}(x_{k};x_{k})$, and using $(\ref{inx})$ and the contraction property of the projection, we have that

Combining this expression with $(\ref{lpre})$, we obtain $(\ref{ldec})$ for

Note that $\gamma>0$ as $\tau,\lambda,\Lambda>0$ and $\eta\in[0,1)$.

Let us define the search direction as $d=\hat{x}-x_{k}$. We now show that by performing a line search along $d$ we can ensure that the algorithm provides sufficient decrease in the objective function $\phi$, and this will allow us to establish the limit (3.2).

Since $g(x)$ satisfies the Lipschitz condition $(\ref{lips})$, we have

and thus

Recalling the definition of $\ell_{k}$, we have

By convexity of the $\ell_{1}$-norm, we have that

Combining this inequality with (3.4), and recalling that $x+d=\hat{x}$, we obtain for $\theta\in(0,1)$,

provided $\left((1-\theta)\lambda-M\alpha\right)\geq 0$. Therefore, the sufficient decrease condition (2.7) is satisfied for any steplength $\alpha$ satisfying

and if the backtracking line search cuts the steplength in half (say) after each trial, we have that the steplength chosen by the line search satisfies

Thus, from (3.5) and (3.3) we obtain

Since $f$ is assumed to be bounded below, so is the objective function $\phi$, and given that the decrease in $\phi$ is proportional to $\|F(x_{k})\|$ we obtain the limit $(\ref{limit})$. $\Box$

We note that to establish this convergence result it was not necessary to assume convexity of $f$.

To analyze the local convergence rate of the successive quadratic approximation method, we use the theory developed in Chapter 7 of Facchinei and Pang [11]. To do this, we first show that, if $x^{*}$ is a nonsingular minimizer of $\phi$, then the functions $F_{q}(x;\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ are a family of uniformly Lipschitzian nonsingular homeomorphisms for all $x$ in a neighborhood of $x^{*}$.

Proof. It is straightforward to show that for any scalars $a\neq b$ and interval $[-\mu,\mu]$,

Therefore for any vectors $y$ and $z$, and for any index $i\in\{1,\cdots,n\}$ we have

where $\bar{d}_{i}\in[0,1]$ is a scalar implied by (4.2). This implies that

where $D={\rm diag}(\bar{d}_{i})$. Hence

Since the right hand side is a quadratic form, we symmetrize the matrix, and if we let $w=z-y$, the right side is

To show that the symmetric matrix inside the square brackets is positive definite, we note that since $(\tau H-D)^{T}(\tau H-D)=\tau^{2}H^{2}-\tau(HD+DH)+D^{2}$ is positive semi-definite, we have that

Substituting this into (4.4) yields

since $D-D^{2}/2$ is positive semi-definite given that the elements of the diagonal matrix $D$ are in $[0,1]$. If $\lambda_{i}$ is an eigenvalue of $H$, the corresponding eigenvalue of the matrix $H-{\tau\over 2}H^{2}$ is $\lambda_{i}-\tau\lambda_{i}^{2}/2\geq\lambda_{i}/2$ since our assumption on $\tau$ implies $1>\tau\|H\|\geq\tau\lambda_{i}$. Therefore, we have from (4.3) that

$\Box$

Inequality (4.1) establishes that $F_{q}(x;\cdot)$ is strongly monotone. Next we show that, when $H$ is defined as the Hessian of $f$, the functions $F_{q}(x;\cdot)$ are homeomorphisms and that they represent an accurate approximation to the function $F$ defined in (2.4).

Proof. Since $\nabla^{2}f(x)$ is continuous, there is a neighborhood ${\cal N}$ of $x^{*}$ and a positive constant $\lambda$ such that $\lambda_{min}(\nabla^{2}f(x))\geq\lambda>0$ and $\tau\|\nabla^{2}f(x)\|<1$, for all $x\in\cal{N}$. It follows from Lemma 4.1 that for any such $x$, the function $F_{q}(x;y)$ given by (4.5) is strongly (or uniformly) monotone with constant greater than $\lambda/2$. We now invoke the Uniform Monotonicity Theorem (see e.g. Theorem 6.4.4 in [19]), which states that if a function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is continuous and uniformly monotone, then $F$ is a homeomorphism of $\mathbb{R}^{n}$ onto itself. We therefore conclude that $F_{q}(x;y)$ is a homeomorphism.

In addition, we have from $(\ref{smono})$ and the Cauchy-Schwartz inequality that

which implies Lipschitz continuity of $F_{q}^{-1}(x;\cdot)$ with constant $2/\lambda$. To establish $(\ref{fqerror})$, note that

and thus

by the non-expansiveness of a projection onto a convex set and Taylor’s theorem. $\Box$

Theorem 4.2 shows that $F_{q}(x;y)$ defines a strong nonsingular Newton approximation in the sense of Definition 7.2.2 of Pang and Facchinei [11]. This implies quadratic convergence for the (exact) successive quadratic approximation (SQA) method.

Proof. By Theorem 4.2, $F_{q}(x_{k};y)$ satisfies the definition of a nonsingular strong Newton approximation of $F$ at $x^{*}$, given by Facchinei and Pang ([11], 7.2.2) and thus by Theorem 7.2.5 of that book the local convergence is quadratic. $\Box$

Now we consider the inexact SQA algorithm that, at each step, computes a point $y$ satisfying

where $r_{k}$ is a vector such that $\|r_{k}\|\leq\eta_{k}\|F(x_{k})\|$ with $\eta_{k}<1$; see (2.6). We obtain the following result for a method that sets $x_{k+1}=y$.

Proof. By Theorem 4.2, the iteration described in the statement of the theorem satisfies all the conditions of Theorem 7.2.8 of [11]. The results then follow immediately from that theorem. $\Box$

We have shown above the the inexact successive quadratic approximation (SQA) method with $\alpha_{k}=1$ yields a fast rate of convergence. We now show that this inexact SQA algorithm will select the steplength $\alpha_{k}=1$ in a neighborhood of the solution. In order to do so, we strengthen the inexactness conditions (2.6) slightly so that they read

where $\eta_{k}<1$, $\zeta\in(\theta,1/2)$ and $\theta$ is the input parameter of Algorithm 2.1 used in (2.7). Thus, instead of simple decrease, we now impose sufficient decrease in $q_{k}$.

Proof. If we denote by $\bar{y}$ the (exact) minimizer of $q_{k}$, we claim that

To see this, note that $q_{k}(y)=\ell_{k}(y)+{1\over 2}(y-x_{k})^{T}H_{k}(y-x_{k})$, and since $\ell_{k}$ and $q_{k}$ are convex and $\bar{y}$ minimizes $q_{k}$, there exists a vector $v\in\partial\ell_{k}(\bar{y})$ such that $v+H_{k}(\bar{y}-x_{k})=0$. By convexity

Therefore,

which proves (4.9).

Now consider the continuous function $q_{k}(x)-q_{k}(x_{k})-\zeta(\ell_{k}(x)-\ell_{k}(x_{k}))$. By $(\ref{masia})$ its value at $x=\bar{y}$ is

where the last inequality follows from (4.10). Therefore by continuity, the value of this function for any $x$ in some neighborhood of $\bar{y}$ is negative, implying that $(\ref{inxx})$ is satisfied by any approximate solution $\hat{x}$ sufficiently close to $\bar{y}$. $\Box$