Practical Inexact Proximal Quasi-Newton Method with Global Complexity Analysis

Introduction

In this paper, we are interested in the following popular convex optimization problem:

Here $\|y\|^{2}_{H}$ denotes $y^{\top}Hy$. Clearly, the computational cost of approximately solving (1.2) depends on the choice of matrix $H$ and the solution approach.

We are particularly interested in the case of sparse optimization, where $g(x)=\lambda\|x\|_{1}$. While the theory we present here applies to the general form (1.1), the efficient method for solving (1.2) that we consider in this paper is designed with $g(x)=\lambda\|x\|_{1}$ example in mind. In this case problem (1.2) takes a form of an unconstrained Lasso problem . We consider matrices $H$ which are a sum of a diagonal matrix and a low-rank matrix and we apply randomized coordinate descent to solve (1.2) approximately. An extension to the group sparsity term $g(x)=\lambda\sum\|x_{i}\|_{2}$ , is rather straightforward.

Problems of the form (1.1) with $g(x)=\lambda\|x\|_{1}$ have been the focus of much research lately in the fields of signal processing and machine learning. This form encompasses a variety of machine learning models, in which feature selection is desirable, such as sparse logistic regression , sparse inverse covariance selection and unconstrained Lasso , etc. These settings often present common difficulties to optimization algorithms due to their large scale. During the past decade most optimization effort aimed at these problems focused on development of efficient first-order methods, such as accelerated proximal gradient methods , block coordinate descent methods and alternating directions methods . These methods enjoy low per-iteration complexity, but typically have slow local convergence rates. Their performance is often hampered by small step sizes. This, of course, has been known about first-oder methods for a long time, however, due to the very large size of these problems, second-order methods are often not a practical alternative. In particular, constructing and storing a Hessian matrix, let alone inverting it, is prohibitively expensive for values of $n$ larger than $10000$, which often makes the use of the Hessian in large-scale problems prohibitive, regardless of the benefits of fast local convergence rate.

Nevertheless, recently several new methods were proposed for sparse optimization which make careful use of second-order information . These new methods are designed to exploit the special structure of the Hessian of specific functions to improve efficiency of solving (1.2). Several successful methods employ coordinate descent to approximately solve (1.2). While other approaches to solve Lasso subproblem were considered in , none generally outperform coordinate descent, which is well suited when special structure of the Hessian approximation, $H$, can be exploited and when low accuracy of the subproblem solutions is sufficient. In particular, proposes a specialized GLMNET implementation for sparse logistic regression, where coordinate descent method is applied to the unconstrained Lasso subproblem constructed using the Hessian of $f(x)$. The special structure of the Hessian is used to reduce the complexity cost of each coordinate step so that it is linear in the number of training instances, and a two-level shrinking scheme proposed to focus the minimization on smaller subproblems. Similar ideas are used in in a specialized algorithm called QUIC for sparse inverse covariance selection, where the Hessian of $f(x)$ also has a favorable structure for solving Lasso subproblems. Another specialized method for graphical Markov random fields was recently proposed in . This method also exploits special Hessian structure to improve coordinate descent efficiency.

There are other common features shared by the methods described above. These methods are often referred to as proximal Newton-type methods. The overall algorithmic framework can be described as follows:

At each iteration $k$ the smooth function $f(x)$ is approximated near the current iterate $x^{k}$ by a convex quadratic function $q^{k}(x)$.

A working subset of coordinates (elements) of $x$ is selected for subproblem optimization.

Then $l(k)$ passes of coordinate descent are applied to optimize (approximately) the function $q^{k}(x)+g(x)$ over the working set, which results in a trial point. Here $l(k)$ is some linear function of $k$.

The trial point is accepted as the new iterate if it satisfies some sufficient decrease condition (to be specified).

Otherwise, a line search is applied to compute a new trial point.

In this paper we do not include the theoretical analysis of various working set selection strategies. Some of these have been analyzed in the prior literature (e.g., see ). Combining such existing analysis with the rate of convergence results in this paper is a subject of a future study.

This paper contains the following three main results.

We discuss theoretical properties of the above framework in terms of global convergence rates. In particular, we show that if we replace the line search by a prox-parameter update mechanism, we can derive sublinear global convergence results for the above methods under mild assumptions on Hessian approximation matrices, which can include diagonal, quasi-Newton and limited memory quasi-Newton approximations. We also provide the convergence rate for the case of inexact subproblem optimization. It turns out that standard global convergence analysis of proximal gradient methods (see ) does not extend in a natural way to proximal quasi-Newton frameworks, hence we use a different technique derived for smooth optimization in , in a novel way, to obtain the global complexity result.

The heuristic of applying $l(k)$ passes of coordinate descent to the subproblem is very useful in practice, but has not yet been theoretically justified, due to the lack of known complexity estimates. Here we use probabilistic complexity bounds of randomized coordinate descent to show that this heuristic is indeed well justified theoretically. In particular, it guarantees the sufficiently rapid decrease of the expectation of the error in the subproblems and hence allows for sublinear global convergence rate to hold for the entire algorithm (again, in expectation). This gives us the first complete global convergence rate result for the algorithmic schemes for practical (inexact) proximal Newton-type methods. Moreover, using the new analysis from we are able to provide lower overall complexity bound than the one that follows from

Finally, we propose an efficient general purpose algorithm that uses the same theoretical framework, but which does not rely on the special structure of the Hessian, and yet in our tests compares favorably with the state-of-the-art, specialized methods such as QUIC and GLMNET. We replace the exact Hessian computation by the limited memory BFGS Hessian approximations (LBFGS) and exploit their special structure within a coordinate descent approach to solve the subproblems.

Let us elaborate a bit further on the new approaches and results developed in this paper and discuss related prior work.

In Byrd et al. propose that the methods in the framework described above should be referred to as sequential quadratic approximation (SQA) instead of proximal Newton methods. They reason that there is no proximal operator or proximal term involved in this framework. This is indeed the case, if a line search is used to ensure sufficient decrease. Here we propose to consider a prox term as a part of the quadratic approximation. Instead of a line search procedure, we update the prox term of our quadratic model, which allows us to extend global convergence bounds of proximal gradient methods to the case of proximal (quasi-)Newton methods. The criteria for accepting a new iteration is based on sufficient decrease condition (much like in trust region methods, and unlike that in proximal gradient methods). We show that our mechanism of updating the prox parameter, based on sufficient decrease condition, leads to an improvement in performance and robustness of the algorithm compared to the line search approach as well as enabling us to develop global convergence rates.

Convergence results for the proximal Newton method have been shown in and more recently in (with the same sufficient decrease condition as ours, but applied within a line search). These papers also demonstrate super linear local convergence rate of the proximal Newton and a proximal quasi-Newton method. Thus it is confirmed in that using second order information is as beneficial for problems of the form (1.1) as it is for the smooth optimization problems. These results apply to our framework when exact Hessian (or a quasi-Newton approximation) of $f(x)$ is used to construct $q(x)$ (and if the matrices have bounded eigenvalues). However, theory in does not provide global convergence rates for these methods, and just as in the case on smooth optimization, the super linear local convergence rates, generally, do not apply in the case of LBFGS Hessian approximations.

The convergence rate that we show is sub linear, which is generally the best that can be expected from a proximal (quasi-)Newton method with no assumptions on the accuracy of the Hessian approximations. Practical benefits of using LBFGS Hessian approximations is well known for smooth optimization and have been exploited in many large scale applications. In this paper we demonstrate this benefit in the composite optimization setting (1.1). Some prior work showing benefit of limiter memory quasi-Newton method in proximal setting include . We also emphasize in our theoretical analysis the potential gain over proximal gradient methods in terms of constants occurring in the convergence rate.

To prove the sub linear rate we borrow a technique from . The technique used in and for the proof of convergence rates of the (inexact) proximal gradient method do not seem to extend to general positive definite Hessian approximation matrix. In a related work the authors analyze global convergence rates of an accelerated proximal quasi-Newton method, as an extension of FISTA method . The convergence rate they obtain match that of accelerated proximal gradient methods, hence it is a faster rate than that of our method presented here. However, they have to impose much stricter conditions on the Hessian approximation matrix, in particular they require that the difference between any two consecutive Hessian approximations (i.e., $H_{k}-H_{k+1}$) is positive semidefinite. This is actually in contradiction to FISTA’s requirement that the prox parameter is never increased. Such a condition is very restrictive as is impractical. In this paper we briefly show how results in and can be extended under some (also possibly strong) assumptions of the Hessian approximations, to give a simple and natural convergence rate analysis. We then present an alternative analysis, which only requires the Hessian approximations to have bounded eigenvalues. Investigating accelerated version of our approach without restrictive assumptions on the Hessian approximations and with the use of randomized coordinate descent is a subject of future research.

Finally, we use the complexity analysis of randomized coordinate descent in to provide a simple and efficient stopping criterion for the subproblems and thus derive the total complexity of proximal (quasi-)Newton methods based on randomized coordinate descent to solve Lasso subproblems.

The paper is organized as follows: in Section 2 we describe the algorithmic framework. Then, in Section 3 we present some of the assumptions and discussions and, in Section 3.1, convergence rate analysis based on and . In Section 4 we show the new convergence rate analysis using for exact and inexact version of our framework. We then extend the analysis for the cases where inexact solution to a subproblem is random in Section 5 and in particular to the randomized coordinate descent in Section 5.1. Brief description of the details of our proposed algorithm are in Section 6 and computational results validating the theory are presented in Section 7.

Basic algorithmic framework and theoretical analysis

The following function is used throughout the paper as an approximation of the objective function $F(x)$.

For a fixed point $\bar{x}$, the function $Q(H,x,\bar{x})$ serves as an approximation of $F(x)$ around $\bar{x}$. Matrix $H$ controls the quality of this approximation. In particular, if $f(x)$ is smooth and $H=\frac{1}{\mu}I$, then $Q(H,x,\bar{x})$ is a sum of the prox-gradient approximation of $f(x)$ at $\bar{x}$ and $g(x)$. This particular form of $H$ plays a key role in the design and analysis of proximal gradient methods (e.g., see ) and alternating direction augmented Lagrangian methods (e.g, see ). If $H=\nabla^{2}f(\bar{x})$, then $Q(H,x,\bar{x})$ is a second order approximation of $F(x)$ . In this paper we assume that $H$ is a positive definite matrix such that $\sigma I\preceq H\preceq MI$ for some positive constants $M$ and $\sigma$.

Minimizing the function $Q(H,u,v)$ over $u$ reduces to solving problem (1.2). We will use the following notation to denote the accurate and approximate solutions of (1.2).

The method that we consider in this paper computes iterates by (approximately) optimizing $Q(H,u,v)$ with respect to $u$ using some particular $H$ which is chosen at each iteration. The basic algorithm is described in Algorithms 1 and 2.

Algorithm 2 chooses Hessian approximations of the form $H_{k}=\frac{1}{\mu_{k}}I+G_{k}$. However, it is possible to consider any procedure of choosing positive definite $H_{k}$ which ensures $MI\succeq H_{k}\succeq\sigma I$ and $F(p_{H_{k}}(x))-F(x)\leq\rho(Q(H_{k},p_{H_{k}}(x),x)-F(x))$, for a given $0<\rho\leq 1$, - a step acceptance condition which is a relaxation of conditions used in and .

An inexact version of Algorithm 1 is obtained by simply replacing $p_{H_{k}}$ by $p_{H_{k},\phi_{k}}$ in both Algorithms 1 and 2 for some sequence of $\phi_{k}$ values.

Basic results, assumptions and preliminary analysis

ISTA is a particular case of Algorithm 1 with $G_{k}=0$, for all $k$, and $\rho=1$. In this case, the value of $\mu_{k}$ is chosen so that the conditions of Lemma 2 hold with $\epsilon=0$. In other words, the reduction achieved in the objective function $F(x)$ is at least the amount of reduction achieved in the model $Q(\mu_{k},p(x^{k}),x^{k})$. It is well known that as long as $\mu_{k}\leq 1/L(f)$ (recall that $L(f)$ is the Lipschitz constant of the gradient) then condition (3.8) holds with $\rho=1$. Relaxing condition (3.8) by using $\rho<1$ allows us to accept larger values of $\mu_{k}$, which in turn implies larger steps taken by the algorithm. This basic idea is the cornerstone of step size selection in most nonlinear optimization algorithms. Instead of insisting on achieving ”full” predicted reduction of the objective function (even when possible), a fraction of this reduction is usually sufficient. In our experiments small values of $\rho$ provided much better performance than values close to $1$.

In the next three sections we present the analysis of convergence rate of Algorithm 1 under different scenarios. Recall that we assume that $f(x)$ is convex and smooth, in other words $\|\nabla f(x)-\nabla f(y)\|\leq L(f)\|x-y\|$ for all $x$ and $y$ in the domain of interest, while $g(x)$ is simply convex. In Section 5.1 we assume that $g(x)=\lambda\|x\|_{1}$. Note that we do not assume that $f(x)$ is strongly convex or that it is twice continuously differentiable, because we do not rely on any accurate second order information in our framework. We only assume that the Hessian approximations are positive definite and bounded, but their accuracy can be arbitrary, as long as sufficient decrease condition holds. Hence we only achieve sublinear rate of convergence. To achieve higher local rates of convergence stronger assumptions on $f(x)$ and on the Hessian approximations have to be made, see for instance, and for related local convergence analysis.

First we present a helpful lemma which is a simple extension of Lemma 2 in to the case of general positive definite Hessian estimate. This lemma establishes some simple properties of an $\phi$-optimal solution to the proximal problem (1.2). It uses the concept of the $\phi$-subdifferential of a convex function $a$ at $x$, $\partial_{\phi}a(x)$, which is defined as the set of vectors $y$ such that $a(x)-y^{T}x\leq a(t)-y^{T}t+\phi$ for all $t$.

where $z=v-H^{-1}\nabla f(v)$. Then there exists $\eta$ such that $\frac{1}{2}\|\eta\|^{2}_{H^{-1}}\leq\phi$ and

Proof. (3.1) indicates that $p_{\phi}(v)$ is an $\phi$-minimizer of the convex function $a(x):=\frac{1}{2}\|x-z\|_{H}^{2}+g(x)$. If we let $a_{1}(x)=\frac{1}{2}\|x-z\|_{H}^{2}$ and $a_{2}(x)=g(x)$, then this is equivalent to

Then (3.2) follows using $z=v-H^{-1}\nabla f(v)$.

We will find the following bound useful on the norm of $\eta$ which follows from the above lemma.

where $\lambda_{max}$ is the largest eigenvalue of $H$ or its upper bound.

Below are the assumptions made in our analysis.

(i) The set of optimal solutions of (1.1), $X^{*}$, is nonempty and $x^{*}$ is any element of that set.

Without loss of generality, we restrict our discussions below to the level set ${\cal X}_{0}:={\cal X}_{F}(x^{0})$ given by some $x^{0}\in\operatorname{dom}(F)$, e.g., the initial iterate of the Algorithm 1.

(iii) $g$ is convex and Lipschitz continuous with constant $L_{g}$ for all $x,y\in{\cal X}_{0}$:

(iv) There exists positive constants $M$ and $\sigma$ such that for all $k\geq 0$, at the $k$-th iteration of Algorithm 1:

(v) There exists a positive constant $D_{{\cal X}_{0}}$ such that for all iterates $\{x^{k}\}$ of Algorithm 1:

In the analysis of the proximal gradient methods Assumption 1 is removed by directly establishing a uniform bound on $\|x^{k}-x^{*}\|$ (rf. ). In the next section we show an outline of a proximal-gradient type analysis where this assumption is imposed for simplicity of the presentation. It is possible to establish a similar, but more complex bound without this assumption. The proximal gradient approach in the next subsection, however, requires another, stronger, assumption on $H_{k}$. In the alternative analysis that follows we impose Assumption 1 but relax the assumption on $H_{k}$ matrices. Note that all iterates $\{x^{k}\}$ fall into the level set ${\cal X}_{0}$, due to the sufficient decrease condition (3.8) that demands a monotonic decrease on the objective values. Assumption 1 thus follows straightforwardly if the level set is bounded, which often holds in real-world problems or can be easily imposed.

In this section we extend the analysis in and to our framework under additional assumptions on the Hessian approximations $H_{k}$. The following lemma, is a generalization of Lemma 2.3 in and of a similar lemma in . This lemma serves to provide a bound on the change in the objective function $F(x)$.

Given $\epsilon$, $\phi$ and $H$ such that

Proof. The proof is an easy extension of that in .

Note that if $\phi=0$, that is the subproblems are solved accurately, then we have $2(F(u)-F(p(v)))\geq\|p(v)-u\|_{H}^{2}-\|v-u\|_{H}^{2}-2\epsilon$.

and Lemma 2 holds at each iteration $k$ of Algorithm 1 with $\epsilon_{k}=-\frac{1-\rho}{\rho}(F(x^{k+1})-F(x^{k}))$.

We now establish the sub linear convergence rate of Algorithm 1 under the following additional assumption.

Let $\{x^{k}\}$ be the sequence of iterates generated by Algorithm 1, then there exists a constant $M_{H}$ such that

where $M_{i}$ is the upper bound on the largest eigenvalue of $H_{i}$ as defined in Assumption 1.

The assumption above is not verifiable, hence it does not appear to be very useful. However, it is easy to see that a condition $\frac{H_{i+1}}{M_{i+1}}\preceq\frac{H_{i}}{M_{i}}$, for all $i$, easily implies (3.9) with $M_{H}=0$. This condition, in turn, is trivially satisfied if $H_{i+1}$ is a multiple of $H_{i}$ for all $i$, as it is in the case of ISTA algorithm. It is also clear that Assumption 2 is a lot weaker than the condition that $H_{i+1}$ is a multiple of $H_{i}$ for all $i$. For example, it also hold if $H_{i+1}$ is a multiple of $H_{i}$ for all $i$ except for a finite number of iterations. It also holds if $H_{i+1}$ converges to $H_{i}$ (or to its multiple) sufficiently rapidly. It is also weaker than the assumption in that $H_{i+1}\preceq{H_{i}}$. Exploring different conditions on the Hessian approximations that ensure Assumption 2 is a subject of a separate study. Below we show how under this condition sub linear convergence rate is established.

Suppose that Assumptions 1 and 2 hold. Assume that all iterates $\{x^{k}\}$ of inexact Algorithm 1 are generated with some $\phi_{k}\geq 0$ (cf. (2.3)), then

where $k_{m}=\sum_{i=0}^{k-1}M_{i}^{-1}$.

Proof. Let us apply Lemma 2, sequentially, with $u=x^{*}$, $p_{\phi}(v)=x^{i}$ and subproblem minimization residual $\phi_{i}$ for $i=0,\ldots,k-1$. Adding up resulting inequalities, we obtain

From (3.5) and definition of $M_{i}$ we know that $\|\eta_{i}\|\leq\sqrt{2M_{i}\phi_{i}}$. Using the already established bound $\sum_{i=0}^{k-1}\epsilon_{i}\leq\frac{(1-\rho)(F(x^{0})-F(x^{*}))}{\rho}$ we obtain

Let us consider the term $\frac{1}{k_{m}}=\frac{1}{\sum_{i=0}^{k-1}M_{i}^{-1}}$. From earlier discussions, we can see that $\frac{1}{k_{m}}\leq\frac{M}{k}$. Moreover, if $H_{k}$ are diagonal matrices, then $M=L(f)$ and, hence $\frac{1}{k_{m}}\leq\frac{L(f)}{k}$, which established a bound similar to that of proximal gradient methods. The role of $M_{i}$ is to show that if most of these values are much smaller than the global Lipschitz constant $L(f)$, then the constant involved in the sub linear rate can be much smaller than that of the proximal gradient methods. This is well known effect of using partial second order information and it is observed in our computational results.

We conclude that under Assumption 2 Algorithm 1 converges at the rate of $O(1/k)$ if $\sum_{i=0}^{k-1}\sqrt{\frac{2\phi_{i}}{M_{i}}}$ is bounded for all $k$. This result in similar to those obtained in . In it is shown how randomized block coordinate descent and other methods can be utilized to optimize subproblems $\min_{u}Q(H_{i},u,x^{i})$ so that $\sqrt{\frac{2\phi_{i}}{M_{i}}}$ decays sufficiently fast to guarantee such a bound (possibly in expectation).

In this paper, however, we focus on a different derivation of the sub linear convergence rate, which results in a different bound on $\phi_{k}$ and different, more complex, dependence on the constants, but, on the other hand, does not require Assumption 2 and results in a weaker assumption on $\phi_{i}$.

Analysis of sub linear convergence

In our analysis below we will use another known technique for establishing sub linear convergence of gradient descent type methods on smooth convex functions . However, due to the non smooth nature of our function the analysis requires significant extensions, especially in the inexact case. Moreover, it does not apply to the line-search algorithm, we will rely on the fact that a proximal quasi-Newton method is used in that each new iteration $x^{k+1}$ is an approximate minimizer of the function $Q(H_{k},u,x^{k})$. Our analysis, hence, also applies to proximal gradient methods.

First we prove the following simple result.

Consider $F(\cdot)$ defined in (1.1). Let Assumptions 1 hold. Then for any three points $u,v,w\in\operatorname{dom}(F)$, we have

where $\gamma_{g,\phi}^{v}\in\partial_{\phi}g(v)$ is any $\phi$-subgradient of $g(\cdot)$ at point $v$.

Proof. From convexity of $f$ and $g$ and the definition of $\phi$-subgradient, it follows that for any points $u,w$ and $v$,

Here we applied (4.2), (4.3) and (4.4) to get (4.5). Using Assumption 1 to bound the term $g(2v-u)-g(v)+g(u)-g(v)$ we can easily derive (4.1).

We now consider the exact version of Algorithm 1, i.e., $\phi_{k}=0$ for all $k$. We have the following lemma.

Moreover, there exists a vector $\gamma_{g}^{k+1}\in\partial g(x^{k+1})$ such that the following bounds hold:

Proof. The proof is a special case of Lemma 7, proved below.

Let Assumptions 1 hold for all $k$. Then the iterates $\{x^{k}\}$ generated by Algorithm 1 satisfy

Proof. We will denote $F(x^{k})-F^{*}$ by $\Delta F_{k}$. Our goal is to bound $\Delta F_{k}$ from above in terms of $1/k$ which we will achieve by deriving a lower bound on $\frac{1}{\Delta F_{k}}$ in terms of $k$.

for some constant $c_{k}$, which depends on iteration $k$, but will be lower bounded by a uniform constant.

This follows simply from (4.1) with $u=x^{k},w=x^{*},v=x^{k+1}$ and $\phi=0$,

Substituting the bounds on $\|x^{k+1}-x^{k}\|$ (cf. (4.7)) and $\|x^{k}-x^{*}\|$ (cf. Assumptions 1) in (4.11) we get the desired bound (4.10).

Indeed $F(x^{k})-F(x^{k+1})\geq\rho(Q(H_{k},x^{k},x^{k})-Q(H_{k},x^{k+1},x^{k}))$. And a bound on the reduction in $Q$ can be established by combining (4.6) and (4.7),

Finally, combining the lower bound on $F(x^{k+1})-F(x^{k})$ together with the upper bound on $\Delta F_{k}^{2}$ we can conclude that

which establishes (4.28) with $c_{k}=\frac{\rho\sigma_{k}^{3}}{2M_{k}^{2}(D_{{\cal X}_{0}}\sigma_{k}+2L_{g})^{2}}$.

Dividing both sides of the inequality above by $\Delta F_{k+1}\Delta F_{k}$ we have

Summing the above expression for $i=0,\ldots,k-1$ we have

Let us note that if $H_{k}=L(f)I$ for all $k$, as in standard proximal gradient methods, where $L(f)$ is the Lipschitz constant of $\nabla f(x)$, then the bound becomes

if $L_{g}\ll D_{{\cal X}_{0}}L(f)$. This bound is similar to $\frac{2\|x^{0}-x^{*}\|^{2}L(f)}{k}$ established for proximal gradient methods, assuming that $D_{{\cal X}_{0}}$ is comparable to $\|x^{0}-x^{*}\|$.

2 The inexact case

We now analyze Algorithm 1 in the case when the computation of $p_{H}(v)$ is performed inexactly. In other words, we consider the version of Algorithm 1 (and 2) where we compute $x^{k+1}:=p_{H_{k},\phi_{k}}(x^{k})$ and $\phi_{k}$ can be positive for any $k$. The analysis is similar to that of the exact case, with a few additional terms that need to be bounded. We begin by extending Lemma 5.

Moreover there exists a vector $\gamma_{g,\phi}^{k+1}\in\partial g_{\phi_{k}}(x^{k+1})$ such that the following bounds hold:

Proof. Recall Lemma 1, from (3.2), there exists a vector, which we will refer to as $\gamma_{g,\phi}^{k+1}$, such that

with $\frac{1}{2}\|\eta_{k}\|^{2}_{H_{k}^{-1}}\leq\phi_{k}$, which, in turn, implies $\|\eta_{k}\|\leq\sqrt{2M_{k}\phi_{k}}$. The following inequality follows from the definition of $\phi$-subdifferential,

Using (4.16) and (4.17), we derive (4.13) as follows,

Unlike the exact case, the inquality $\frac{1}{\Delta F_{k+1}}-\frac{1}{\Delta F_{k}}\geq c_{k}$ can no longer be guaranteed to hold on each iteration. The convergence rate is obtained by observing that when this desired inequality fails another inequality always holds, which bounds $\Delta F_{k}$ in terms of $\phi_{k}$. Specifically, we have the following theorem.

Consider $k$th iteration of the inexact Algorithm 1 with $0\leq\phi_{k}\leq 1$. Let $\Delta F_{k}:=F(x^{k})-F(x^{*})$. Then there exists large enough positive constant $\theta>0$, such that one of the following two cases must hold,

where $b_{k}$ and $c_{k}$ are given below,

Proof. First, applying (4.1) with $u=x^{k},w=x^{*}$ and $v=x^{k+1}$ we obtain,

We will consider two cases that are possible at each iteration $k$ for some fixed constant $\theta>1$ which we will specify later.

Let us assume that Case 1 holds, then from (4.14) and (4.22) it simply follows that

Using (4.22), (4.21), the bound on $\|x^{k+1}-x^{k}\|$ from (4.24) together with the bound on $\|x^{k}-x^{*}\|$ from Assumptions 1 we get

We now consider Case 2, where (4.23) along with (4.14) from Lemma 7 imply

Substituting into (4.21) the upper bound on $\|x^{k+1}-x^{k}\|$ from (4.26) and the bound on $\|x^{k}-x^{*}\|$ from Assumptions 1 we now get

From $\phi_{k}\leq 1$ it follows that $\|\nabla f(x^{k})+\gamma_{g,\phi}^{k+1}\|\geq\theta\sqrt{2M_{k}\phi_{k}}\geq\theta\sqrt{2M_{k}}\phi_{k}$. Hence we obtain

We will show that in this case, as in the exact case, we have

for some constant $c_{k}$ (different from that in the exact case). Towards that goal we will establish a lower bound on $F(x^{k})-F(x^{k+1})$ in terms of $\|\nabla f(x^{k})+\gamma_{g,\phi}^{k+1}\|^{2}$ as before. We still have $F(x^{k})-F(x^{k+1})\geq\rho(Q(H_{k},x^{k},x^{k})-Q(H_{k},x^{k+1},x^{k}))$. We now use the bounds (4.13) from Lemma 7, (4.23) and (4.26) and obtain

By selecting a sufficiently large $\theta$ we can ensure that

for $t_{k}=\frac{\sigma_{k}}{2}(\frac{\theta-1}{\theta M_{k}})^{2}-\frac{1+\theta}{\theta^{2}\sigma_{k}}-\frac{1}{2\theta^{2}M_{k}}>0$. Finally, combining the lower bound (4.29) on $F(x^{k+1})-F(x^{k})$ together with the upper bound (4.27) on $\Delta F_{k}^{2}$ we can conclude that (4.28) holds with

Finally, (4.19) follows from (4.28) divided by $\Delta F_{k}\Delta F_{k+1}$ and using the fact that $\frac{\Delta F_{k}}{\Delta F_{k+1}}\geq 1$.

Let us discuss the result of the above lemma. The lemma applies for any value of $\theta$ for which $t_{k}$, and hence, $c_{k}$ is positive for all $k$. It is easy to see that large values of $\theta$ imply large values of $c_{k}$. On the other hand, large $\theta$ is likely to cause Case 1 to hold (i.e., (4.18)) instead of Case 2 (i.e., (4.19)) on any given iteration, with $b_{k}$ also growing with the size of $\theta$. As we will show below the overall rate of convergence of the algorithm is derived using the two bounds - (4.18), where the rate is controlled by the rate of $\phi_{k}\to 0$ and (4.19), which is similar to the bound in the exact case. The overall bound, thus, will depend on the upper bound on $b_{k}$’s and the inverse of the lower bound on $c_{k}$’s. If, again, we assume that $\sigma_{k}=M_{k}=L(f)$ for all $k$, then $\theta=O(\sqrt{L(f)})$ is sufficient to ensure that $c_{k}>0$ and this results in $b_{k}\leq O(D_{{\cal X}_{0}}L(f))$ and $1/c_{k}\geq O(D_{{\cal X}_{0}}^{2}L(f))$, thus again, we obtain a bound which is comparable to that of proximal gradient methods, although with more complex constants.

We now derive the overall convergence rate, under the assumption that $\phi_{k}$ decays sufficiently fast.

Suppose that Assumption 1 holds. Assume that all iterates $\{x^{k}\}$ of inexact Algorithm 1 are generated with some $\phi_{k}\geq 0$ that satisfy

Let $\theta$ be chosen as specified in Lemma 8. Then for any $k$

Proof. Consider all iterations before a particular iteration $k$. From Lemma 8, it follows that either (4.18) or (4.19) must hold for each prior iteration. Let $k_{1}<k$ denote the index of the last iteration before $k$, for which (4.18) holds. If no such $k_{1}$ exists, then (4.19) holds for all $k$ and without loss of generality we can consider $k_{1}=0$. From (4.18) and from the fact that the function value never increases

For each iteration from $k_{1}+1$ to $k-1$, (4.19) gives

Summing up the above inequalities and using (4.34) we obtain the following bound on $\Delta F_{k}$,

To derive a simple uniform bound on $\Delta F_{k}$ we will use $b$ and $c$ - uniform upper and lower bounds, respectively, for $b_{k}$ and $c_{k}$ given in (4.20), i.e., $b_{k}\leq b$, $c_{k}\geq c,\forall k\geq k_{0}$. From Assumptions 1 we can derive the expressions for $c$ as follows,

Bound $b$ can be obtained in a similar fashion.

Substituting bounds $b$ and $c$ in (5.5) we get

It follows that the inexact version of Algorithm 1 has sublinear convergence rate if $\phi_{i}\leq a^{2}/i^{2}$ for some $a<1$ and all iterations $i=0,\ldots,k$. In contrast, the bounds in and in Section 3.1 require that $\sum_{i=0}^{\infty}\sqrt{\phi_{i}}$ is bounded. This bound on the overall sequence is clearly stronger than $\phi_{i}\leq a^{2}/i^{2}$, since $\sum_{i=0}^{\infty}\frac{a}{i}=\infty$. On the other hand, it does not impose any particular requirement on any given iteration, except that each $\phi_{i}$ is finite, which our bound on $\phi_{i}$ is assumed to hold at each iteration, so far.

3 Complexity in terms of subproblem solver iterations

Let us discuss conditions on the sequence of $\phi_{i}$ established above and how they can be ensured. Firstly, let us note that condition $a<1$ in Theorem 10 can easily be removed. We introduced it for the sake of brevity, to ensure that $\phi_{i}\leq 1$ on each iteration. Clearly, an arbitrarily large $a$ can be used and in that case $\phi_{i}\leq 1$ for all $i\geq 1/a$. Moreover, the condition $\phi_{i}\leq 1$ is only needed to replace $\phi_{i}$ with $\sqrt{\phi_{i}}$ in Lemma 8 in inequality (4.25). Instead we can use bound (4.23) and upper bounds on $\nabla f(x^{i})$ and $\gamma^{i+1}_{g,\phi}$ to replace $\phi_{k}$ with a constant multiple of $\sqrt{\phi_{i}}$. In conclusion, it is sufficient to solve the subproblem on the $i$-th iteration to accuracy $O(1/i^{2})$.

The question now is: what method and what stopping criterion should be used for subproblem optimization, so that sufficient accuracy is achieved and no excessive computations are performed, in other word, how can we guarantee the bound on $\phi_{i}$, while maintaining efficiency of the subproblem optimization? It is possible to consider terminating the optimization of the $i$-th subproblem once the duality gap is smaller than the required bound on $\phi_{i}$. However, checking duality gap can be computationally very expensive. Alternatively one can use an algorithm with a known convergence rate. This way it can be determined apriori how many iterations of such an algorithm should be applied to the $i$-th subproblem to achieve the desired accuracy. In particular, we note that the objective functions in our subproblems are all $\sigma$-strongly convex, so a simple proximal gradient method, or some of its accelerated versions, will enjoy linear convergence rates when applied to these subproblems. Hence, after $i$ iterations of optimizing $Q_{i}$, such a method will achieve accuracy $\phi_{i}$ that decays geometrically, i.e., $\phi_{i}=C\delta^{i}$, for some constants $C>0$ and $0<\delta<1$, hence $\sum_{i=0}^{\infty}\sqrt{\phi_{i}}$ is bounded. Note that the same property holds for any linearly convergent method, such as the proximal gradient or a semi-smooth Newton method. Also, it is easy to see that $\phi_{i}\leq a^{2}/i^{2}$ holds for some $a>0$ for all $i$. One can also can use FISTA to optimize $Q_{i}$ which will ensure $\phi_{i}\leq a^{2}/i^{2}$ for some $a>0$ but will not guarantee $\sum_{i=0}^{\infty}\sqrt{\phi_{i}}$. The advantage of using FISTA and its resulting rate is that ist does not depend on the strong convexity constant, hence the subproblem complexity does not depend on $\sigma$ - the lower bound on the smallest eigenvalues of the Hessian approximations. In conclusion, we have the following overall bounds.

Suppose that Assumptions 1 hold and that at the $k$-th iteration of inexact Algorithm 1 function $Q(H_{k},u,x^{k})$ is approximately minimized, to obtain $x^{k+1}$ by applying $l(k)=\alpha k+\beta$ steps of any algorithm which guarantees that $Q(H_{k},x^{k+1},x^{k})\leq Q(H_{k},x^{k},x^{k})$ and whose convergence rate ensures the error bound $\phi_{k}\leq a^{2}/(\alpha k+\beta)^{2}$ for some $a>0$. Then accuracy $F(x^{k})-F(x^{*})\leq\epsilon$ is achieved after at most

inner iterations (of the chosen algorithm), with $b,c$ given in Theorem 10.

Proof. A proof follows trivially from Theorem 10.

Suppose that Assumptions 1 hold and that at the $k$-th iteration of inexact Algorithm 1 function $Q(H_{k},u,x^{k})$ is approximately minimized, to obtain $x^{k+1}$ by applying $l(k)$ steps of an algorithm, which guarantees that $Q(H_{k},x^{k+1},x^{k})\leq Q(H_{k},x^{k},x^{k})$ and whose convergence rate ensures the error bound $\phi_{k}\leq\delta^{l(k)}M_{Q}$, for some constants $0<\delta<1$ and $M_{Q}>0$. Then, by setting $l_{k}=2log_{\frac{1}{\delta}}(k)$, accuracy $F(x^{k})-F(x^{*})\leq\epsilon$ is achieved after at most

inner iterations (of the chosen algorithm), with $t=\lceil\frac{\max\{ba,\frac{1}{c}\}}{\epsilon}+1\rceil$ and $b,c$ given in Theorem 10.

Proof. A proof follows trivially from Theorem 10.

The total complexity in terms of the inner iterations should not be viewed as a summary of the whole complexity of Algorithm 1. A key step of the algorithm is the computation of $F(x^{k})$ and $\nabla f(x^{k})$ at each iteration. In big data applications this is often the most extensive step, hence the main complexity is defined by the number of function and gradient computations. Due to backtracking via proximal parameter update to satisfy sufficient decrease condition, the number of function and gradient computation may be larger than the number of iterations of Algorithm 1, however, it does not exceed this number by more than a logarithmic factor. In practice, only several initial iterations contain backtracking steps, hence Theorem 10 provides the bound on the complexity in terms of function and gradient computations.

In the next section we extend our convergence rate results to the case of solving subproblems via randomized coordinate descent, where $\phi$ is random and hence does not satisfy required bounds on each iteration.

Analysis of the inexact case under random subproblem accuracy

As we pointed out in the introduction, the most efficient practical approach to subproblem optimization, in the case when $g(x)=\lambda\|x\|_{1}$, seems to be the coordinate descent method. One iteration of a coordinate descent step can be a lot less expensive than that of a proximal gradient or a Newton method. In particular, if matrix $H$ is constructed via the LBFGS approach, then one step of a coordinate decent method takes a constant number of operations, $m$ (the memory size of LBFGS, which is typically 10-20). On the other hand, one step of proximal gradient takes $O(mn)$ operations and Newton method takes $O(nm^{2})$.

Unfortunately, cyclic (Gauss-Seidel) coordinate descent, does not have deterministic complexity bounds, hence it is not possible to know when the work on a particular subproblem can be terminated to guarantee the desired level of accuracy. However, a randomized coordinate descent has probabilistic complexity bounds, which can be used to demonstrate the linear rate of convergence in expectation.

We have the following probabilistic extension of Theorem 10.

Suppose that Assumption 1 holds. Assume that for all $k$ iterates $\{x^{k}\}$ of inexact Algorithm 1 are generated with some $\phi_{k}\geq 0$ that satisfy

conditioned on the past. Let $\theta$, $b$ and $c$ be as specified in Theorem 10. Then for any $k$

Proof. As in the proof of Theorem 10 consider all iterations before a particular iteration $k$. From Lemma 8, it follows that either (4.18) or (4.19) must hold for each prior iteration. Let $k_{1}<k$ denote the index of the last iteration before $k$, for which (4.18) holds. Let $k_{2}$ denote the index of the second to last such iteration and so on, hence $k_{i}$ is the index of the last iteration such that there are exactly $i$ iterations between $k_{i}$ and $k-1$ for which (4.18) holds. Without loss of generality we can assume that $k_{i}$ exists for each $i$, because if it does not - we can set $k_{i}=0$ and obtain a better bound. Let us now assume for a given $i$ that $\phi_{k_{i}}\leq\frac{a^{2}}{{k_{i}}^{2}}$ holds, but that $\phi_{k_{j}}>\frac{a^{2}}{{k_{j}}^{2}}$ for all $j=1,\ldots,i-1$ (if $i=1$ we have the case analyzed in the proof of Theorem 10). The analysis in the proof of Theorem 10 extends easily to the case $i>1$ by observing that

holds for any $k_{i}+1\leq l\leq k-1$, $l\neq k_{j},\,j=1\ldots,i-1$, that is all the iterations for which or (4.18) does not hold. We also have

for $k_{i}+1\leq l\leq k-1$, $l=k_{j},\,j=1\ldots,i-1$, simply from the fact that function values never increase. Summing up the above inequalities and using the fact that $\phi_{k_{i}}\leq\frac{a^{2}}{{k_{i}}^{2}}$ we obtain $\Delta F_{k}$,

Now, recall that $P\{\phi_{k_{i}}\leq\frac{a^{2}}{{k_{i}}^{2}}\}\geq 1-p$, for any iteration $i$, independently of the other iteration. This means that the probability that $\{\phi_{k_{i}}\leq\frac{a^{2}}{{k_{i}}^{2}}\}$ and $\{\phi_{k_{j}}>\frac{a^{2}}{{k_{j}}^{2}}\}$ for all $j=1\ldots,i-1$ is $(1-p)p^{i}$. This implies

To bound the term $\sum_{i=1}^{k-1}(p^{i-1}+\frac{i-1}{k-i}p^{i-1})$ observe that $\frac{i-1}{k-i}p^{i-1}\leq(i-1)p^{i-1}$ and hence

which gives us the final bound on the expected error.

We note that the (2-p) factor is the result of an overestimate of the weighted geometric series and a tighter bound should be possible to obtain.

Below we show that randomized coordinate descent can guarantee sufficient accuracy for subproblem solutions and hence maintain the sub linear convergence rate (in expectation) of Algorithm 1. Moreover, we show in Section 7, that the randomized coordinate descent is as efficient in practice as the cyclic one.

In randomized coordinate descent the model function $Q(\cdot)$ is iteratively minimized over one randomly chosen coordinate, while the others remain fixed. The method is presented in Algorithm 3 and is applied for $l$ steps, with $l$ being an input parameter.

Here we will show how properties of coordinate descent can be used together with the analysis in Section 4. Combination of coordinate descent with the waker analysis presented in Section 3.1 can be found in .

Our analysis is based on Richtarik and Takac’s results on iteration complexity of randomized coordinate descent . In particular, we make use of Theorem 7 in , which we restate below without proof, while adapting it to our context.

where $\mu(H)$ is a constant that measures conditioning of $H$ along the coordinate directions and in the worst case is at most $M/\sigma$ - the condition number of $H$.

Let us now state the version of Algorithms 1 and 2 which is close to what implement in practice and discuss in the following sections and for which the complexity bound can be applied.

The conclusion of Lemma 16 in application to Algorithms 3-5 is that when applying $l(k)$ steps of randomized coordinate descent approximately to optimize $Q(H_{k},u,x^{k})$, $\phi_{k}\leq M_{Q}\delta^{l(k)}$, with probability $p$, where $M_{Q}$ is an upper bound on $\frac{Q(H_{k},x^{k},x^{k})-Q(H_{k},p_{H}(x^{k}),x^{k})}{p}$ and $\delta=e^{-\frac{1}{n(1+\mu(H_{k}))}}$. This, together with Theorem 15 implies that to solve subproblem on iteration $k$ it is sufficient to set $l(k)=O(n(1+\mu(H))\log(kp/M_{Q}))$ and the resulting convergence rate will then obey Theorems 10 and 13. However, it is necessary to know constants $M_{Q}$ and $\mu(H)$ to be able to construct efficient expression $l(k)$. In practice, a successful strategy is to select a slow growing linear function of $k$, $l(k)=ak+b$. This certainly guarantees convergence rate of the outer iteration as in Theorem 15. In terms of overall rate this gives inferior complexity, however, we believe that the real difference in terms of the workload appears only in the limit, while in most cases the algorithm successfully terminated before our practical formula $l(k)=ak+b$, described in the next section, significantly exceeds, the theoretical bound $O(n(1+\mu(H))log(kp/M_{Q}))$ with appropriate constants. Moreover, as noted earlier, the number of function and gradient evaluations may be the dominant complexity in the big data cases, hence it may be worthwhile to increase workload of coordinate descent in order to reduce the constants in the bound in Theorem 15.

Finally, we note that when using ISTA method for the subproblem, instead of randomized coordinate descent, the number of inner iteration does not need to depend on the dimension $n$. However, it depends on $\sigma_{k}/M_{k}$ and the cost per iteration is roughly $n$ times bigger than that of coordinate descent (with LBFGS matrices). Hence the overall complexity of using coordinate descent is better than that of ISTA if $\mu(H_{k})\ll\sigma_{k}/M_{k}$. This indeed happens often in practical problems, as is discussed in and other works on coordinate descent.

Optimization Algorithm

In this section we briefly describe the specifics of the general purpose algorithm that we propose within the framework of Algorithms 4, 5 and 3 and that takes advantage of approximate second order information while maintaining low complexity of subproblem optimization steps. The algorithm is designed to solve problems of the form (1.1) with $g(x)=\lambda\|x\|_{1}$, but it does not use any special structure of the smooth part of the objective, $f(x)$.

At iteration $k$ a step $d_{k}$ is obtained, approximately, as follows

with $H_{k}=G_{k}+\frac{1}{2\mu_{k}}I$ - a positive definite matrix and ${\cal A}_{k}$ - a set of coordinates fixed at the current iteration.

The positive definite matrix $G_{k}$ is computed by a limited memory BFGS approach. In particular, we use a specific form of Hessian estimate, (see e.g. ),

where $Q$, $\gamma_{k}$ and $R$ are defined below,

Note that there is low-rank structure present in $G_{k}$, the matrix given by $Q\hat{Q}$, which we can exploit, but $G_{k}$ itself by definition is always positive definite. Let $m$ be a small integer which defines the number of latest BFGS updates that are ”remembered” at any given iteration (we used $10-20$). Then $S_{k}$ and $T_{k}$ are the $p\times m$ matrices with columns defined by vector pairs $\{s_{i},t_{i}\}_{i=k-m}^{k-1}$ that satisfy $s_{i}^{T}t_{i}>0,s_{i}=x^{i+1}-x_{i}$ and $t_{i}=\nabla f(x^{i+1})-\nabla f(x^{i})$, $M_{k}$ and $D_{k}$ are the $k\times k$ matrices

The particular choice of $\gamma_{k}$ is meant to promote well-scaled quasi-Newton steps, so that less time is spent on line search or updating of prox parameter $\mu_{k}$ . In fact instead of updating and maintaining $\mu_{k}$, exactly as described in Algorithm 2 we simply double $\gamma_{k}$ in (6.1) at each backtracking step. This can be viewed as choosing $\mu_{k}=\infty$ for the first step of backtracking, and $\mu_{k}=1/(2^{i-1}-1)\gamma_{k}$ for the $i$-th backtracking step, when $i>1$. As long as $G_{K}$ in (6.1), is positive definite, with smallest eigenvalue bounded by $\sigma>0$, our theory applies to this particular backtracking procedure.

An active-set selection strategy maintains a sequence of sets of indices ${\cal A}_{k}$ that iteratively estimates the optimal active set ${\cal A}^{*}$ which contains indices of zero entries in the optimal solution $x^{*}$ of (1.1). We introduce this strategy as a heuristic aiming to improve the efficiency of the implementation and to make it comparable with state-of-the-art methods, which also use active set strategies. A theoretical analysis of the effects of these strategies is a subject of future study. The complement set of ${\cal A}_{k}$ is ${\cal I}_{k}=\{i\in{\cal P}~{}|~{}i\notin{\cal A}_{k}\}$. Let $(\partial F(x^{k}))_{i}$ be the $i$-th component of a subgradient of $F(x)$ at $x^{k}$. We define two sets,

As is done in and we select ${\cal I}_{k}$ to include the entire set ${\cal I}^{(2)}_{k}$ and the entire set ${\cal I}^{(1)}_{k}$. We also tested a strategy which includes only a small subset of indices from ${\cal I}^{(1)}_{k}$ for which the corresponding elements $|(\partial F(x^{k}))_{i}|$ are the largest. This strategy resulted in a smaller size of subproblems (6) at the early stages of the algorithm, but did not appear to improve the overall performance of the algorithm.

2 Solving the inner problem via coordinate descent

We apply coordinate descent method to the piecewise quadratic subproblem (6) to obtain the direction $d_{k}$ and exploit the special structure of $H_{k}$. Suppose $j$-th coordinate in $d$ is updated, hence $d^{\prime}=d+ze_{j}$ ($e_{j}$ is the $j$-th vector of the identity). Then $z$ is obtained by solving the following one-dimensional problem

which has a simple closed-form solution .

The most costly step of an iteration of the coordinate descent method is computing or maintaining vector $H_{k}d$. Naively, or in the case of general $H_{k}$, this step takes $O(n)$ flops, since the vector needs to be updated at the end of each iteration, when one of the coordinates of vector $d$ changes. The special form of $G_{k}$ in $H_{k}=G_{k}+\frac{1}{\mu}I=\gamma_{k}I-Q\hat{Q}$ gives us an opportunity to accelerate this step, reducing the complexity from problem-dependent $O(n)$ to $O(m)$ with $m$ chosen as a small constant. In particular we only store the diagonal elements of $G_{k}$, $(G_{k})_{ii}=\gamma_{k}-q_{i}^{T}\hat{q}_{i}$, where $q_{i}$ is the $i$th row of the matrix $Q$ and $\hat{q}_{i}$ is the $i$th column vector of the matrix $\hat{Q}$. We compute $(G_{k}d)_{i}$, whenever it is needed, by maintaining a $2m$ dimensional vector $v:=\hat{Q}d$, which takes $O(2m)$ flops, and using $(G_{k}d)_{i}=\gamma_{k}d_{i}-q_{i}^{T}v$. After each coordinate step $v$ is updated by $v\leftarrow v+z_{i}\hat{q}_{i}$, which costs $O(m)$. We also need to use extra memory for caching $\hat{Q}$ and $\hat{d}$ which takes $O(2mp+2m)$ space. With the other $O(2p+2mn)$ space for storing the diagonal of $G_{k}$, $Q$ and $d$, altogether we need $O(4mp+2n+2m)$ space, which is essentially $O(4mn)$ when $n\gg m$.

Computational experiments

The aim of this section is to provide validation for our general purpose algorithm, but not to conduct extensive comparison of various inexact proximal Newton approaches. In particular, we aim to demonstrate a) that using the exact Hessian is not necessary in these methods, b) that backtracking using prox parameter, based on sufficient decrease condition, which our theory uses, does in fact work well in practice and c) that randomized coordinate descent is at least as effective as the cyclic one, which is standardly used by other methods.

QUIC: the quadratic inverse covariance algorithm for solving SICS described in .

LIBLINEAR: an improved version of GLMNET for solving SLR described in .

Note that both of these packages have been shown to be the state-of-the-art solvers in their respective categories (see e.g. ).

Both QUIC and LIBLINEAR adopt line search to ensure function reduction. We have implemented line search in LHAC as well to see how it compares to the updating of prox parameter proposed in Algorithm 2. In all the experiments presented below use the following notation.

LHAC: Algorithm 4 with backtracking on prox parameter.

LHAC-L: Algorithm 4 with Armijo line search procedure described below in (7.4).

For all of the experiments we choose the initial point $x_{0}=\mathbf{0}$, and we report running time results in seconds, plotted against log-scale relative objective function decrease given by

where $F^{*}$ is the optimal function value. Since $F^{*}$ is not available, we compute an approximation by setting a small optimality tolerance, specifically $10^{-7}$, in QUIC and LIBLINEAR. All the experiments are executed through the MATLAB mex interface. We also modify the source code of LIBLINEAR in both its optimization routine and mex gateway function to obtain the records of function values and the running time. We note that we simply store, in a double array, and pass the function values which the algorithm already computes, so this adds little to nothing to LIBLINEAR’s computational costs. We also adds a function call of clock() at every iteration to all the tested algorithms, except QUIC, which includes a “trace” mode that returns automatically the track of function values and running time, by calling clock() iteratively. For both QUIC and LIBLINEAR we downloaded the latest versions of the publicly available source code from their official websites, compiled and built the software on the machine on which all experiments were executed, and which uses 2.4GHz quad-core Intel Core i7 processor, 16G RAM and Mac OS.

The optimal objective values $F^{*}$ obtained approximately by QUIC and LIBLINEAR are later plugged in LHAC and LHAC-L to terminate the algorithm when the following condition is satisfied

In LHAC we chose $\bar{\mu}=1,\beta=1/2$ and $\rho=0.01$ for sufficient decrease (see Algorithm 5), and for LBFGS we use $m=10$.

When solving the subproblems, we terminate the RCD procedure whenever the number of coordinate steps exceeds

where $|{\cal I}_{k}|$ denotes the number of coordinates in the current working set. Condition (7.3) indicates that we expect to update each coordinate in ${\cal I}_{k}$ only once when $k<m$, and that when $k>m$ we increase the number of expected passes through ${\cal I}_{l}$ by 1 every $m$ iterations, i.e., after LBFGS receives a full update. The idea is not only to avoid spending too much time on the subproblem especially at the beginning of the algorithm when the Hessian approximations computed by LBFGS are often fairly coarse, but also to solve the subproblem more accurately as the iterate moves closer to the optimality. Note that in practice when $|{\cal I}_{k}|$ is large, the value of (7.3) almost always dominates $k$, hence it can be lower bounded by $l(k)=ak+b$ with some reasonably large values of $a$ and $b$, which, as we analyzed in Section 5.1, guarantees the sub linear convergence rate. We also find that (7.3) works quite well in practice in preventing from “over-solving” the subproblems, particularly for LBFGS type algorithms. In Figures 2 we plot the data with respect to the number of RCD iterations. In particular Figures 2(a) and 2(c) show the number of RCD steps taken at the $k$-th iteration, as a function of $k$. Figures 2(b) and 2(d) show convergence of the objective function to its optimal value as a function of the total number of RCD steps taken so far (both values are plotted in logarithmic scale). Note that RCD steps are not the only component of the CPU time of the algorithms, since gradient computation has to be performed at least once per iteration.

In LHAC-L, a line search procedure is employed, as is done in QUIC and LIBLINEAR, for the convergence to follow from the framework by . In particular, the Armijo rule chooses the step size $\alpha_{k}$ to be the largest element from $\{\beta^{0},\beta^{1},\beta^{2},...\}$ satisfying

where $0<\beta<1,0<\sigma<1$, and $\Delta_{k}:=\nabla f_{k}^{T}d_{k}+\lambda\|x_{k}+d_{k}\|_{1}-\lambda\|x_{k}\|_{1}$. In all the experiments we chose $\beta=0.5,\sigma=0.001$ for LHAC-L.

2 Sparse Inverse Covariance Selection

The sparse inverse covariance selection problem is defined by

For SICS we report results on four real world data sets, denoted as ER_692, Arabidopsis, Leukemia and hereditarybc, which are preprocessed from breast cancer data and gene expression networks. We refer to for detailed information about those data sets.

We set the regularization parameter $\lambda=0.5$ for all experiments as suggested in . The plots presented in Figure 1 show that LHAC and LHAC-L is almost twice as fast as QUIC, in the two largest data sets Leukemia and hereditarybc (see Figure 1(c) and 1(d)). In the other two smaller data sets the results are less clear-cut, but all of the methods solve the problems very fast and the performance of LHAC is comparable to that of QUIC. The performances of LHAC and LHAC-L are fairly similar in all experiments. Again we should note that with the sufficient decrease condition proposed in Algorithm 2 we are able to establish the global convergence rate, which has not been shown in the case of Armijo line search.

3 Sparse Logistic Regression

The objective function of sparse logistic regression is given by

We report results of SLR on four data sets downloaded from UCI Machine Learning repository , whose statistics are summarized in Table 1. In particular, the first data set is the well-known UCI Adult benchmark set a9a used for income classification, determining whether a person makes over $50K/yr or not, based on census data; the second one we use in the experiments is called epsilon, an artificial data set for PASCAL large scale learning challenge in 2008; the third one, slices, contains features extracted from CT images and is often used for predicting the relative location of CT slices on the human body; and finally we consider gisette, a handwritten digit recognition problem from NIPS 2003 feature selection challenge, with the feature set of size 5000 constructed in order to discriminate between two confusable handwritten digits: the four and the nine.

The results are shown in Figure 3. In most cases LHAC and LHAC-L outperform LIBLINEAR. On data set slice, LIBLINEAR experiences difficulty in convergence which results in LHAC being faster by an order of magnitude. On the largest data set epsilon, LHAC and LHAC-L is faster than LIBLINEAR by about one third and reaches the same precision. Finally we note that the memory usage of LIBLINEAR is more than doubled compared to that of LHAC and LHAC-L, as we observed in all the experiments and is particularly notable on the largest data set epsilon.

Conclusion

In this paper we presented analysis of global convergence rate of inexact proximal quasi-Newton framework, and showed that randomized coordinate descent, among other subproblem methods, can be used effectively to find inexact quasi-Newton directions, which guarantee sub linear convergence rate of the algorithm, in expectation. This is the first global convergence rate result for an algorithm that uses coordinate descent to inexactly optimize subproblems at each iteration. Moreover, we improve upon results for inexact proximal gradient method in in that our requirements on the error in subproblem solution are weaker and the resulting bound on the total number of inner solver iterations is smaller.

Our framework does not rely on or exploit the accuracy of second order information, and hence we do not obtain fast local convergence rates. We also do not assume strong convexity of our objective function, hence a sublinear conference rate is the best global rate we can hope to obtain. In an accelerated scheme related to our framework is studied and an optimal sublinear convergence rate is shown, but the assumptions on the Hessian approximations are a lot stronger in than in our paper, hence the accelerated method is not as widely applicable. The framework studied by us in this paper covers several existing efficient algorithms for large scale sparse optimization. However, to provide convergence rates we had to depart from some standard techniques, such as line-search, replacing it instead by a prox-parameter updating mechanism with a trust-region-like sufficient decrease condition for acceptance of iterates. We also use randomized coordinate descent instead of a cyclic one. We demonstrated that this modified framework is, nevertheless, very effective in practice and is competitive with state-of-the-art specialized methods.

References