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

Introduction

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” , 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 , which is a first-order method, and an orthant-based method . 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. , where (1.2) is solved using a coordinate descent method, and which focuses on the inverse covariance selection method; which also employs coordinate descent but uses a different working set identification than , and makes use of a quasi-Newton model; and Olsen et al. , where the inner solver is fista. None of these papers address convergence for inexact solutions of the subproblem. Recently Lee, Sun and Saunders 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 , and more recently by Dontchev and Rockafellar . Special cases of the general methods described in those papers result in inexact sequential quadratic approximation algorithms. Patricksson 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 is very general, but it too does not cover the condition presented in this paper. Our rule, and those presented in , is inspired by the classical inexactness condition proposed by Dembo et al. , 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 . The inexactness conditions applied in those papers require that the subproblem objective function value be $\epsilon$-close to the optimal subproblem objective , or that the approximate solution be exact with respect to an $\epsilon$-perturbed subdifferential , 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.

Notation. In the remainder, we let $g(x_{k})=\nabla f(x_{k})$, and let $\|\cdot\|$ denote any vector norm. We sometimes abbreviate successive quadratic approximation method as “SQA method”, and note that this algorithm is often referred to in the literature as the “proximal Newton method”.

The Algorithm

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 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 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}$.

Inexactness Conditions. A point $\hat{x}$ is considered an acceptable approximate solution of subproblem (2.1) if

for some parameter $0\leq\eta_{k}<1$, where $\|\cdot\|$ is any norm. (Note that $F_{q}(x_{k};x_{k})=F(x_{k})$, so that the first condition can also be written as $\|F_{q}(x_{k};\hat{x})\|\leq\eta_{k}\|F(x_{k})\|$.)

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. 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.)

Global Convergence

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.

Suppose that $f$ is a smooth function that is bounded below and that satisfies $(\ref{lips})$. Let $\{x_{k}\}$ be the sequence of iterates generated by Algorithm 2.1, and suppose that there exist constants $0<\lambda\leq\Lambda$ such that the sequence $\{H_{k}\}$ satisfies

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\}$

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

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

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$.

Local Convergence

If $H$ is a symmetric positive definite matrix with smallest eigenvalue $\lambda>0$, then the function of $y$ given by

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$ is a scalar implied by (4.2). This implies that

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

since $D-D^{2}/2$ is positive semi-definite given that the elements of the diagonal matrix $D$ are in $$. 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

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).

If $\nabla^{2}f(x^{*})$ is positive definite and $\tau<1/\|\nabla^{2}f(x^{*})\|$, then there is a neighborhood $\cal{N}$ of $x^{*}$ such that for all $x\in\cal{N}$ the functions of $y$ given by

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

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 . This implies quadratic convergence for the (exact) successive quadratic approximation (SQA) method.

If $\nabla^{2}f(x)$ is Lipschitz continuous and positive definite at $x^{*}$, and $\tau<1/\|\nabla^{2}f(x^{*})\|$, then there is a neighborhood of $x^{*}$ such that, if $x_{0}$ lies in that neighborhood, the iteration that defines $x_{k+1}$ as the unique solution to

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 (, 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$.

Suppose that $\nabla^{2}f(x)$ is Lipschitz continuous and positive definite at $x^{*}$, $\tau<1/\|\nabla^{2}f(x^{*})\|$, and that $x_{k+1}$ is computed by solving

Proof. By Theorem 4.2, the iteration described in the statement of the theorem satisfies all the conditions of Theorem 7.2.8 of . 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}$.

If $H_{k}$ is positive definite, the inexactness condition $(\ref{inxx})$ is satisfied by any sufficiently accurate solution to $(\ref{lassop})$.

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

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$

Suppose that $H_{k}=\nabla^{2}f(x_{k})$ in Algorithm 2.1, and that we modify Step 2 in that algorithm to require that the approximate solution $\hat{x}$ satisfies (4.8) instead of (2.6). If we assume that $\nabla^{2}f(x)$ is Lipschitz continuous, then for all $k$ sufficiently large we have $\alpha_{k}=1$.

Proof. Given that $\hat{x}=x_{k}+d_{k}$ satisfies $(\ref{inxx})$, it follows from Taylor’s theorem, the Lipschitz continuity of $\nabla^{2}f(x)$, and equation $(\ref{lpre})$ that for some constant $\rho>0$

if $\|d_{k}\|\leq(\zeta-\theta)\lambda/2\rho$ . Since the global convergence analysis implies $\|d_{k}\|\rightarrow 0$, we have from (2.7) that eventually the steplength $\alpha_{k}=1$ is accepted and used. $\Box$

We note that (4.8) is stronger than (2.6), and therefore, all the results presented in this and the previous section apply also to the strengthened condition (4.8). Theorem 4.6 implies that if Algorithm 2.1 is run with the strengthened accuracy condition $(\ref{inxx})$, and $H_{k}=\nabla^{2}f(x_{k})$, then once the iterates are close enough to a nonsingular minimizer $x^{*}$, the iterates have the linear, superlinear or quadratic convergence rates described in Theorem 4.4 if $\eta_{k}$ is chosen appropriately.

Numerical Results

To study this question, we explore various algorithmic options within the successive quadratic approximation method, and evaluate their performance using data sets with different characteristics. One of the data sets concerns the covariance selection problem (where the unknown is a matrix), and the other involves a logistic objective function (where the unknown is a vector). Our benchmark is fista applied directly to problem (1.1). fista enjoys convergence guarantees when applied to problem (1.1), and is generally regarded as an effective method.

The methods employed in our numerical tests are as follows.

FISTA. This is the fista algorithm applied to the original problem (1.1). We used the implementation from the TFOCS package, called N83 . This implementation differs from the (adaptive) algorithm described by Beck and Teboulle in the way the Lipschitz parameter is updated, and performed significantly better in our test set than the method in .

PNOPT. This is the sequential quadratic approximation (proximal Newton) method of Lee, Sun and Saunders . The Hessian $H_{k}$ in the subproblem (2.1) is updated using the limited memory BFGS formula, with a (default) memory of 50. (The pnopt package also allows for the use of the exact Hessian, but since this matrix must be formed and factored at each iteration, its use is impractical.) The subproblem (2.1) is solved using the N83 implementation of fista mentioned above. pnopt provides the option of using sparsa instead of N83 as an inner solver, but the performance of sparsa was not robust in our tests, and we will not report results with it.

SQA. Is the sequential quadratic approximation method described in Algorithm 2. We implemented 3 variants that differ in the method used to solve the subproblem (2.1).

SQA-FISTA. This is an sqa method using fista-n83 to solve the subproblem (2.1). The matrix $H_{k}$ is the exact Hessian $\nabla^{2}f(x_{k})$; each inner fista iteration requires two multiplications with $H_{k}$.

SQA-OBM-CG. This is an sqa method that employs an orthant based method to solve the subproblem (2.1). The obm method performs the subspace minimization step using a Newton-CG iteration. The number of CG iterations varies during the course of the (outer) iteration according to the rule $\min\{3,1+\lfloor k/10\rfloor\}$, where $k$ is the outer iteration number.

SQA-OBM-QN. This is an sqa method where the inner solver is an obm method in which the subspace phase consists of a limited memory BFGS step, with a memory of 50. The correction pairs used to update the quasi-Newton matrix employ gradient differences from the outer iteration (as in pnopt).

The initial point was set to the zero vector in all experiments, and the iteration was terminated if $\|F(x_{k})\|_{\infty}\leq 10^{-5}$, where $F$ is defined in (2.4). The maximum number of outer iterations for all solvers was $3000$. In the SQA method, the parameter $\eta_{k}$ in the inexactness condition (2.6) was defined as $\eta_{k}=\max\{1/k,0.1\}$, and we set $\theta=0.1$ in (2.7). For pnopt we set ‘ftol’=$1e-16$, and ‘xtol’=$1e-16$ (so that those two tests do not terminate the iteration prematurely), and chose ‘Lbfgs_mem’=50.

We noted above that Algorithm 2 can employ the inexactness conditions (2.6) or (4.8). We implemented both conditions, with $\zeta=\theta=0.1$, and obtained identical results in all our runs.

We now describe the numerical tests performed with these methods.

The task of estimating a high dimensional sparse inverse covariance matrix is closely tied to the topic of Gaussian Markov random fields , and arises in a variety of recognition tasks. This model can be used to recover a sparse social or genetic network from user or experimental data.

where $S$ is a given sample covariance matrix, $P$ denotes the unknown inverse covariance matrix, $\mu$ is the regularization parameter, and $\|P\|\stackrel{{\scriptstyle\rm def}}{{=}}\|vec(P)\|_{1}$. We note that the Hessian of the first two terms in (5.1) has a very special structure: it is given by $P^{-1}\otimes P^{-1}$.

Since the objective is not defined when $\det(P)\leq 0$, we define it as $+\infty$ in that case to ensure that all iterates remain positive definite. Such a strategy could, however, be detrimental to a solver like fista, and to avoid this we selected the starting point so that the condition $\det(P)\leq 0$ did not occur.

We employ three data sets: the well-known Estrogen and Leukemia test sets , and the problem given in Olsen et al. , which we call OONR. The characteristics of the data sets are given in Table 1, where nnz($P^{\ast}_{\mu}$) denotes the number of nonzeros in the solution.

The performance of the algorithms on these three test problems is given in Tables 2, 3 and 4. We note that fista does not perform inner iterations since it is applied directly to the original problem (1.1), and that pnopt-fista does not compute Hessian vector products because the matrix $H_{k}$ in the model (2.1) is defined by quasi-Newton updating. Each inner iteration of sl-obm-qn performs a Hessian-vector multiplication to compute the subproblem objective, and a multiplication of the inverse Hessian times a vector to compute the unconstrained minimizer on the active orthant face — we report these as two Hessian-vector products in Tables 2, 3 and 4.

Table 2. ESTROGEN; $\mu=0.5$, optimality tolerance = $10^{-5}$

Table 3. LEUKEMIA; $\mu=0.5$, optimality tolerance = $10^{-5}$

Table 4. OONR; $\mu=0.5$, optimality tolerance = $10^{-5}$

We now comment on the results given in Tables 2-4. For the inverse covariance selection problem (5.1), Hessian-vector products are not as expensive as for other problems (c.f. Tables 5-6) — in fact, these products are not much costlier than computations with the limited memory BFGS matrix. This fact, combined with the effectiveness of the obm method, makes sqa-obm-cg the most efficient of the methods tested. obm is a good subproblem solver due to its ability to estimate the set of zero variables quickly, so that the subspace step is computed in a small reduced space (the density of $P_{\mu}^{\ast}$ is less than $2.5\%$ for the three test problems.) In addition, the obm-cg method can decrease $\|F_{q}\|$ drastically in a single iteration, often yielding a high quality sqa step and thus a low number of outer iterations.

We note that the quasi-Newton algorithms sl-obm-qn and pnopt are different methods because of the subproblem solvers they employ. sl-obm-qn uses the two-phase obm method in which the quasi-Newton step is computed in a subspace, whereas pnopt applies the fista iteration to subproblem (1.2) where $H_{k}$ is a quasi-Newton matrix. Although the number of outer iterations of both methods is comparable for problems Estrogen and OONR, there is a large difference in the number of inner iterations due to power of the obm approach.

Note that the number of inner fista iterations in sqa-fista is always smaller than for fista. We repeated the experiment with problem OONR using looser optimality tolerances (TOL); the total number of fista is given in Table 5.

Table 5. Effect of convergence tolerance TOL; OONR

These results are typical for the covariance selection problems, where the sqa-fista is clearly more efficient than fista; we will see that this is not the case for the problems considered next.

2 Logistic Regression Problems

We employed the data given in Table 6, which was downloaded from the SVMLib repository. The values of the regularization parameter $\mu$ were taken from Lee et al. .

Table 6. Test problems for logistic regression tests Data set $N$ number of features $\mu$ nnz($x^{\ast}_{\mu}$) Gisette (scaled) 6,000 5,000 6.67e-04 482 (9.64%) RCV1 (binary) 20,242 47,236 3.31e-04 140 (0.30%)

Table 7. GISETTE; $\mu=6.67e-04$, optimality tolerance = $10^{-5}$

Table 8. RCV1; $\mu=3.31e-04$, optimality tolerance = $10^{-5}$

For the logistic regression problems, Hessian-vector products are expensive, particularly for gisette, where the data set is dense. As a result, the obm variant that employs quasi-Newton approximations, namely sqa-obm-qn, performs best (even though sqa-obm-cg requires a smaller number of outer iterations). Note that sqa-fista is not efficient; in fact it requires a much larger number of inner iterations than the total number of iterations in fista. In Table 9 we observe the effect of the optimality tolerance, on these two methods, using problem gisette.

Table 9. Effect of convergence tolerance TOL ; Gisette

We observe from Table 9, that fista requires a smaller number of iterations; it is only for a very high accuracy of $10^{-6}$ that sqa-fista becomes competitive. This is in stark contrast with Table 5.

In summary, for the logistic regression problems the advantage of the sqa method is less pronounced than for the inverse covariance estimation problems, and is achieved only through the appropriate choice of model Hessian $H_{k}$ (quasi-Newton) and the appropriate choice of inner solver (active set obm method).

3 Description of the orthant based method (OBM)

We conclude this section by describing the orthant-based method used in our experiments to solve the subproblem (2.1). We let $t$ denote the iteration counter of the obm method, and let $z_{t}$ denote its iterates.

where $v_{t}$ is the minimum norm subgradient of $q_{k}$ computed at $z_{t}$, i.e.,

Defining $\Omega_{t}$ in this manner was proposed, among others, by Andrew and Gao . In the relative interior of $\Omega_{t}$, the model function $q_{k}$ is differentiable. The active set in the orthant-based method, defined as $A^{k}=\{i:\omega^{k}_{i}=0\}$, determines the variables that are kept at zero, while the rest of the variables are chosen to minimize a (smooth) quadratic model. Specifically, the search direction $d_{t}$ of the algorithm is given by $d_{t}=\hat{z}-z_{t}$, where $\hat{z}$ is a solution of

Note that $\psi(z)=f(x_{k})+(g(x_{k})+\omega_{t}\mu)^{T}(z-x_{k})+\frac{1}{2}(z-x_{k})^{T}H_{k}(z-x_{k})$.

In the obm-cg variant, we set $H_{k}=\nabla^{2}f(x_{k})$, and perform an approximate minimization of this problem using the projected conjugate gradient iteration . In the obm-qn version, $H_{k}$ is a limited memory BFGS matrix and $\hat{z}$ is the exact solution of (5.5). This requires computation of the inverse reduced Hessian $R_{k}=(Z_{k}^{T}\nabla^{2}H_{k}Z_{k})^{-1}$, where $Z_{k}$ is a basis for the space defined by (5.5). The matrix $R_{k}$ can be updated using the compact representations of quasi-Newton matrices . After the direction $d_{t}=\hat{z}-z_{t}$ has been computed, the obm method performs a line search along $d_{t}$, projecting the iterate back onto the orthant face $\Omega_{t}$, until a sufficient reduction in the function $q_{k}$ has been obtained. Although this algorithm performed reliably in our tests, its convergence has not been proved (to the best of our knowledge) because the orthant face identification procedure (5.2)-(5.5) can lead to arbitrarily small steps.

Our obm-qn algorithm differs from the owl method in two respects: it does not realign the direction $z-z_{t}$ so that the sign of its components match those of $v_{t}$, and it performs the minimization of the model exactly, while the owl method computes only an approximate solution – defined by computing the reduced inverse Hessian $Z_{k}^{T}\nabla^{2}H_{k}^{-1}Z_{k}$, instead of the inverse of the reduced Hessian $R_{k}$.

Final Remarks

One of the key ingredients in making the successive quadratic approximation (or proximal Newton) method practical for problem (1.1) is the ability to terminate the inner iteration as soon as a step of sufficiently good quality is computed. In this paper, we have proposed such an inexactness criterion; it employs an optimality measure that is tailored to the structure of the problem. We have shown that the resulting algorithm is globally convergent, that its rate of convergence can be controlled through an inexactness parameter, and that the inexact method will naturally accept unit step lengths in a neighborhood of the solution. We have also argued that our inexactness criterion is preferable to the one proposed by Lee et al. .

The method presented in this paper can use any algorithm for the inner minimization of the subproblem (1.2). In particular, all the results are applicable to the case when this inner minimization is performed using a coordinate descent algorithm . In our numerical tests we employed fista and an orthant-based method as the inner solvers, and found the latter method to be particularly effective. The efficacy of the successive quadratic approximation approach depends of the choice of matrix $H_{k}$ in (1.2), which is problem dependent: when Hessian-vector products are expensive to compute, then a quasi-Newton approximation is most efficient; otherwise defining $H_{k}$ as the exact Hessian and implementing a Newton-CG iteration is likely to give the best results.

Acknowledgement. The authors thank Jong-Shi Pang for his insights and advice throughout the years. The theory presented by Facchinei and Pang in the context of variational inequalities was used in our analysis, showing the power and generality that masterful book.

References