Proximal Quasi-Newton Methods for Regularized Convex Optimization with Linear and Accelerated Sublinear Convergence Rates

Introduction

We address the convex optimization problem of the form

where $L$ is the global Lipschitz constant of the gradient $\nabla f$. This class of problems, when $g(x)=\lambda\|x\|_{1}$, contains some of the most common machine learning models such as sparse logistic regression lin ; shalev , sparse inverse covariance selection hsieh ; olsen ; rish , and unconstrained Lasso tibshirani .

The Proximal Gradient Algorithm (PGA) is a variant of the proximal methods and is a well-known first-order method for solving optimization problem (1). Although classical subgradient methods can be applied to problem (1) when $g$ is nonsmooth, they can only achieve the rate of convergence of $\mathcal{O}(1/\sqrt{k})$ nesterov2 , while PGA converges at a rate of $\mathcal{O}(1/k)$ in both smooth and nonsmooth cases nesterov4 ; beck . In order to improve the global sublinear rate of convergence of PGA further, the Accelerated Proximal Gradient Algorithm (APGA) has been originally proposed by Nesterov in nesterov3 , and refined by Beck and Teboulle in beck . It has been shown that the APGA provides a significant improvement compared to PGA, both theoretically, with a rate of convergence of $\mathcal{O}({1}/{k^{2}})$, and practically beck . This rate of convergence is known to be the best that can be obtained using only first-order information yudin ; nesterov2 ; goldfarb , causing APGA to be known as an optimal first-order method. The class of accelerated methods contains many variants that share the same convergence rates and use only first-order information nesterov2 ; nesterov5 ; tseng . The main known drawback of most of the variants of APGA is that the sequence of the step-size $\{\mu_{k}\}$ has to be nonincreasing. This theoretical restriction sometimes has a big impact on the performance of this algorithm in practice. In goldfarb , in order to overcome this difficulty, a new version of APGA has been proposed. This variant of APGA allows to increase step-sizes from one iteration to the next, but maintain the same rate of convergence of $\mathcal{O}({1}/{k^{2}})$. In particular, the authors have shown that a full backtracking strategy can be applied in APGA and that the resulting complexity of the algorithm depends on the average value of step-size parameters, which is closely related to local Lipschitz constants, rather than the global one.

To make PGA and APGA practical, for some complicated instances of (1), one needs to allow for inexact computations in the algorithmic steps. In mark , inexact variants of PGA and APGA have been analyzed with two possible sources of error: one in the calculation of the gradient of the smooth term and the other in the proximal operator with respect to the nonsmooth part. The convergence rates are preserved if the sequence of errors converges to zero sufficiently fast. Moreover, it has been shown that both of these algorithms obtain a linear rate of convergence, when the smooth term $f$ is strongly convexFor APGA a different variant is analyzed in the case of strong convexity.. Recently, in dmitri , the linear convergence of PGA has been shown under the quadratic growth condition, which is weaker than a strong convexity assumption. In particular, their analysis relies on the fact that PGA linearly bounds the distance to the solution set by the step lengths. This property, called an error bound condition, has been proved to be equivalent to the standard quadratic growth condition. More precisely, a strong convexity assumption is a sufficient, but not a necessary condition for this error bound property.

While PGA and APGA can be efficient in solving (1), it has been observed that using (partial) second-order information often significantly improves the performance of the algorithms. Hence, Newton type proximal algorithms, also known as the proximal Newton methods, have become popular sra ; olsen ; lee ; nocedal3 ; tang and are often the methods of choice. When accurate (or nearly accurate) second-order information is used, the method no longer falls in the first-order category and faster convergence rates are expected, at least locally. Indeed, the global convergence and the local superlinear rate of convergence of the Proximal Quasi-Newton Algorithm (PQNA) are presented in lee and nocedal3 , respectively for both the exact and inexact settings. However, in the case of limited memory BFGS method nocedal1 ; nocedal2 , the method is still essentially a first-order method. While practical performance may be by far superior, the rates of convergence are at best the same as those of the pure first-order counterparts. In tang , an inexact PQNA with global sublinear rate of $\mathcal{O}({1}/{k})$ is proposed. While the algorithm can use any positive definite Hessian estimates, as long as their eigenvalues are uniformly bounded above and away from zero, the practical implementation proposed in tang used a limited memory BFGS Hessian approximation. The inexact setting of the algorithm allows for a relaxed sufficient decrease condition as well as an inexact subproblem optimization, for example via coordinate descent.

In this work, we show that PQNA, as proposed in tang , using general Hessian estimates $H_{k}$, has the linear convergence rate in the case of strongly convex smooth term $f$. Moreover, we consider an inexact variant, similar to the ones in nocedal3 ; tang , allowing inexact subproblem solutions as well as a relaxed sufficient decrease condition. In order to control the errors in the inexact subproblem optimization, we establish a simple stopping criterion for the subproblem solver, based on the iteration count, which guarantees that the inner subproblem is solved to the required accuracy. In contrast, in related works toh ; villa , it is assumed that an approximate subproblem solver yields an approximate subdifferential, which is a strong assumption on the subproblem solver which also does not clearly result in a simple stopping criterion.

Next, we apply Nesterov’s acceleration scheme to PQNA as proposed in tang , with a view of developing a version of this algorithm with faster convergence rates in the general convex case. In toh , the authors have introduced the Accelerated Proximal Quasi-Newton Algorithm (APQNA) with rate of convergence of $\mathcal{O}({1}/{k^{2}})$. However, this rate of convergence is achieved under condition $0\prec H_{k}\preceq H_{k-1}$, on the Hessian estimate $H_{k}$, at each iteration $k$. At the same time, this sequence of matrices has to be chosen so that $H_{k}$ is sufficiently positive definite to provide an overapproximation of $f$. Hence, these two conditions may contradict with each other unless the sequence of $\{H_{k}\}$ consists of unnecessarily large matrices. Moreover, in a particular case, when $H_{k}$ is set to be a scalar multiple of the identity, i.e., $H_{k}=\frac{1}{\mu_{k}}I$, then assumption $0\prec H_{k}\preceq H_{k-1}$ enforces $\mu_{k}\geq\mu_{k-1}$, implying nondecreasing step-size parameters, which contradicts the standard condition of APGA, which is $\mu_{k}\leq\mu_{k-1}$.

In this work, we introduce a new variant of APQNA, where we relax the restrictive assumptions imposed in toh . We use the scheme, originally introduced in goldfarb , which allows for the increasing and decreasing step-size parameters. We show that our version of APQNA achieves the convergence rate of $\mathcal{O}(1/k^{2})$ under some assumptions on the Hessian estimates. While we show that this assumption is rather strong and may not be satisfied by general matrices, it is not contradictory. Firstly, our result applies under the same restrictive condition from toh . We also show that our condition on the matrices holds in the case when the approximate Hessian at each iteration is a scaled version of the same “fixed” matrix $H$, which is a generalization of APGA. We investigate the performance of this algorithm in practice, and discover that this restricted version of a proximal quasi-Newton method is quite effective in practice. We also demonstrate that the general L-BFGS based PQNA does not benefit from the acceleration, which supports our analysis of the theoretical limitations.

This paper is organized as follows. In the next section, we describe the basic definitions and existing algorithms, PGA, APGA and PQNA, that we refer to later in the paper. In Section 3, we analyze PQNA in the strongly convex case. In Section 4, we propose, state and analyze a general APQNA and its simplified version. We present computational results in Section 5. Finally, we state our conclusions in Section 6.

Notation and Preliminaries

The proximal mapping of a convex function $g$ at a given point $v$, with parameter $\mu$ is defined as

The proximal mapping is the base operation of the proximal methods. In order to solve the composite problem (1), each iteration of PGA computes the proximal mapping of the function $g$ at point $x_{k}-\mu_{k}\nabla f(x_{k})$ as follows:

We will call the objective function, that is minimized in (3), a composite quadratic approximation of the convex function $F(x):=f(x)+g(x)$. This approximation at a given point $v$, for a given $\mu$ is defined as

Then, the proximal operator can be written as

Using this notation we first present the basic PGA framework with backtracking over $\mu$ in Algorithm 1. The simple backtracking scheme enforces that the sufficient decrease condition

holds. This condition is essential in the convergence rate analysis of PGA and is easily satisfied when $\mu\leq{1}/{L}$. The backtracking is used for two reasons–because the constant $L$ may not be known and because $\mu\leq{1}/{L}$ may be too pessimistic, i.e., condition (5) may be satisfied for much larger values of $\mu$ allowing for larger steps.

We now present the accelerated variant of PGA stated as APGA, where at each iteration $k$, instead of constructing $Q_{\mu_{k}}$ at the current minimizer $x_{k}$, it is constructed at a different point $y_{k}$, which is chosen as a specific linear combination of the latest two or more minimizers, e.g.

where the sequence $\{\alpha_{k}\}$ is chosen in such a way to guarantee an accelerated convergence rate as compared to the original PGA. Algorithm 2 is a variant of APGA framework, often referred to as FISTA, presented in beck , where $\alpha_{k}=({t_{k}-1})/({t_{k+1}})$. In this work, we choose to focus on FISTA algorithm specifically. The choice of the accelerated parameter $t_{k+1}$ in (6a) is dictated by the analysis of the complexity of FISTA beck and the condition $\mu_{k+1}\leq\mu_{k}$ that is imposed by the initialization of the backtracking procedure with $\mu_{k+1}^{0}:=\mu_{k}$. In goldfarb , the definition of $t_{k+1}$ was generalized to allow $\mu_{k+1}^{0}>\mu_{k}$, while retaining the convergence rate. We will use a similar technique in our proposed APQNA.

In this work, we are interested in the extensions of the above proximal methods, which utilize an approximation function $Q_{\mu}$, using partial second-order information about $f$. These quasi-Newton type proximal algorithms use a generalized form of the proximal operator (2), known as the scaled proximal mapping of $g$, which are defined for a given point $v$ as

where matrix $H$ is a positive definite matrix. In particular, the following operator

computes the minimizer, over $u$, of the following composite quadratic approximation of function $F$

when $v=x_{k}$. Matrix $H$ is the approximate Hessian of $f$ and its choice plays the key role in the quality of this approximation. Clearly, when $H=\frac{1}{\mu}I$, approximation (8) converts to (4), which is used throughout PGA. If we set $H=\nabla^{2}f(v)$, then (8) is the second-order Taylor approximation of $F$. At each iteration of PQNA the following optimization problem needs to be solved

where $z:=v-H^{-1}\nabla f(v)$. Obtaining such an inexact solution can be achieved by applying any linearly convergent algorithm, as will be discussing in detail at the end of this section.

In addition, for a given $\eta\in(0,1]$, the typical condition (5), used in beck and bach , is relaxed by using a trust-region like sufficient decrease condition

This relaxed condition was proposed and tested in tang for PQNA and was shown to lead to superior numerical performance, saving multiple backtracking steps during the earlier iterations of the algorithm. Note that, one can obtain the exact version of Algorithm 3 by replacing $p_{H_{k},\epsilon_{k}}$ with $p_{H_{k}}$, and setting $\eta=1$.

Throughout our analysis, we make the following assumptions.

$f$ is convex with Lipschitz continuous gradient with constant $L$.

$g$ is a lower semi-continuous proper convex function.

There exists positive constants $m$ and $M$ such that, for all $k>0$,

In Algorithm 3, as long as the sequence of positive definite matrices $G_{k}$ has uniformly bounded eigenvalues, condition (12) is satisfied. In fact, since the sufficient decrease condition in Step 3 is satisfied for $H_{k}\succeq LI$, then it is satisfied when $\mu_{k}\leq{1}/{L}$. Hence, at each iteration we have a finite and bounded number of backtracking steps and the resulting $H_{k}$ has bounded eigenvalues. The lower bound on the eigenvalues of $H_{k}$ is simply imposed either by choosing a positive definite $G_{k}$ or bounding $\mu_{k}^{0}$ from above.

In the next section, we analyze the convergence properties of PQNA when $f$ in (1) is strongly convex.

Analysis of the Proximal Quasi-Newton Algorithm under Strong Convexity

In this section, we analyze the convergence properties of PQNA to solve problem (1), in the case when the smooth function $f$ is $\gamma-$strongly convex. In particular, the following assumption is made throughout this section.

To establish a linear convergence rate of PQNA, we consider extending two different approaches used to show a similar result for PGA. The first approach we consider can be found in phd , and is based on the proof techniques used in dmitri for PGA. The reason we chose the approach in dmitri is due to the fact that the linear rate of convergence is shown under the quadratic growth condition, which is a relaxation of the strong convexity. Hence, extending this analysis to PQNA, as a subject of a future work, may allow us to relax the strong convexity assumption for this algorithm as well. However, there appears to be some limitations in the extension of this analysis phd , in particular in the inexact case. This observation motivates us to present the approach used in nesterov4 to analyze convergence properties of inexact PQNA. As we see below, this analysis readily extends to our case and allows us to establish simple rules for subproblem solver termination to achieve the desired subproblem accuracy.

Let us consider Algorithm 3 for which (10) holds for some sequence of errors $\epsilon_{k}\geq 0$. The relaxed sufficient decrease condition

for a given $\eta\in(0,1]$, can be written as

where the sequence of the errors $\xi_{k}$ is defined as

In particular, setting $\eta=1$ results in $\xi_{k}=0$, for all $k$ and enforces the algorithm to accept only those steps that achieve full (predicted) reduction. However, using $\eta<1$ allows the algorithm to take steps satisfying only a fraction of the predicted reduction, which may lead to larger steps and faster progress.

Under the above inexact condition, we can show the following result.

Suppose that Assumptions 2.1 and 3.1 hold. At each iteration of the inexact PQNA, stated in Algorithm 3, we have

when $\rho=1-{(\eta\gamma)}/{(\gamma+M)}$, and

Applying (14), with $v=x_{k}$ and consequently $p_{H_{k},\epsilon_{k}}(x_{k})=x_{k+1}$, we have

Now, by substituting the expression for $\xi_{k}$, as stated in (15), we will have

where $\rho=1-\eta t^{\prime}$. Now, we can conclude the final result as

where $\rho=1-{(\eta\gamma)}/{(\gamma+M)}$. ∎

In Theorem 3.2, by setting $\epsilon_{k}=0$ and $\eta=1$, which implies $\xi_{k}=0$, we achieve the linear convergence rate of the exact variant of PQNA.

We have shown that the linear rate of PQNA is $\rho=1-{(\eta\gamma)}/{(\gamma+M)}$. As argued in Remark 1, $M$ is of the same order as $L$ in the worst case, hence in that case the linear rate of PQNA is the same as that of the simple PGA. However, it is easy to see that in the proof of Theorem 3.2, the linear rate is derived using the upper bound on $\|x_{*}-x_{k}\|_{H_{k}}^{2}$, where $H_{k}$ is the approximate Hessian on step $k$. Clearly, the idea of using the partial second-order information is to reduce the worst case bound of $H_{k}$ in general and consequently on $\|x_{*}-x_{k}\|_{H_{k}}^{2}$. In particular, obtaining a smaller bound $M_{k}$ on each iteration yields a larger convergence coefficient $\rho_{k}=1-{(\eta\gamma)}/{(\gamma+M_{k})}$. While for general $H_{k}$, we do not expect to improve upon the regular PGA in theory, this remark serves to explain the better performance of PQNA in practice.

where $\alpha\in(0,1)$. Our goal is to ensure that $\epsilon_{k}\leq C\rho^{k\cdot\delta}$, which can be achieved by applying sufficient number of iterations of the subproblem algorithm. To be specific, the following theorem characterizes this required bound on the number of inner iterations.

Suppose that at the $k$-th iteration of Algorithm 3, after applying the subproblem solver satisfying (17) for $r(k)$ iterations, starting with $u_{0}=x_{k}$, we obtain solution $x_{k+1}=u_{r(k)}$. Let $r(k)$ satisfy

for some $\delta>1$, and $\rho$ defined in Theorem 3.2. Then

holds for all $k$, with $C$ being the uniform bound on $Q_{H_{k}}(x_{k},x_{k})-Q_{k}^{*}$, and the linear convergence of Algorithm 3 is achieved.

First, assume that at the $k-$th iteration we have applied the subproblem solver for ${r(k)}$ iterations to minimize strongly convex function $Q_{H_{k}}$. Now, by combining $Q_{H_{k}}(u_{0},u_{0})-Q_{k}^{*}\leq C$ and (17), we can conclude the following upper bound, so that

Now, if $\alpha^{r(k)}C\leq\epsilon_{k}$, we can guarantee that $u_{r(k)}$ is an $\epsilon_{k}$-solution of the $k$-th subproblem, so that $Q_{H_{k}}(u_{r(k)},u_{0})\leq Q_{k}^{*}+\epsilon_{k}$. Now, assuming that $\rho$ is known, we can set the error rate of the $k$-th iteration as $\epsilon_{k}\leq C\rho^{k\cdot\delta}$, for a fixed $\delta>1$. In this case, the number of inner iterations which guarantees the $\epsilon_{k}$-minimizer will be

Since subproblems are strongly convex, the required linear convergence rate for the subproblem solver, stated in (17), can be guaranteed via some basic first-order algorithms or their accelerated variants. However, one difficulty in obtaining lower bound (18) is that it depends on the prior knowledge of $\rho$ and $\alpha$. Consider the following simple modification of Theorem 3.3; instead of condition (18), consider $r(k)$ satisfying

In the next subsection, we extend our analysis to the case of solving subproblems via the randomized coordinate descent, where at each iteration the desired error bound related to $\epsilon_{k}$ is only satisfied in expectation.

2 Solving Subproblems via Randomized Coordinate Descent

As we mentioned before, in order to achieve linear convergence rate of the inexact PQNA, any simple first-order method (such as PGA) can be applied. However, as discussed in tang , in the case when $g(x)=\lambda\|x\|_{1}$ and $H_{k}$ is sum of a diagonal and a low rank matrix, as in the case of L-BFGS approximations, the coordinate descent method is the most efficient approach to solve the strongly convex quadratic subproblems. In this case, each iteration of coordinate descent has complexity of $\mathcal{O}(m)$, where $m$ is the memory size of L-BFGS, which is usually chosen to be less than $20$, while each iteration of a proximal gradient method has complexity of $\mathcal{O}(nm)$ and each iteration of the Newton type proximal method has complexity of $\mathcal{O}(nm^{2})$. While more iterations of coordinate descent may be required to achieve the same accuracy, it tends to be the most efficient approach. To extend our theory of the previous section and to establish the bound on the number of coordinate descent steps needed to solve each subproblem, we utilize convergence results for the randomized coordinate descent martin , as is done in tang .

Algorithm 4 shows the framework of the randomized coordinate descent method, which can be used as a subproblem solver of Algorithm 3 and is identical to the method used in tang . In Algorithm 4, function $Q_{H}$ is iteratively minimized over a randomly chosen coordinate, while the other coordinates remain fixed.

In what follows, we restate Theorem 6 in martin , which establishes linear convergence rate of the randomized coordinate descent algorithm, in expectation, to solve strongly convex problems.

Suppose we apply randomized coordinate descent for $r$ iterations, to minimize the $m$-strongly convex function $Q$ with $M$-Lipschitz gradient, to obtain the random point $u_{r}$.

Now, we want to analyze how we can utilize the result of Theorem 3.4 to achieve the linear convergence rate of inexact PQNA, in expectation. Toward this end, first we need the following theorem as the probabilistic extension of Theorem 3.2.

Suppose that Assumptions 2.1 and 3.1 hold. At each iteration $k$ of the inexact PQNA, stated in Algorithm 3, assume that the error $\epsilon_{k}$ is a nonnegaitve random variable defined on some probability space with an arbitrary distribution. Then, we have

when $\rho=1-{(\eta\gamma)}/{(\gamma+M)}$, and

The proof is a trivial modification of that of Theorem 3.2. ∎

In what follows, we describe how the randomized coordinate descent method ensures the required accuracy of subproblems and consequently guarantees linear convergence of the inexact PQNA.

Suppose that at the $k$-th iteration of Algorithm 3, after applying Algorithm 4 for $r(k)$ iterations, starting with $u_{0}=x_{k}$, we obtain solution $x_{k+1}=u_{r(k)}$. If

Accelerated Proximal Quasi-Newton Algorithm

We now turn to an accelerated variant of PQNA. As we described in the introduction section, the algorithm proposed in toh is a proximal quasi-Newton variant of FISTA, described in Algorithm 2. In toh , the convergence rate of $\mathcal{O}({1}/{k^{2}})$ is shown under the condition that the Hessian estimates satisfy $0\prec H_{k}\preceq H_{k-1}$, at each iteration. On the other hand, the sequence $\{H_{k}\}$ is chosen so that the quadratic approximation of $f$ is an over approximation. This leads to an unrealistic setting where two possible contradictory conditions need to be satisfied and as mentioned earlier, this condition contradicts the assumptions of the original APGA, stated in Algorithm 2. We propose a more general version, henceforth referred to as APQNA, which allows a more general sequence of $H_{k}$ and is based on the relaxed version of FISTA, proposed in goldfarb , which does not impose monotonicity of the step-size parameters. Moreover, our algorithm allows more general Hessian estimates as we explain below.

The main framework of APQNA as stated in Algorithm 5 is similar to that of Algorithm 2, where the simple composite quadratic approximation $Q_{\mu}$ was replaced by the scaled version $Q_{H}$, as is done in Algorithm 3, using (partial) Hessian information. As in the case of Algorithm 3, we assume that the approximate Hessian $H_{k}$ is a positive definite matrix such that $mI\preceq H_{k}\preceq MI$, for some positive constants $m$ and $M$. As discussed in Remark 1, it is simple to show that this condition can be satisfied for any positive $m$ and for any large enough $M$. Here, however, we will need additional much stronger assumptions on the sequence $\{H_{k}\}$. The algorithm, thus, has some additional steps compared to Algorithm 2 and the standard FISTA type proximal quasi-Newton algorithm proposed in toh . Below, we present Algorithm 5 and discuss the steps of each iteration in detail.

The key requirement imposed by Algorithm 5 on the sequence $\{H_{k}\}$ is that $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$, while $\theta_{k}:={\sigma_{k}}/{\sigma_{k+1}}$ is used to evaluate the accelerated parameter $t_{k+1}$ through (24a). During Steps 4 and 5 of iteration $k$, initial guesses for $\sigma_{k+1}^{0}$ and $H_{k+1}$ are computed and used to define $\theta_{k}$, which is then used to compute $t_{k+1}$ and $y_{k+1}$. Since the approximate Hessian $H_{k+1}$ may change during Step 2 of iteration $k+1$, $\sigma_{k+1}$ may need to change as well in order to satisfy condition $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$. In particular, we may shrink the value of $\sigma_{k+1}$ and consequently will need to recompute $\theta_{k}$ and, thus, $t_{k+1}$ and $y_{k+1}$. Therefore, the backtracking process in Step 2 of Algorithm 5 involves a loop which may require repeated computations of $y_{k}$ and hence $\nabla f(y_{k})$.

We do not specify how to compute $H_{k}$ in Algorithm 5, as long as it satisfies (12) and condition $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$. Note that Algorithm 5 does not allow the use of exact Hessian information at $y_{k+1}$, i.e., $H_{k+1}=\nabla^{2}f(y_{k+1})$, because it is assumed that $H_{k+1}$ is computed before $y_{k+1}$ (since $y_{k+1}$ uses the value of $\sigma_{k+1}$, whose value may have to be dependent on $H_{k+1}$). However, it is possible to use $H_{k+1}=\nabla^{2}f(x_{k})$ in Algorithm 5. To use $H_{k+1}=\nabla^{2}f(y_{k+1})$, one would need to be able to compute $\sigma_{k+1}$ before $H_{k+1}$ and somehow ensure that condition $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$ is satisfied. This condition can eventually be satisfied by applying similar technique to Step 2, but in that case $H_{k+1}$ will not be equal to the Hessian, but to some multiple of the Hessian, i.e., $\frac{1}{\beta^{i}}\nabla^{2}f(y_{k+1})$, for some $i$.

In our numerical results, we construct $H_{k}$ via L-BFGS and ignore condition $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$, since enforcing it in this case causes a very rapid decrease in $\sigma$. It is unclear, however, if a practical version of Algorithm 5, based on L-BFGS Hessian approximation can be derived, which may explain why the accelerated version of our algorithm does not represent any significant advantage.

One trivial choice of the matrix sequence is $H_{k}=\frac{1}{\mu_{k}}I$. In this case, the sequence of scalars $\sigma_{k}=\mu_{k}$, satisfies $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$, for all $k$. This choice of Hessian reduces Algorithm 5 to the version of APGA with full backtracking of the step-size parameters, proposed in goldfarb , hence Algorithm 5 is the generalization of that algorithm. Another choice for the matrix sequence is $H_{k}=\frac{1}{\sigma_{k}}H$, where the matrix $H$ is any fixed positive definite matrix. This setting of $H_{k}$ automatically satisfies condition $\sigma_{k+1}H_{k+1}\preceq\sigma_{k}H_{k}$, and Algorithm 5 reduces to the simplified version stated below in Algorithm 6.

Note that, by the same logic that was used in Remark 1, the number of backtracking steps at each iteration of Algorithm 6 is uniformly bounded. Thus, as long as the fixed approximate Hessian $H$ is positive definite, a Hessian estimate $H_{k}=\frac{1}{\sigma_{k}}H$ has positive eigenvalues bounded from above and below. In our implementation, we compute a fixed matrix $H$ by applying L-BFGS for a fixed number of iterations and then apply Algorithm 6.

In the next section, we analyze the convergence properties of Algorithm 5, where the approximate Hessian $H_{k}$ is produced by some generic unspecified scheme. The motivation is to be able to apply the analysis to popular and efficient Hessian approximation methods, such as L-BFGS. However, in the worst case for general $H_{k}$, a positive lower bound for $\{\sigma_{k}\}$ can not be guaranteed for such a generic scheme. This observation motivates the analysis of Algorithm 6, as a simplified version of Algorithm 5. It remains to be seen if some bound on $\{\sigma_{k}\}$ may be derived for matrices arising specifically via L-BFGS updates.

2 Convergence Analysis

In this section, we prove that if the sequence $\{\sigma_{k}\}$ is bounded away from zero, Algorithm 5 achieves the same rate of convergence as APGA, i.e., $\mathcal{O}({1}/{k^{2}})$. First, we state a simple result based on the optimality of $p_{H}$.

The proof is followed immediately from the optimality condition of the convex optimization problem (9). ∎

Now, we can show the following lemma, which bounds the change in the objective function $F$ and is a simple extension of Lemma 2.3 in beck .

Now, based on the convexity of functions $f$ and $g$, we have

where $\nu_{g}(p_{H}(v))$ is defined in Lemma 1. Summing the above inequalities yields

The following result is a simple corollary of Lemma 2.

The result immediately follows by applying the following identity

The next lemma states the key properties which are used in the convergence analysis.

At each iteration of Algorithm 5, the following relations hold

The proof follows trivially from the conditions in Algorithm 5 and the fact that $\theta_{k}\leq{\sigma_{k}}/{\sigma_{k+1}}$. ∎

Now, using this lemma and previous results we derive the key property of the iterations of APQNA.

For all $k\geq 1$ for Algorithm 5 we have

where $v_{k}=F(x_{k})-F(x_{*})$ and $u_{k}=t_{k}x_{k}-(t_{k}-1)x_{k-1}-x_{*}$.

In (29), by setting $v=y_{k+1}$, $p_{H}(v)=x_{k+1}$, $H=H_{k+1}$, and $x=x_{k}$ and then by multiplying the resulting inequality by $\sigma_{k+1}(t_{k+1}-1)$, we will have

On the other hand, in (29), by setting $x=x_{*}$ and multiplying it by $\sigma_{k+1}$, we have

By adding these two inequalities, we have

Multiplying the last inequality by $t_{k+1}$ and applying inequality (31b) give

Hence, by using the definition of $y_{k+1}$ and $u_{k}$, we have

Now, we are ready to state and prove the convergence rate result.

The sequence of iterates $x_{k}$, generated by Algorithm 5, satisfies

By setting $t_{1}=1$, using the definition of $u_{k}$ at $k=1$, which is $u_{1}=x_{1}-x_{*}$, and also considering the positive definiteness of $H_{k}$ for all $k\geq 1$, it follows from Lemma 4 that

Setting $x=x_{*}$, $v=y_{1}=x_{0}$, $p_{H}(v)=x_{1}$, $t_{1}=1$, and $H=H_{1}$ in (29) implies

Multiplying the above by $\sigma_{1}$ gives

Finally, by setting $\sigma_{1}=1$ and $H_{1}=I$, we obtain

Now, based on the result of Theorem 4.1, in order to obtain the rate of convergence of $\mathcal{O}({1}/{k^{2}})$ for Algorithm 5, it is sufficient to show that

for some constant $\psi>0$. The next result is a simple consequence of the relation (31b), or equivalently (24a).

The sequence $\{\sigma_{k}\}$ generated by Algorithm 5 satisfies

We can prove this lemma by using induction. Trivially, for $k=1$, since $t_{1}=1$, the inequality holds. As the induction assumption, assume that for $k>1$, we have $\sigma_{k}t_{k}^{2}\geq\left(\frac{\sum_{i=1}^{k}\sqrt{\sigma_{i}}}{2}\right)^{2}$. Since (24a) holds for all $k$, it follows that

Multiplying by $\sqrt{\sigma_{k+1}}$ implies

Finally, by using induction assumption, we will have have

Hence, if we assume that the sequence $\{\sigma_{k}\}$ is bounded below by a positive constant $\underline{\sigma}$, i.e., $\sigma_{k}\geq\underline{\sigma}$, we can establish the desired bound on $\sigma_{k}t_{k}^{2}$, as stated in the following theorem.

If for all iterations of Algorithm 5 we have $\sigma_{k}\geq\underline{\sigma}$, then for all $k$

Under the assumption $\sigma_{k}\geq\underline{\sigma}$, we will have

and consequently, by using Lemma 5, we obtain

Then, by using Theorem 4.1, we have the desired rate of convergence of $\mathcal{O}({1}/{k^{2}})$ as stated in (33). ∎

The assumption of the existence of a bounded sequence $\{\sigma_{k}\}$ such that $\sigma_{k}\geq\underline{\sigma}$ and (31b) holds may not be satisfied when we use a general approximate Hessian. To illustrate this, consider the following simple sequence of matrices:

Clearly, $\sigma_{2k+1}\leq\sigma_{2k}/10$ and $\sigma_{2k}\leq\sigma_{2k-1}/10$, and hence $\sigma_{k}\leq 10^{-k}$. In this case, based on the result of Theorem 4.1, we cannot guarantee any convergence result. Some convergence result can still be attained, when $\sigma_{k}\to 0$, for example, if $\sigma_{k}\geq{\underline{\sigma}}/{k}$, as we show in the following relaxed version of Theorem 4.2.

If for all iterations of Algorithm 5 we have $\sigma_{k}\geq{\underline{\sigma}}/{k}$, then for all $k$

From $\sigma_{k}\geq{\underline{\sigma}}/{k}$, we will have

and consequently, by using Lemma 5, we obtain

Then, by using Theorem 4.1, we have (34). ∎

The above theorem shows that if $\sigma_{k}$ converges to zero, but not faster than $1/k$, then our APQNA method may loose its accelerated rate of convergence, but still converges at least at the same rate as PQNA. Establishing lower bounds of $\sigma_{k}$ for different choices of Hessian estimates is a nontrivial task and is the subject of future research. As we will demonstrate in our computational section, APQNA with L-BFGS Hessian approximation does not seem to have any practical advantage over its nonaccelerated counterpart, however it is clearly convergent.

We can establish the accelerated rate of Algorithm 6, since in this case we can guarantee a lower bound on $\sigma_{k}$, due to the restricted nature of the $H_{k}$ matrices.

In Algorithm 6, let $mI\preceq H$, then $\sigma_{k}\geq{\beta m}/{L}$ and hence the convergence rate of $\mathcal{O}({1}/{k^{2}})$ is achieved.

In Algorithm 6, we define $H_{k}=\frac{1}{\sigma_{k}}H$. The sufficient decrease condition $F(p_{H_{k}}(y_{k}))\leq Q_{H_{k}}(p_{H_{k}}(y_{k}),y_{k})$, is satisfied for any $H_{k}\succeq LI$, hence it is satisfied for any $H_{k}=\frac{1}{\sigma_{k}}H$ with $\sigma_{k}\leq{m}/{L}$. By the mechanism of Step 3 in Algorithm 6, we observe that for all $k$, we have $\sigma_{k}\geq{\beta m}/{L}$. Let us note now that Algorithm 6 is a special case of Algorithm 5, hence all the above results, in particular Theorem 4.1 and Lemma 5 hold. Consequently, based on Theorem 4.2, the desired convergence rate of $\mathcal{O}({1}/{k^{2}})$ for Algorithm 6 is obtained. ∎

We have studied only the exact variant of APQNA in this section. Incorporating inexact subproblem solutions, as was done for APQNA in the previous section, is relatively straightforward following the techniques for inexact APGA, mark . It is easy to show that if the exact algorithm has the accelerated convergence rate, then the inexact counterpart, with subproblems solved by a linearly convergent method, such as randomized coordinate descent, inherits this convergence rate. However, using the relaxed sufficient decrease condition does not apply here as it does not preserve the accelerated convergence rate.

In the next section, we present the numerical results comparing the performance of Algorithm 5 and Algorithm 6 to their nonaccelerated counterparts, to see how much practical acceleration is achieved.

Numerical Experiments

In this section, we investigate the practical performance of several algorithms discussed in this work, applied to the sparse logistic regression problem

The algorithms that we compare here are as follows:

Accelerated Proximal Gradient Algorithm (APGA), proposed in beck , (also known as FISTA),

Proximal Quasi-Newton Algorithm with Fixed Hessian approximation (call it PQNA-FH),

Accelerated Proximal Quasi-Newton Algorithm with Fixed Hessian approximation (call it APQNA-FH),

Proximal Quasi-Newton Algorithm with L-BFGS Hessian approximation (call it PQNA-LBFGS), proposed in tang , and

Accelerated Proximal Quasi-Newton Algorithm with L-BFGS Hessian approximation (call it APQNA-LBFGS).

In PQNA-FH and APQNA-FH, we set $H_{k}=\frac{1}{\sigma_{k}}H$, where $H$ is a positive definite matrix computed via applying L-BFGS updates over the first few iterations of the algorithm which then is fixed for all remaining iterations. On the other hand, PQNA-LBFGS and APQNA-LBFGS employ the L-BFGS updates to compute Hessian estimates throughout the algorithm. In all of the above algorithms, we use the coordinate descent scheme, as described in tang , to solve the subproblems inexactly. According to the theory in tang , PQNA-FH and PQNA-LBFGS converge at the rate of $\mathcal{O}({1}/{k})$. If $f$ is strongly convex (which depends on the problem data), then according to Theorem 3.2,PQNA-FH and PQNA-LBFGS converge at a linear rate. By Lemma 6, in APQNA-FH, condition $\sigma_{k}H_{k}\preceq\sigma_{k-1}H_{k-1}$ holds automatically and the algorithm converges at the rate of $\mathcal{O}({1}/{k^{2}})$. On the other hand, for APQNA-LBFGS, condition $\sigma_{k}H_{k}\preceq\sigma_{k-1}H_{k-1}$ has to be enforced. We have tested various implementations that ensure this condition and none have produced a practical approach. We then chose to set $\theta_{k}=1$ and relax the condition $\sigma_{k}H_{k}\preceq\sigma_{k-1}H_{k-1}$. The resulting algorithm is practical and is empirically convergent, but as we will see does not provide an improvement over PQNA-LBFGS.

The algorithms are implemented in MATLAB R2014b and computations were performed on the COR@L computational cluster of the ISE department at Lehigh, consisting of 16-cores AMD Operation, 2.0 GHz nodes with 32 Gb of memory.

First, in order to demonstrate the effect of using even limited Hessian information within an accelerated method, we compared the performance of APQNA-FH and APGA, both in terms of the number of iterations and the total solution time, see the results in Table 2.

Based on the results shown in Table 2, we conclude that APQNA-FH consistently dominates the APGA, both in terms of the number of function evaluations and also in terms of the total solution time.

It is worth mentioning that although in terms of computational effort, each iteration of APGA is cheaper than each iteration of APQNA-FH, the total solution time of APQNA-FH is significantly less than APGA, due to the smaller number of iterations of APQNA-FH compared to APGA.

The next experiment is to compare the performance of APQNA-FH and PQNA-FH to observe the effect of acceleration in the fixed matrix setting. This comparison is done in terms of the number of iterations and the number of function evaluations, and is shown in Figure 1 and Figure 2, respectively. The subproblem solution time is the same for both algorithms. As we can see in Figure 1, in terms of the number of iterations, APQNA-FH dominates PQNA-FH, for a range of memory sizes of L-BFGS which have been used to compute matrix $H$. Moreover, as is seen in Figure 2, APQNA-FH dominates PQNA-FH, in terms of the number of function evaluations, even though each iteration of APQNA-FH requires two function evaluations, because of the nature of the accelerated scheme. This shows that APQNA-FH achieves practical acceleration compared to PQNA-FH, as supported by the theory in the previous section.

Next, we compare the performance of APQNA-FH versus APQNA-LBFGS to compare the effect of using the fixed approximate Hessian $H_{k}=\frac{1}{\sigma_{k}}H$, which satisfies condition $\sigma_{k}H_{k}\preceq\sigma_{k-1}H_{k-1}$ versus using variable Hessian estimates computed via L-BFGS method at each iteration, while relaxing condition $\sigma_{k}H_{k}\preceq\sigma_{k-1}H_{k-1}$. Table 3 shows the results of this comparison, obtained based on the best choices of memory size for L-BFGS, in particular $\bar{k}=8$ and $\bar{k}=9$, respectively. As we can see, these two algorithms are competitive both in terms of the number of iterations and also the total solution time. Since APQNA-FH does not use the local information of function $f$ to approximate $H_{k}$, it often takes more iterations than APQNA-LBFGS, which constantly updates $H_{k}$ matrices. On the other hand, since APQNA-FH does not require additional computational effort to evaluate $H_{k}$, hence one iteration of this algorithm is cheaper than one iteration of APQNA-LBFGS, which causes the competitive total solution time.

Finally, we compare APQNA-LBFGS and PQNA-LBFGS, to demonstrate the effect of using an accelerated scheme in the quasi-Newton type proximal algorithms. The results of this comparison are shown in Figure 3 and Figure 4 in terms of the number of iterations and the number of function evaluations, respectively, for different memory sizes of L-BFGS Hessian approximation.

Clearly, as we can see in Figure 3 and 4, not only does the accelerated scheme not achieve practical acceleration compared to PQNA-LBFGS in terms of the number of iterations, but it is also inferior in terms of the number of function evaluations, since every iteration requires two function evaluations. Thus, we believe that the practical experiments support our theoretical analysis in that applying acceleration scheme in the case of variable Hessian estimates may not result in a faster algorithm.

Conclusion

In this paper, we established a linear convergence rate of PQNA proposed in tang under a strong convexity assumption. To our knowledge, this is the first such result, for proximal quasi-Newton type methods, which have lately been popular in the literature. We also show that this convergence rate is preserved when subproblems are solved inexactly. We provide a simple and practical rule for the number of inner iterations that guarantee sufficient accuracy of subproblem solutions. Moreover, we allow a relaxed sufficient decrease condition during backtracking, which preserves the convergence rate, while it is known to improve the practical performance of the algorithm.

Furthermore, we presented a variant of APQNA as an extension of PQNA. We have shown that this algorithm has the convergence rate of $\mathcal{O}({1}/{k^{2}})$ under a strong condition on the Hessian estimates, which can not always be guaranteed in practice. We have shown that this condition holds when Hessian estimates are a multiple of a fixed matrix, which is computationally less expensive than the more common methods, such as the L-BFGS scheme. Although, this proposed algorithm has the same rate of convergence as the classic APGA, it is significantly faster in terms of the final number of iterations and also the total solution time. Based on the theory, using L-BFGS Hessian approximation, may result in worse convergence rate, however, our computational results show that the practical performance is about the same as that while the fixed matrix. On the other hand, although in these two algorithms, we are applying the accelerated scheme, their practical performances are inferior to that of PQNA-LBFGS, which does not use any accelerated scheme and potentially has a slower sublinear rate of convergence in the absence of strong convexity. We conclude that using variable Hessian estimates is the most efficient approach, and will result in the linear convergence rate in the presence of strong convexity, but that a standard accelerated scheme is not useful in this setting. Exploring other possibly more effective accelerated schemes for the proximal quasi-Newton methods is the subject of future research.

References