Accelerating Greedy Coordinate Descent Methods

Haihao Lu, Robert M. Freund, Vahab Mirrokni

Introduction: Related Work and Accelerated Coordinate Descent Framework

Coordinate descent methods have received much-deserved attention recently due to their capability for solving large-scale optimization problems (with sparsity) that arise in machine learning applications and elsewhere. With inexpensive updates at each iteration, coordinate descent algorithms obtain faster running times than similar full gradient descent algorithms in order to reach the same near-optimality tolerance; indeed some of these algorithms now comprise the state-of-the-art in machine learning algorithms for loss minimization.

Most recent research on coordinate descent has focused on versions of randomized coordinate descent, which can essentially recover the same results (in expectation) as full gradient descent, including obtaining “accelerated” (i.e., O(1/k2)O(1/k^{2})) convergence guarantees. On the other hand, in some important machine learning applications, greedy coordinate methods demonstrate superior numerical performance while also delivering much sparser solutions. For example, greedy coordinate descent is one of the fastest algorithms for the graphical LASSO implemented in DP-GLASSO . And sequence minimization optimization (SMO) (a variant of greedy coordinate descent) is widely known as the best solver for kernel SVM and is implemented in LIBSVM and SVMLight.

In general, for smooth convex optimization the standard first-order methods converge at a rate of O(1/k)O(1/k) (including greedy coordinate descent). In 1983, Nesterov proposed an algorithm that achieved a rate of O(1/k2)O(1/k^{2}) – which can be shown to be the optimal rate achievable by any first-order method . This method (and other similar methods) is now referred to as Accelerated Gradient Descent (AGD).

However, there has not been much work on accelerating the standard Greedy Coordinate Descent (GCD) due to the inherent difficulty in demonstrating O(1/k2)O(1/k^{2}) computational guarantees (we discuss this difficulty further in Section 4.1). One work that might be close is , which updates the zz-sequence using the full gradient and thus should not be considered as a coordinate descent method in the standard sense. There is a very related concurrent work and we will discuss the connections to our results in Section 4.1.

In this paper, we study ways to accelerate greedy coordinate descent in theory and in practice. We introduce and study two algorithms: Accelerated Semi-Greedy Coordinate Descent (ASCD) and Accelerated Greedy Coordinate Descent (AGCD). While ASCD takes greedy steps in the xx-updates and randomized steps in the zz-updates, AGCD is a straightforward extension of GCD that only takes greedy steps. On the theory side, our main results are for ASCD: we show that ASCD achieves O(1/k2)O(1/k^{2}) convergence, and it also achieves accelerated linear convergence when the objective function is furthermore strongly convex. However, a direct extension of convergence proofs for ARCD does not work for ASCD due to the different coordinates we use to update xx-sequence and zz-sequence. Thus, we present a new proof technique – which shows that a greedy coordinate step yields better objective function value than a full gradient step with a modified smoothness condition.

On the empirical side, we first note that in most of our experiments ASCD outperforms Accelerated Randomized Coordinate Descent (ARCD) in terms of running time. On the other hand, we note that AGCD significantly outperforms the other accelerated coordinate descent methods in all instances, in spite of a lack of theoretical guarantees for this method. To complement the empirical study of AGCD, we present a Lyapunov energy function argument that points to an explanation for why a direct extension of the proof for AGCD does not work. This argument inspires us to introduce a technical condition under which AGCD is guaranteed to converge at an accelerated rate. Interestingly, we confirm that technical condition holds in a variety of instances in our empirical study, which in turn justifies our empirical observation that AGCD works very well in practice.

Coordinate Descent. Coordinate descent methods have a long history in optimization, and convergence of these methods has been extensively studied in the optimization community in the 1980s-90s, see , , and . There are roughly three types of coordinate descent methods depending on how the coordinate is chosen: randomized coordinate descent (RCD), cyclic coordinate descent (CCD), and greedy coordinate descent (GCD). RCD has received much attention since the seminal paper of Nesterov . In RCD, the coordinate is chosen randomly from a certain fixed distribution. provides an excellent review of theoretical results for RCD. CCD chooses the coordinate in a cyclic order, see for basic convergence results. More recent results show that CCD is inferior to RCD in the worst case , while it is better than RCD in certain situations . In GCD, we select the coordinate yielding the largest reduction in the objective function value. GCD usually delivers better function values at each iteration in practice, though this comes at the expense of having to compute the full gradient in order to select the gradient coordinate with largest magnitude. The recent work shows that GCD has faster convergence than RCD in theory, and also provides several applications in machine learning where the full gradient can be computed cheaply. A parallel GCD method is proposed in and numerical results show its advantage in practice.

Accelerated Randomized Coordinate Descent. Since Nesterov’s paper on RCD there has been significant focus on accelerated versions of RCD. In particular, developed the first accelerated randomized coordinate gradient method for minimizing unconstrained smooth functions. present a sharper convergence analysis of Nesterov’s method using a randomized estimate sequence framework. proposed the APPROX (Accelerated, Parallel and PROXimal) coordinate descent method and obtained an accelerated sublinear convergence rate, and developed an efficient implementation of ARCD.

2 Accelerated Coordinate Descent Framework

Let ⟨x,y⟩L:=∑i=1nLixiyi=⟨x,Ly⟩=⟨Lx,y⟩\langle x,y\rangle_{L}:=\sum_{i=1}^{n}L_{i}x_{i}y_{i}=\langle x,{\bf{L}}y\rangle=\langle{\bf{L}}x,y\rangle denote the LL-inner product. Define the norm ∥x∥L:=⟨x,Lx⟩\|x\|_{L}:=\sqrt{\langle x,{\bf{L}}x\rangle}. Letting L−1{\bf{L}}^{-1} denote the inverse of L{\bf{L}}, we will also use the norm ∥⋅∥L−1\|\cdot\|_{L^{-1}} defined by ∥v∥L−1:=⟨v,L−1v⟩=∑i=1nLi−1vi2\|v\|_{L^{-1}}:=\sqrt{\langle v,{\bf{L}}^{-1}v\rangle}=\sqrt{\sum_{i=1}^{n}{\bf{L}}_{i}^{-1}v_{i}^{2}}.

Algorithm 1 presents a generic framework for accelerated coordinate descent methods that is flexible enough to encompass deterministic as well as randomized methods. One specific case is the standard Accelerated Randomized Coordinate Descent (ARCD). In this paper we propose and study two other cases. The first is Accelerated Greedy Coordinate Descent (AGCD), which is a straightforward extension of greedy coordinate descent to the acceleration framework and which, surprisingly, has not been previously studied (that we are aware of). The second is a new algorithm which we call Accelerated Semi-Greedy Coordinate Descent (ASCD) that takes greedy steps in the xx-updates and randomized steps in the zz-update.

In the framework of Algorithm 1 we choose a coordinate jk1j_{k}^{1} of the gradient ∇f(yk)\nabla f(y^{k}) to perform the update of the xx-sequence, and we choose (a possibly different) coordinate jk2j_{k}^{2} of the gradient ∇f(yk)\nabla f(y^{k}) to perform the update of the zz-sequence. Herein we will study three different rules for choosing the coordinates jk1j_{k}^{1} and jk2j_{k}^{2} which then define three different specific algorithms:

ARCD (Accelerated Randomized Coordinate Descent): use the rule

AGCD (Accelerated Greedy Coordinate Descent): use the rule

ASCD (Accelerated Semi-Greedy Coordinate Descent): use the rule

In ARCD a random coordinate jk1j_{k}^{1} is chosen at each iteration kk, and this coordinate is used to update both the xx-sequence and the zz-sequence. ARCD is well studied, and is known to have the following convergence guarantee in expectation (see for details):

where the expectation is on the random variables used to define the first kk iterations.

In AGCD we choose the coordinate jk1j_{k}^{1} in a “greedy” fashion, i.e., corresponding to the maximal (weighted) absolute value coordinate of the the gradient ∇f(yk)\nabla f(y^{k}). This greedy coordinate is used to update both the xx-sequence and the zz-sequence. As far as we know AGCD has not appeared in the first-order method literature. One reason for this is that while AGCD is the natural accelerated version of greedy coordinate descent, the standard proof methodologies for establishing acceleration guarantees (i.e., O(1/k2)O(1/k^{2}) convergence) fail for AGCD. Despite this lack of worst-case guarantee, we show in Section 5 that AGCD is extremely effective in numerical experiments on synthetic linear regression problems as well as on practical logistic regression problems, and dominates other coordinate descent methods in terms of numerical performance. Furthermore, we observe that AGCD attains O(1/k2)O(1/k^{2}) convergence (or better) on these problems in practice. Thus AGCD is worthy of significant further study, both computationally as well as theoretically. Indeed, in Section 4 we will present a technical condition that implies O(1/k2)O(1/k^{2}) convergence when satisfied, and we will argue that this condition ought to be satisfied in many settings.

ASCD, which we consider to be the new theoretical contribution of this paper, combines the salient features of AGCD and ARCD. In ASCD we choose the greedy coordinate of the gradient to perform the xx-update, while we choose a random coordinate to perform the zz-update. In this way we achieve the practical advantage of greedy xx-updates, while still guaranteeing O(1/k2)O(1/k^{2}) convergence in expectation by virtue of choosing the random coordinate used in the zz-update, see Theorem 2.1. And under strong convexity, ASCD achieves linear convergence as will be shown in Section 3.

The paper is organized as follows. In Section 2 we present the O(1/k2)O(1/k^{2}) convergence guarantee (in expectation) for ASCD. In Section 3 we present an extension of the accelerated coordinate descent framework to the case of strongly convex functions, and we present the associated linear convergence guarantee for ASCD under strong convexity. In Section 4 we study AGCD; we present a Lyapunov energy function argument that points to why standard analysis of accelerated gradient descent methods fails in the analysis of AGCD. In Section 4.2 we present a technical condition under which AGCD will achieve O(1/k2)O(1/k^{2}) convergence. In Section 5, we present results of our numerical experiments using AGCD and ASCD on synthetic linear regression problems as well as practical logistic regression problems.

Accelerated Semi-Greedy Coordinate Descent (ASCD)

In this section we present our computational guarantee for the Accelerated Semi-Greedy Coordinate Descent (ASCD) method in the non-strongly convex case. Algorithm 1 with rule (5) presents the Accelerated Semi-Greedy Coordinate Descent method (ASCD) for the non-strongly convex case. At each iteration kk the ASCD method choose the greedy coordinate jk1j_{k}^{1} to do the xx-update, and chooses a randomized coordinate jk2∼U[1,⋯ ,n]j_{k}^{2}\sim{\mathcal{U}}[1,\cdots,n] to do the zz-update. Unlike ARCD where the same randomized coordinate is used in both the xx-update and the zz-update – in ASCD jk1j_{k}^{1} is chosen in a deterministic greedy way, jk1j_{k}^{1} and jk2j_{k}^{2} are likely to be different.

At each iteration kk of ASCD the random variable jk2j_{k}^{2} is introduced, and therefore xkx^{k} depends on the realization of the random variable

For convenience we also define ξ0:=∅\xi_{0}:=\emptyset.

The following theorem presents our computational guarantee for ASCD for the non-strongly convex case:

Consider the Accelerated Semi-Greedy Coordinate Descent method (Algorithm 1 with rule (5)). If f(⋅)f(\cdot) is coordinate-wise LL-smooth, it holds for all k≥1k\geq 1 that:

In the interest of both clarity and a desire to convey some intuition on proofs of accelerated methods in general, we will present the proof of Theorem 2.1 after first establishing some intermediary results along with some explanatory comments. We start with the “Three-Point Property” of Tseng . Given a differentiable convex function h(⋅)h(\cdot), the Bregman distance for h(⋅)h(\cdot) is Dh(y,x):=h(y)−h(x)−⟨∇h(x),y−x⟩D_{h}(y,x):=h(y)-h(x)-\langle\nabla h(x),y-x\rangle. The Three-Point property can be stated as follows:

(Three-Point Property (Tseng )) Let ϕ(⋅)\phi(\cdot) be a convex function, and let Dh(⋅,⋅)D_{h}(\cdot,\cdot) be the Bregman distance for h(⋅)h(\cdot). For a given vector zz, let

with equality holding in the case when ϕ(⋅)\phi(\cdot) is a linear function and h(⋅)h(\cdot) is a quadratic function. ∎

Also, it follows from elementary integration and the coordinate-wise Lipschitz condition (2) that

At each iteration k=0,1,…k=0,1,\ldots of ASCD, notice that xk+1x^{k+1} is one step of greedy coordinate descent from yky^{k} in the norm ∥⋅∥L\|\cdot\|_{L}. Now define sk+1:=yk−1nL−1∇f(yk)s^{k+1}:=y^{k}-\tfrac{1}{n}{\bf{L}}^{-1}\nabla f(y^{k}), which is a full steepest-descent step from yky^{k} in the norm ∥⋅∥nL\|\cdot\|_{nL} . We first show that the greedy coordinate descent step yields a good objective function value as compared to the quadratic model that yields sk+1s^{k+1}.

where the first inequality of (9) derives from the smoothness of f(⋅)f(\cdot), and is a simple instance of (8) using x=ykx=y^{k}, i=jk1i=j_{k}^{1}, and h=−1Ljk1∇jk1f(yk)h=-\tfrac{1}{L_{j_{k}^{1}}}\nabla_{j_{k}^{1}}f(y^{k}). The second inequality of (9) follows from the definition of jk1j_{k}^{1} which yields:

The last equality of (9) follows by using the definition of sk+1s^{k+1} and rearranging terms. ∎

Utilizing the interpretation of sk+1s^{k+1} as a gradient descent step from yky^{k} but with a larger smoothness descriptor (nLnL as opposed to LL), we can invoke the standard proof for accelerated gradient descent derived in for example. We define tk+1:=zk−1nθkL−1∇f(yk)t^{k+1}:=z^{k}-\tfrac{1}{n\theta_{k}}{\bf{L}}^{-1}\nabla f(y^{k}), or equivalently we can define tk+1t^{k+1} by:

(which corresponds to zk+1z^{k+1} in for standard accelerated gradient descent). Then we have:

where the first equality of (12) utilizes sk+1−yk=θk(tk+1−zk)s^{k+1}-y^{k}=\theta_{k}(t^{k+1}-z^{k}). The second equality of (12) follows as an application of the Three-Point-Property (Lemma 2.1) together with (10), where we set ϕ(x)=⟨∇f(yk),x−zk⟩\phi(x)=\langle\nabla f(y^{k}),x-z^{k}\rangle and h(x)=nθk2∥x∥L2h(x)=\tfrac{n\theta_{k}}{2}\|x\|_{L}^{2} (whereby Dh(x,v)=nθk2∥x−v∥L2D_{h}(x,v)=\tfrac{n\theta_{k}}{2}\|x-v\|_{L}^{2}). The third equality of (12) is derived from yk=(1−θk)xk+θkzky^{k}=(1-\theta_{k})x^{k}+\theta_{k}z^{k} and rearranging the terms. And the last inequality of (12) is an application of the gradient inequality at yky^{k} applied to xkx^{k} and also to x∗x^{*}. ∎

Notice that tk+1t^{k+1} is an all-coordinate update of zkz^{k}, and computing tk+1t^{k+1} can be very expensive. Instead we will use zk+1z^{k+1} to replace tk+1t^{k+1} in (11) by using the equality in the next lemma.

where the first and third equations above are straightforward arithmetic rearrangements, and the second equation follows from the two easy-to-verify identities tk+1−zk=nEjk2[zk+1−zk]t^{k+1}-z^{k}=nE_{j_{k}^{2}}\left[z^{k+1}-z^{k}\right] and ∥tk+1−zk∥L2=nEjk2[∥zk+1−zk∥L2]\left\|t^{k+1}-z^{k}\right\|_{L}^{2}=nE_{j_{k}^{2}}\left[\left\|z^{k+1}-z^{k}\right\|_{L}^{2}\right] . ∎

We now have all of the ingredients needed to prove Theorem 2.1.

Proof of Theorem 2.1 Substituting (13) into (11), we obtain:

Rearranging and substituting 1−θk+1θk+12=1θk2\tfrac{1-\theta_{k+1}}{\theta_{k+1}^{2}}=\tfrac{1}{\theta_{k}^{2}}, we arrive at:

Taking the expectation over the random variables j12,j22,…,jk2j_{1}^{2},j_{2}^{2},\ldots,j_{k}^{2}, it follows that:

Applying the above inequality in a telescoping manner for k=1,2,…k=1,2,\ldots, yields:

Note from an induction argument that θi≤2i+2\theta_{i}\leq\tfrac{2}{i+2} for all =0,1,…=0,1,\ldots, whereby the above inequality rearranges to:

Accelerated Coordinate Descent Framework under Strong Convexity

We begin with the definition of strong convexity as developed in :

Note that μ\mu can be viewed as an extension of the condition number of f(⋅)f(\cdot) in the traditional sense since μ\mu is defined relative to the coordinate smoothness coefficients through ∥⋅∥L\|\cdot\|_{L}, see . Algorithm 2 presents the generic framework for accelerated coordinate descent methods in the case when f(⋅)f(\cdot) is μ\mu-strongly convex for known μ\mu.

Just as in the non-strongly convex case, we extend the three algorithms ARCD, AGCD, and ASCD to the strongly convex case by using the rules (3), (4), and (5) in Algorithm 2. The following theorem presents our computational guarantee for ASCD for strongly convex case:

Consider the Accelerated Semi-Greedy Coordinate Descent method for strongly convex case (Algorithm 2 with rule (5)). If f(⋅)f(\cdot) is coordinate-wise LL-smooth and μ\mu-strongly convex with respect to ∥⋅∥L\|\cdot\|_{L}, it holds for all k≥1k\geq 1 that:

We provide a concise proof of Theorem 3.1 in the Appendix.

Accelerated Greedy Coordinate Descent

In this section we discuss accelerated greedy coordinate descent (AGCD), which is Algorithm 1 with rule (4). In the interest of clarity we limit our discussion to the non-strongly convex case. We present a Lyapunov function argument which shows why the standard type of proof of accelerated gradient methods fails for AGCD, and we propose a technical condition under which AGCD is guaranteed to have an O(1/k2)O(1/k^{2}) accelerated convergence rate. Although there are no guarantees that the technical condition will hold for a given function f(⋅)f(\cdot), we provide intuition as to why the technical condition ought to hold in most cases.

The mainstream research community’s interest in Nesterov’s accelerated method started around 15 years ago; and yet even today most researchers struggle to find basic intuition as to what is really going on in accelerated methods. Indeed, Nesterov’s estimation sequence proof technique seems to work out arithmetically but with little fundamental intuition. There are many recent work trying to explain this acceleration phenomenon . A line of recent work has attempted to give a physical explanation of acceleration techniques by studying the continuous-time interpretation of accelerated gradient descent via dynamical systems, see , , and . In particular, introduced the continuous-time dynamical system model for accelerated gradient descent, and presented a convergence analysis using a Lyapunov energy function in the continuous-time setting. studied discretizations of the continuous-time dynamical system, and also showed that Nesterov’s estimation sequence analysis is equivalent to the Lyapunov energy function analysis in the dynamical system in the discrete-time setting. And presented an energy dissipation argument from control theory for understanding accelerated gradient descent.

In the discrete-time setting, one can construct a Lyapunov energy function of the form :

where AkA_{k} is a parameter sequence with Ak∼O(k2)A_{k}\sim O(k^{2}), and one shows that EkE_{k} is non-increasing in kk, thereby yielding:

The proof techniques of acceleration methods such as , and , as well as the recent proof techniques for accelerated randomized coordinate descent (such as , , and ) can all be rewritten in the above form (up to expectation) each with slightly different parameter sequences {Ak}\{A_{k}\}.

Now let us return to accelerated greedy coordinate descent. Let us assume for simplicity that L1=⋯=LnL_{1}=\cdots=L_{n} (as we can always do rescaling to achieve this condition). Then the greedy coordinate jk1j_{k}^{1} is chosen as the coordinate of the gradient with the largest magnitude, which corresponds to the coordinate yielding the greatest guaranteed decrease in the objective function value. However, in the proof of acceleration using the Lyapunov energy function, one needs to prove a decrease in EkE_{k} (17) instead of a decrease in the objective function value f(xk)f(x^{k}). The coordinate jk1j_{k}^{1} is not necessarily the greedy coordinate for decreasing the energy function EkE_{k} due to the presence of the second term ∥x∗−zk∥L2\|x^{*}-z^{k}\|_{L}^{2} in (17). This explains why the greedy coordinate can fail to decrease EkE_{k}, at least in theory. And because x∗x^{*} is not known when running AGCD, there does not seem to be any way to find the greedy descent coordinate for the energy function EkE_{k}.

That is why in ASCD we use the greedy coordinate to perform the xx-update (which corresponds to the fastest coordinate-wise decrease for the first term in energy function), while we choose a random coordinate to perform the zz-update (which corresponds to the second term in the energy function); thereby mitigating the above problem in the case of ASCD.

In a concurrent paper , the authors develop computational theory for matching pursuit algorithms, which can be viewed as a generalized version of greedy coordinate descent where the directions do not need to form an orthogonal basis. The paper also develops an accelerated version of the matching pursuit algorithms, which turns out to be equivalent to the algorithm ASCD discussed here in the special case where the chosen directions are orthogonal. Although the focus in and in our paper are different – is more focused on (accelerated) greedy direction updates along a certain linear subspace whereas our focus is on when and how one can accelerate greedy coordinate updates – both of the works share a similar spirit and similar approaches in developing accelerated greedy methods. Moreover, both works use a decoupling of the coordinate update for the {xk}\{x^{k}\} sequence (with a greedy rule) and the {zk}\{z^{k}\} sequence (with a randomized rule). In fact, is consistent with the argument in our paper as to why one cannot accelerate greedy coordinate descent in general.

2 How to make AGCD work (in theory)

Here we propose the following technical condition under which the proof of acceleration of AGCD can be made to work.

There exists a positive constant γ\gamma and an iteration number KK such that for all k≥Kk\geq K it holds that:

where ji=arg⁡max⁡i1Li∣∇if(yk)∣j_{i}=\arg\max_{i}\tfrac{1}{\sqrt{L_{i}}}|\nabla_{i}f(y^{k})| is the greedy coordinate at iteration ii.

One can show that this condition is sufficient to prove an accelerated convergence rate O(1/k2)O(1/k^{2}) for AGCD. Therefore let us take a close look at Technical Condition 4.1. The condition considers the weighted sum (with weights 1θi∼O(i2)\frac{1}{\theta_{i}}\sim O(i^{2})) of the inner product of ∇f(yk)\nabla f(y^{k}) and zk−x∗z^{k}-x^{*}, and the condition states that the inner product corresponding to the greedy coordinate (the right side above) is larger than the average of all coordinates in the inner product, by a factor of γ\gamma. In the case of ARCD and ASCD, it is easy to show that Technical Condition 4.1 holds automatically up to expectation, with γ=1\gamma=1.

Here is an informal explanation of why Technical Condition 4.1 ought to hold for most convex functions and most iterations of AGCD. When kk is sufficiently large, the three sequence {xk}\{x^{k}\}, {yk}\{y^{k}\} and {zk}\{z^{k}\} ought to all converge to x∗x^{*} (which always happens in practice though lack of theoretical justification), whereby zkz^{k} is close to yky^{k}. Thus we can instead consider the inner product ⟨∇f(yk),yk−x∗⟩\langle\nabla f(y^{k}),y^{k}-x^{*}\rangle in (18). Notice that for any coordinate jj it holds that ∣yjk−xj∗∣≥1Lj∣∇jf(yk)∣|y^{k}_{j}-x^{*}_{j}|\geq\frac{1}{L_{j}}|\nabla_{j}f(y^{k})|, and therefore ∣∇jf(yk)⋅(yjk−xj∗)∣≥1Lj∣∇jf(jk)∣2|\nabla_{j}f(y^{k})\cdot(y^{k}_{j}-x^{*}_{j})|\geq\frac{1}{L_{j}}|\nabla_{j}f(j^{k})|^{2}. Now the greedy coordinate is chosen by ji:=arg⁡max⁡j1Lj∣∇jf(jk)∣2j_{i}:=\arg\max_{j}\frac{1}{L_{j}}|\nabla_{j}f(j^{k})|^{2}, and therefore it is reasonably likely that in most cases the greedy coordinate will yield a better product than the average of the components of the inner product.

The above is not a rigorous argument, and we can likely design some worst-case functions for which Technical Condition 4.1 fails. But the above argument provides some intuition as to why the condition ought to hold in most cases, thereby yielding the observed improvement of AGCD as compared with ARCD that we will shortly present in Section 5, where we also observe that Technical Condition 4.1 holds empirically on all of our problem instances.

With a slight change in the proof of Theorem 2.1, we can show the following result:

Consider the Accelerated Greedy Coordinate Descent (Algorithm 1 with rule (4)). If f(⋅)f(\cdot) is coordinate-wise LL-smooth and satisfies Technical Condition 4.1 with constant γ≤1\gamma\leq 1 and iteration number KK, then it holds for all k≥Kk\geq K that:

We note that if γ<1\gamma<1 (which we always observe in practice), then AGCD will have a better convergence guarantee than ARCD.

The arguments in Section 4.1 and Section 4.2 also work for strongly convex case, albeit with suitable minor modifications.

Numerical Experiments

We consider solving synthetic instances of the linear regression model with least-squares objective function:

using ASCD, ARCD and AGCD, where the mechanism for generating the data (y,X)(y,X) and the algorithm implementation details are described in the supplementary materials. Figure 1 shows the optimality gap versus time (in seconds) for solving different instances of linear regression with different condition numbers of the matrix XTXX^{T}X using ASCD, ARCD and AGCD. In each plot, the vertical axis is the objective value optimality gap f(βk)−f∗f(\beta^{k})-f^{*} in log scale, and the horizontal axis is the running time in seconds. Each column corresponds to an instance with the prescribed condition number κ\kappa of XTXX^{T}X, where κ=∞\kappa=\infty means that the minimum eigenvalue of XTXX^{T}X is 00. The first row of plots is for Algorithm Framework 1 which is ignorant of any strong convexity information. The second row of plots is for Algorithm Framework 2, which uses given strong convexity information. And because the linear regression optimization problem is quadratic, it is straightforward to compute κ\kappa as well as the true parameter μ\mu for the instances where κ>0\kappa>0. The last column of the figure corresponds to κ=∞\kappa=\infty, and in this instance we set μ\mu using the smallest positive eigenvalue of XTXX^{T}X, which can be shown to work in theory since all relevant problem computations are invariant in the nullspace of XX.

Here we see in Figure 1 that AGCD and ASCD consistently have superior performance over ARCD for both the non-strongly convex case and the strongly convex case, with ASCD performing almost as well as AGCD in most instances.

We remark that the behavior of any convex function near the optimal solution is similar to the quadratic function defined by the Hessian at the optimum, and therefore the above numerical experiments show promise that AGCD and ASCD are likely to outperform ARCD asymptotically for any twice-differentiable convex function.

2 Logistic Regression

Here we consider solving instances of the logistic regression loss minimization problem:

using ASCD, ARCD and AGCD, where {xi,yi}\{x_{i},y_{i}\} is the feature-response pair for the ii-th data point and yi∈{−1,1}y_{i}\in\{-1,1\}. Although the loss function f(β)f(\beta) is not in general strongly convex, it is essentially locally strongly convex around the optimum but with unknown strong convexity parameter μˉ\bar{\mu}. And although we do not know the local strong convexity parameter μˉ\bar{\mu}, we can still run the strongly convex algorithm (Algorithm Framework 2) by assigning a value of μˉ\bar{\mu} that is hopefully close to the actual value. Using this strategy, we solved a large number of logistic regression instances from LIBSVM . Figure 2 shows the optimality gap versus time (in seconds) for solving two of these instances, namely w1a and a1a, which were chosen here because the performance of the algorithms on these two instances is representative of others in LIBSVM. In each plot, the vertical axis is the objective value optimality gap f(βk)−f∗f(\beta^{k})-f^{*} in log scale, and the horizontal axis is the running time in seconds. Each column corresponds to a different assigned value of the local strong convexity parameter μˉ\bar{\mu}. The right-most column in the figure uses the assignment μˉ=0\bar{\mu}=0, in which case the algorithms are implemented as in the non-strongly convex case (Algorithm Framework 1).

Here we see in Figure 2 that AGCD always has superior performance as compared to either ASCD and ARCD. In the relevant range of optimality gaps (≤10−9\leq 10^{-9}), ASCD typically outperforms ARCD for smaller values of the assigned strong convexity parameter μˉ\bar{\mu}. However, the performance of ASCD and ARCD are essentially the same when no strong convexity is presumed.

Last of all, we attempt to estimate the parameter γ\gamma that arises in Technical Condition 4.1 for AGCD in several of the datasets in SVMLIB. Although for small kk, the ratio between ∑i=0k1θi⟨∇f(yi),zi−x∗⟩\sum_{i=0}^{k}\frac{1}{\theta_{i}}\langle\nabla f(y^{i}),z^{i}-x^{*}\rangle and ∑i=0knθi∇jif(yi)(zjii−xji∗)\sum_{i=0}^{k}\frac{n}{\theta_{i}}\nabla_{j_{i}}f(y^{i})(z^{i}_{j_{i}}-x^{*}_{j_{i}}) can fluctuate widely, this ratio stabilizes after a number of iterations in all of our numerical tests. From Technical Condition 4.1, we know that γ\gamma is the upper bound of such ratio for all k≥Kk\geq K for some large enough value of KK. Table 1 presents the observed values of γ\gamma for all K≥Kˉ:=5,000K\geq\bar{K}:=5,000. Recalling from Theorem 4.1 that the γ\gamma value represents how much better AGCD can perform compared with ARCD in terms of computational guarantees, we see from Table 1 that AGCD should outperform ARCD for these representative instances – and indeed this is what we observe in practice in our tests.

Appendix A Appendix

In order to prove Theorem 3.1, we first prove the following three lemmas:

where the second equality utilizes uk=a2a2+bzk+ba2+byku^{k}=\frac{a^{2}}{a^{2}+b}z^{k}+\frac{b}{a^{2}+b}y^{k} and the other equalities are just mathematical manipulations. ∎

Define tk+1:=uk−aa2+b1nL−1∇f(yk)t^{k+1}:=u^{k}-\frac{a}{a^{2}+b}\frac{1}{n}{\bf{L}}^{-1}\nabla f(y^{k}), then

where the second equality is from the relationship of tk+1=uk−aa2+b1nL−1∇f(yk)t^{k+1}=u^{k}-\frac{a}{a^{2}+b}\frac{1}{n}{\bf{L}}^{-1}\nabla f(y^{k}) and zk+1=uk−aa2+b1nLjk2∇jk2f(yk)ejk2z^{k+1}=u^{k}-\frac{a}{a^{2}+b}\frac{1}{nL_{j_{k}^{2}}}\nabla_{j_{k}^{2}}f(y^{k})e_{j_{k}^{2}}, and the first and third equations are just rearrangement. ∎

Proof: Remember that b=μan2b=\frac{\mu a}{n^{2}} and a>0a>0, thus the above inequality is equivalent to

Substituting a=μn+μa=\tfrac{\sqrt{\mu}}{n+\sqrt{\mu}}, the above inequality becomes

We furnish the proof by noting μn−μn+μ=μn(n+μ)≤μn2\tfrac{\sqrt{\mu}}{n}-\tfrac{\sqrt{\mu}}{n+\sqrt{\mu}}=\tfrac{\mu}{n(n+\sqrt{\mu})}\leq\tfrac{\mu}{n^{2}}. ∎

Proof of Theorem 3.1: Recall that tk+1=uk−aa2+b1nL−1∇f(yk)t^{k+1}=u^{k}-\frac{a}{a^{2}+b}\frac{1}{n}{\bf{L}}^{-1}\nabla f(y^{k}), then it is easy to check that

by writing the optimality conditions of the right-hand side.

where the first inequality is due to coordinate-wise smoothness, the first equality utilizes xk+1=yk−1Ljk1∇jk1f(yk)ejk1x^{k+1}=y^{k}-\frac{1}{L_{j_{k}^{1}}}\nabla_{j_{k}^{1}}f(y^{k})e_{j_{k}^{1}}, the second inequality follows from the fact that jk1j_{k}^{1} is the greedy coordinate of ∇f(yk)\nabla f(y^{k}) in the ∥⋅∥L−1\|\cdot\|_{L^{-1}} norm, the third inequality follows from the basic inequality ∥v∥L2+∥w∥L−12≥2⟨v,w⟩\|v\|_{L}^{2}+\|w\|_{L^{-1}}^{2}\geq 2\langle v,w\rangle for all v,wv,w, the second equality is from Three Point Property by noticing

the third equality follows from Lemma 3.1, and the fourth and sixth equalities each utilize Lemma 3.2.

On the other hand, by strong convexity we have

where the second equality uses the fact that yk=(1−a)xk+azky^{k}=(1-a)x^{k}+az^{k} and the last inequality is from the gradient inequality.

Notice that b=μan2b=\frac{\mu a}{n^{2}} and a2≤(1−a)(a2+b)a^{2}\leq(1-a)(a^{2}+b) following from Lemma 3.3. Thus summing up (23) and (25) leads to

which furnishes the proof using a telescoping series. ∎

Appendix B More Material on the Numerical Experiments

To be consistent with the notation in statistics and machine learning we use pp to denote the dimension of the variables in the optimization problems describing linear and logistic regression. Then the per-iteration computation cost of AGCD and ASCD is dominated by three computations: (i) pp-dimensional vector operations (such as in computing yky^{k} using xkx^{k} and zkz^{k}), (ii) computation of the gradient ∇f(⋅)\nabla f(\cdot), and (iii) computation of the maximum (weighted) magnitude coordinate of the gradient ∇f(⋅)\nabla f(\cdot). proposed an efficient way to avoid (i) by changing variables. Distinct from the dual approaches discussed in , , and , here we only consider the primal problem in the regime n>pn>p, and therefore the cost of (ii) dominates the cost of (i) in these cases. For this reason in our numerical experiments we use the simple implementation of ARCD proposed by Nesterov and which we adopt for AGCD and ASCD as well. We note that both the randomized methods and the greedy methods can take advantage of the efficient calculations proposed in as well.

For the linear regression experiments we focused on synthetic problem instances with different condition numbers κ\kappa of the matrix XTXX^{T}X and where XX is dense. In this case the cost for computation (i) is O(p)O(p). And by taking advantage of the coordinate update structure, we can implement (ii) in O(p)O(p) operations by pre-computing and storing XTXX^{T}X in memory, see and for further details. The cost of (iii) is simply O(p)O(p).

For the logistic regression experiments, the cost of (ii) at each iteration of AGCD and ASCD can be much larger than O(p)O(p) because there is no easy way to update the full gradient ∇f(⋅)\nabla f(\cdot). For these problems we have

where XX is the sample matrix with xix_{i} composing the ii-th row, and w(β)i:=11+exp⁡(yiβTxi)w(\beta)_{i}:=\frac{1}{1+\exp(y_{i}\beta^{T}x_{i})}. Notice that calculating w(β)w(\beta) can be done using a rank 11-update with cost O(n)O(n). But calculating the matrix-vector product XTw(β)X^{T}w(\beta) will cost O(np)O(np), which dominates the cost of (i) and/or (iii). However, in the case when XX is a sparse matrix with density ρ\rho, the cost can be decreased to O(ρnp)O(\rho np).

B.2 Comparing the Algorithms using Running Time and the Number of Iterations

Figure 3 shows the optimality gap versus running time (seconds) in the left plot and and versus the number of iterations in the right plot, logistic regression problem using the dataset madelon in LIBSVM , with μˉ=10−7\bar{\mu}=10^{-7}. Here we see that AGCD and ASCD are vastly superior to ARCD in term of the number of iterations, but not nearly as much in terms of running time, because one iteration of AGCD or ASCD can be more expensive than an iteration of ARCD.

B.3 Comparing Accelerated Method with Non-Accelerated Method

Figure 4 shows the optimality gap versus running time (seconds) with GCD, ASCD and AGCD for logistic regression problem using the dataset madelon in LIBSVM , with μˉ=10−6\bar{\mu}=10^{-6}. Here we see that ASCD and AGCD are superior to non-accelerated GCD.

B.4 Numerical Results for Logistic Regression with Other Datasets

We present numerical results for logistic regression problems for several other datasets in LIBSVM solved by ASCD, ARCD and AGCD in Figure 5. Here we see that AGCD always has superior performance as compared to ASCD and ARCD, and ASCD outperforms ARCD in most of the cases.

References