Coordinate Descent Converges Faster with the Gauss-Southwell Rule Than Random Selection

Coordinate Descent Methods

There has been substantial recent interest in applying coordinate descent methods to solve large-scale optimization problems, starting with the seminal work of Nesterov (2012), who gave the first global rate-of-convergence analysis for coordinate-descent methods for minimizing convex functions. This analysis suggests that choosing a random coordinate to update gives the same performance as choosing the “best” coordinate to update via the more expensive Gauss-Southwell (GS) rule. (Nesterov also proposed a more clever randomized scheme, which we consider later in this paper.) This result gives a compelling argument to use randomized coordinate descent in contexts where the GS rule is too expensive. It also suggests that there is no benefit to using the GS rule in contexts where it is relatively cheap. But in these contexts, the GS rule often substantially outperforms randomized coordinate selection in practice. This suggests that either the analysis of GS is not tight, or that there exists a class of functions for which the GS rule is as slow as randomized coordinate descent.

After discussing contexts in which it makes sense to use coordinate descent and the GS rule, we answer this theoretical question by giving a tighter analysis of the GS rule (under strong-convexity and standard smoothness assumptions) that yields the same rate as the randomized method for a restricted class of functions, but is otherwise faster (and in some cases substantially faster). We further show that, compared to the usual constant step-size update of the coordinate, the GS method with exact coordinate optimization has a provably faster rate for problems satisfying a certain sparsity constraint (Section 5). We believe that this is the first result showing a theoretical benefit of exact coordinate optimization; all previous analyses show that these strategies obtain the same rate as constant step-size updates, even though exact optimization tends to be faster in practice. Furthermore, in Section 6, we propose a variant of the GS rule that, similar to Nesterov’s more clever randomized sampling scheme, uses knowledge of the Lipschitz constants of the coordinate-wise gradients to obtain a faster rate. We also analyze approximate GS rules (Section 7), which provide an intermediate strategy between randomized methods and the exact GS rule. Finally, we analyze proximal-gradient variants of the GS rule (Section 8) for optimizing problems that include a separable non-smooth term.

Problems of Interest

The rates of Nesterov show that coordinate descent can be faster than gradient descent in cases where, if we are optimizing $n$ variables, the cost of performing $n$ coordinate updates is similar to the cost of performing one full gradient iteration. This essentially means that coordinate descent methods are useful for minimizing convex functions that can be expressed in one of the following two forms:

where $x_{i}$ is element $i$ of $x$, $f$ is smooth and cheap, the $f_{ij}$ are smooth, $G=\{V,E\}$ is a graph, and $A$ is a matrix. (It is assumed that all functions are convex.)We could also consider slightly more general cases like functions that are defined on hyper-edges (Richtárik and Takáč, 2015), provided that we can still perform $n$ coordinate updates for a similar cost to one gradient evaluation. The family of functions $h_{1}$ includes core machine-learning problems such as least squares, logistic regression, lasso, and SVMs (when solved in dual form) (Hsieh et al., 2008). Family $h_{2}$ includes quadratic functions, graph-based label propagation algorithms for semi-supervised learning (Bengio et al., 2006), and finding the most likely assignments in continuous pairwise graphical models (Rue and Held, 2005).

In general, the GS rule for problem $h_{2}$ is as expensive as a full gradient evaluation. However, the structure of $G$ often allows efficient implementation of the GS rule. For example, if each node has at most $d$ neighbours, we can track the gradients of all the variables and use a max-heap structure to implement the GS rule in $O(d\log n)$ time (Meshi et al., 2012). This is similar to the cost of the randomized algorithm if $d\approx|E|/n$ (since the average cost of the randomized method depends on the average degree). This condition is true in a variety of applications. For example, in spatial statistics we often use two-dimensional grid-structured graphs, where the maximum degree is four and the average degree is slightly less than $4$. As another example, for applying graph-based label propagation on the Facebook graph (to detect the spread of diseases, for example), the average number of friends is around $200$ but no user has more than seven thousand friends.https://recordsetter.com/world-record/facebook-friends The maximum number of friends would be even smaller if we removed edges based on proximity. A non-sparse example where GS is efficient is complete graphs, since here the average degree and maximum degree are both $(n-1)$. Thus, the GS rule is efficient for optimizing dense quadratic functions. On the other hand, GS could be very inefficient for star graphs.

If each column of $A$ has at most $c$ non-zeroes and each row has at most $r$ non-zeroes, then for many notable instances of problem $h_{1}$ we can implement the GS rule in $O(cr\log n)$ time by maintaining $Ax$ as well as the gradient and again using a max-heap (see Appendix A). Thus, GS will be efficient if $cr$ is similar to the number of non-zeroes in $A$ divided by $n$. Otherwise, Dhillon et al. (2011) show that we can approximate the GS rule for problem $h_{1}$ with no $g_{i}$ functions by solving a nearest-neighbour problem. Their analysis of the GS rule in the convex case, however, gives the same convergence rate that is obtained by random selection (although the constant factor can be smaller by a factor of up to $n$). More recently, Shrivastava and Li (2014) give a general method for approximating the GS rule for problem $h_{1}$ with no $g_{i}$ functions by writing it as a maximum inner-product search problem.

Existing Analysis

We are interested in solving the convex optimization problem

where $\nabla f$ is coordinate-wise $L$-Lipschitz continuous, i.e., for each $i=1,\ldots,n$,

where $e_{i}$ is a vector with a one in position $i$ and zero in all other positions. For twice-differentiable functions, this is equivalent to the assumption that the diagonal elements of the Hessian are bounded in magnitude by $L$. In contrast, the typical assumption used for gradient methods is that $\nabla f$ is $L^{f}$-Lipschitz continuous (note that $L\leq L^{f}\leq Ln$). The coordinate-descent method with constant step-size is based on the iteration

The randomized coordinate-selection rule chooses $i_{k}$ uniformly from the set $\{1,2,\dots,n\}$. Alternatively, the GS rule

chooses the coordinate with the largest directional derivative. Under either rule, because $f$ is coordinate-wise Lipschitz continuous, we obtain the following bound on the progress made by each iteration:

We focus on the case where $f$ is $\mu$-strongly convex, meaning that, for some positive $\mu$,

where $x^{*}$ is the optimal solution of (1). This bound is obtained by minimizing both sides of (3) with respect to $y$.

Conditioning on the $\sigma$-field $\mathcal{F}_{k-1}$ generated by the sequence $\{x^{0},x^{1},\ldots,x^{k-1}\}$, and taking expectations of both sides of (2), when $i_{k}$ is chosen with uniform sampling we obtain

Using (4) and subtracting $f(x^{*})$ from both sides, we get

This is a special of case of Nesterov (2012, Theorem 2) with $\alpha=0$ in his notation.

2 Gauss-Southwell

We now consider the progress implied by the GS rule. By the definition of $i_{k}$,

Applying this inequality to (2), we obtain

This is a special case of Boyd and Vandenberghe (2004, §9.4.3), viewing the GS rule as performing steepest descent in the $1$-norm. While this is faster than known rates for cyclic coordinate selection (Beck and Tetruashvili, 2013) and holds deterministically rather than in expectation, this rate is the same as the randomized rate given in (5).

Refined Gauss-Southwell Analysis

The deficiency of the existing GS analysis is that too much is lost when we use the inequality in (6). To avoid the need to use this inequality, we instead measure strong-convexity in the $1$-norm, i.e.,

which is the analogue of (3). Minimizing both sides with respect to $y$, we obtain

which makes use of the convex conjugate $(\frac{\mu_{1}}{2}\|\cdot\|_{1}^{2})^{*}=\frac{1}{2\mu_{1}}\|\cdot\|_{\infty}^{2}$ (Boyd and Vandenberghe, 2004, §3.3). Using (8) in (2), and the fact that $(\nabla_{i_{k}}f(x^{k}))^{2}=\|\nabla f(x^{k})\|_{\infty}^{2}$ for the GS rule, we obtain

It is evident that if $\mu_{1}=\mu/n$, then the rates implied by (5) and (9) are identical, but (9) is faster if $\mu_{1}>\mu/n$. In Appendix B, we show that the relationship between $\mu$ and $\mu_{1}$ can be obtained through the relationship between the squared norms $||\cdot||^{2}$ and $||\cdot||_{1}^{2}$. In particular, we have

Thus, at one extreme the GS rule obtains the same rate as uniform selection ($\mu_{1}\approx\mu/n$). However, at the other extreme, it could be faster than uniform selection by a factor of $n$ ($\mu_{1}\approx\mu$). This analysis, that the GS rule only obtains the same bound as random selection in an extreme case, supports the better practical behaviour of GS.

We illustrate these two extremes with the simple example of a quadratic function with a diagonal Hessian $\nabla^{2}f(x)=\hbox{diag}({\lambda_{1},\ldots,\lambda_{n}})$. In this case,

We prove the correctness of this formula for $\mu_{1}$ in Appendix C. The parameter $\mu_{1}$ achieves its lower bound when all $\lambda_{i}$ are equal, $\lambda_{1}=\cdots=\lambda_{n}=\alpha>0$, in which case

Thus, uniform selection does as well as the GS rule if all elements of the gradient change at exactly the same rate. This is reasonable: under this condition, there is no apparent advantage in selecting the coordinate to update in a clever way. Intuitively, one might expect that the favourable case for the Gauss-Southwell rule would be where one $\lambda_{i}$ is much larger than the others. However, in this case, $\mu_{1}$ is again similar to $\mu/n$. To achieve the other extreme, suppose that $\lambda_{1}=\beta$ and $\lambda_{2}=\lambda_{3}=\cdots=\lambda_{n}=\alpha$ with $\alpha\geq\beta$. In this case, we have $\mu=\beta$ and

If we take $\alpha\to\infty$, then we have $\mu_{1}\to\beta$, so $\mu_{1}\to\mu$. This case is much less intuitive; GS is $n$ times faster than random coordinate selection if one element of the gradient changes much more slowly than the others.

2 ‘Working Together’ Interpretation

In the separable quadratic case above, $\mu_{1}$ is given by the harmonic mean of the eigenvalues of the Hessian divided by $n$. The harmonic mean is dominated by its smallest values, and this is why having one small value is a notable case. Furthermore, the harmonic mean divided by $n$ has an interpretation in terms of processes ‘working together’ (Ferger, 1931). If each $\lambda_{i}$ represents the time taken by each process to finish a task (e.g., large values of $\lambda_{i}$ correspond to slow workers), then $\mu$ is the time needed by the fastest worker to complete the task, and $\mu_{1}$ is the time needed to complete the task if all processes work together (and have independent effects). Using this interpretation, the GS rule provides the most benefit over random selection when working together is not efficient, meaning that if the $n$ processes work together, then the task is not solved much faster than if the fastest worker performed the task alone. This gives an interpretation of the non-intuitive scenario where GS provides the most benefit: if all workers have the same efficiency, then working together solves the problem $n$ times faster. Similarly, if there is one slow worker (large $\lambda_{i}$), then the problem is solved roughly $n$ times faster by working together. On the other hand, if most workers are slow (many large $\lambda_{i}$), then working together has little benefit.

3 Fast Convergence with Bias Term

Consider the standard linear-prediction framework,

where we have included a bias variable $\beta$ (an example of problem $h_{1}$). Typically, the regularization parameter $\sigma$ of the bias variable is set to be much smaller than the regularization parameter $\lambda$ of the other covariates, to avoid biasing against a global shift in the predictor. Assuming that there is no hidden strong-convexity in the sum, this problem has the structure described in the previous section ($\mu_{1}\approx\mu$) where GS has the most benefit over random selection.

Rates with Different Lipschitz Constants

Consider the more general scenario where we have a Lipschitz constant $L_{i}$ for the partial derivative of $f$ with respect to each coordinate $i$,

and we use a coordinate-dependent step-size at each iteration:

By the logic of (2), in this setting we have

Noting that $L=\max_{i}\{L_{i}\}$, we have

Thus, the convergence rate based on the $L_{i}$ will be faster, provided that at least one iteration chooses an $i_{k}$ with $L_{i_{k}}<L$. In the worst case, however, (13) holds with equality even if the $L_{i}$ are distinct, as we might need to update a coordinate with $L_{i}=L$ on every iteration. (For example, consider a separable function where all but one coordinate is initialized at its optimal value, and the remaining coordinate has $L_{i}=L$.) In Section 6, we discuss selection rules that incorporate the $L_{i}$ to achieve faster rates whenever the $L_{i}$ are distinct, but first we consider the effect of exact coordinate optimization on the choice of the $L_{i_{k}}$.

For problems involving functions of the form $h_{1}$ and $h_{2}$, we are often able to perform exact (or numerically very precise) coordinate optimization, even if the objective function is not quadratic (e.g., by using a line-search or a closed-form update). Note that (12) still holds when using exact coordinate optimization rather than using a step-size of $1/L_{i_{k}}$, as in this case we have

which is equivalent to (11). However, in practice using exact coordinate optimization leads to better performance. In this section, we show that using the GS rule results in a convergence rate that is indeed faster than (9) for problems with distinct $L_{i}$ when the function is quadratic, or when the function is not quadratic but we perform exact coordinate optimization.

The key property we use is that, after we have performed exact coordinate optimization, we are guaranteed to have $\nabla_{i_{k}}f(x^{k+1})=0$. Because the GS rule chooses $i_{k+1}=\mathop{\hbox{argmax}}_{i}|\nabla_{i}f(x^{k+1})|$, we cannot have $i_{k+1}=i_{k}$, unless $x^{k+1}$ is the optimal solution. Hence, we never choose the same coordinate twice in a row, which guarantees that the inequality (13) is strict (with distinct $L_{i}$) and exact coordinate optimization is faster. We note that the improvement may be marginal, as we may simply alternate between the two largest $L_{i}$ values. However, consider minimizing $h_{2}$ when the graph is sparse; after updating $i_{k}$, we are guaranteed to have $\nabla_{i_{k}}f(x^{k+m})=0$ for all future iterations $(k+m)$ until we choose a variable $i_{k+m-1}$ that is a neighbour of node $i_{k}$ in the graph. Thus, if the two largest $L_{i}$ are not connected in the graph, GS cannot simply alternate between the two largest $L_{i}$.

By using this property, in Appendix D we show that the GS rule with exact coordinate optimization for problem $h_{2}$ under a chain-structured graph has a convergence rate of the form

where $\rho_{2}^{G}$ is the maximizer of $\sqrt{(1-\mu_{1}/L_{i})(1-\mu_{1}/L_{j})}$ among all consecutive nodes $i$ and $j$ in the chain, and $\rho_{3}^{G}$ is the maximizer of $\sqrt{(1-\mu_{1}/L_{i})(1-\mu_{1}/L_{j})(1-\mu_{1}/L_{k})}$ among consecutive nodes $i$, $j$, and $k$. The implication of this result is that, if the large $L_{i}$ values are more than two edges from each other in the graph, then we obtain a much better convergence rate. We conjecture that for general graphs, we can obtain a bound that depends on the largest value of $\rho_{2}^{G}$ among all nodes $i$ and $j$ connected by a path of length $1$ or $2$. Note that we can obtain similar results for problem $h_{1}$, by forming a graph that has an edge between nodes $i$ and $j$ whenever the corresponding variables are both jointly non-zero in at least one row of $A$.

Rules Depending on Lipschitz Constants

If the $L_{i}$ are known, Nesterov (2012) showed that we can obtain a faster convergence rate by sampling proportional to the $L_{i}$. We review this result below and compare it to the GS rule, and then propose an improved GS rule for this scenario. Although in this section we will assume that the $L_{i}$ are known, this assumption can be relaxed using a backtracking procedure (Nesterov, 2012, §6.1).

Taking the expectation of (11) under the distribution $p_{i}=L_{i}/\sum_{j=1}^{n}L_{j}$ and proceeding as before, we obtain

where $\bar{L}=\frac{1}{n}\sum_{j=1}^{n}L_{j}$ is the average of the Lipschitz constants. This was shown by Leventhal and Lewis (2010) and is a special case of Nesterov (2012, Theorem 2) with $\alpha=1$ in his notation. This rate is faster than (5) for uniform sampling if any $L_{i}$ differ.

Under our analysis, this rate may or may not be faster than (9) for the GS rule. On the one extreme, if $\mu_{1}=\mu/n$ and any $L_{i}$ differ, then this Lipschitz sampling scheme is faster than our rate for GS. Indeed, in the context of the problem from Section 4.1, we can make Lipschitz sampling faster than GS by a factor of nearly $n$ by making one $\lambda_{i}$ much larger than all the others (recall that our analysis shows no benefit to the GS rule over randomized selection when only one $\lambda_{i}$ is much larger than the others). At the other extreme, in our example from Section 4.1 with many large $\alpha$ and one small $\beta$, the GS and Lipschitz sampling rates are the same when $n=2$, with a rate of $(1-\beta/(\alpha+\beta))$. However, the GS rate will be faster than the Lipschitz sampling rate for any $\alpha>\beta$ when $n>2$, as the Lipschitz sampling rate is $(1-\beta/((n-1)\alpha+\beta))$, which is slower than the GS rate of $(1-\beta/(\alpha+(n-1)\beta))$.

2 Gauss-Southwell-Lipschitz Rule

Since neither Lipschitz sampling nor GS dominates the other in general, we are motivated to consider if faster rules are possible by combining the two approaches. Indeed, we obtain a faster rate by choosing the $i_{k}$ that minimizes (11), leading to the rule

which we call the Gauss-Southwell-Lipschitz (GSL) rule. Following a similar argument to Section 4, but using (11) in place of (2), the GSL rule obtains a convergence rate of

where $\mu_{L}$ is the strong-convexity constant with respect to the norm $\|x\|_{L}=\sum_{i=1}^{n}\sqrt{L_{i}}|x_{i}|$. This is shown in Appendix E, and in Appendix F we show that

Thus, the GSL rule is always at least as fast as the fastest of the GS rule and Lipschitz sampling. Indeed, it can be more than a factor of $n$ faster than using Lipschitz sampling, while it can obtain a rate closer to the minimum $L_{i}$, instead of the maximum $L_{i}$ that the classic GS rule depends on.

An interesting property of the GSL rule for quadratic functions is that it is the optimal myopic coordinate update. That is, if we have an oracle that can choose the coordinate and the step-size that decreases $f$ by the largest amount, i.e.,

this is equivalent to using the GSL rule and the update in (10). This follows because (11) holds with equality in the quadratic case, and the choice $\alpha_{k}=1/L_{i_{k}}$ yields the optimal step-size. Thus, although faster schemes could be possible with non-myopic strategies that cleverly choose the sequence of coordinates or step-sizes, if we can only perform one iteration, then the GSL rule cannot be improved.

For general $f$, (15) is known as the maximum improvement (MI) rule. This rule has been used in the context of boosting (Rätsch et al., 2001), graphical models (Della Pietra et al., 1997; Lee et al., 2006; Scheinberg and Rish, 2009), Gaussian processes (Bo and Sminchisescu, 2008), and low-rank tensor approximations (Li et al., 2015). Using an argument similar to (14), our GSL rate also applies to the MI rule, improving existing bounds on this strategy. However, the GSL rule is much cheaper and does not require any special structure (recall that we can estimate $L_{i}$ as we go).

3 Connection between GSL Rule and Normalized Nearest Neighbour Search

Dhillon et al. (2011) discuss an interesting connection between the GS rule and the nearest-neighbour-search (NNS) problem for objectives of the form

This is a special case of $h_{1}$ with no $g_{i}$ functions, and its gradient has the special form

where $r(x)=\nabla f(Ax)$. We use the symbol $r$ because it is the residual vector ($Ax-b$) in the special case of least squares. For this problem structure the GS rule has the form

where $a_{i}$ denotes column $i$ of $A$ for $i=1,\dots,n$. Dhillon et al. (2011) propose to approximate the above $\mathop{\hbox{argmax}}$ by solving the following NNS problem

where $i$ in the range $(n+1)$ through $2n$ refers to the negation $-(a_{i-n})$ of column $(i-n)$ and if the selected $i_{k}$ is greater than $n$ we return $(i-n)$. We can justify this approximation using the logic

Thus, the NNS computes an approximation to the GS rule that is biased towards coordinates where $\|a_{i}\|$ is small. Note that this formulation is equivalent to the GS rule in the special case that $\|a_{i}\|=1$ (or any other constant) for all $i$. Shrivastava and Li (2014) have more recently considered the case where $\|a_{i}\|\leq 1$ and incorporate powers of $\|a_{i}\|$ in the NNS to yield a better approximation.

A further interesting property of the GSL rule is that we can often formulate the exact GSL rule as a normalized NNS problem. In particular, for problem (16) the Lipschitz constants will often have the form $L_{i}=\gamma\|a_{i}\|^{2}$ for a some positive scalar $\gamma$. For example, least squares has $\gamma=1$ and logistic regression has $\gamma=0.25$. When the Lipschitz constants have this form, we can compute the exact GSL rule by solving a normalized NNS problem,

The exactness of this formula follows because

Thus, the form of the Lipschitz constant conveniently removes the bias towards smaller values of $\|a_{i}\|$ that gets introduced when we try to formulate the classic GS rule as a NNS problem. Interestingly, in this setting we do not need to know $\gamma$ to implement the GSL rule as a NNS problem.

Approximate Gauss-Southwell

In many applications, computing the exact GS rule is too inefficient to be of any practical use. However, a computationally cheaper approximate GS rule might be available. Approximate GS rules under multiplicative and additive errors were considered by Dhillon et al. (2011) in the convex case, but in this setting the convergence rate is similar to the rate achieved by random selection. In this section, we give rates depending on $\mu_{1}$ for approximate GS rules.

In the multiplicative error regime, the approximate GS rule chooses an $i_{k}$ satisfying

for some $\epsilon_{k}\in[0,1)$. In this regime, our basic bound on the progress (2) still holds, as it was defined for any $i_{k}$. We can incorporate this type of error into our lower bound (8) to obtain

Thus, the convergence rate of the method is nearly identical to using the exact GS rule for small $\epsilon_{k}$ (and it degrades gracefully with $\epsilon_{k})$. This is in contrast to having an error in the gradient (Friedlander and Schmidt, 2012), where the error $\epsilon$ must decrease to zero over time.

2 Additive Errors

In the additive error regime, the approximate GS rule chooses an $i_{k}$ satisfying

for some $\epsilon_{k}\geq 0$. In Appendix G, we show that under this rule, we have

where $L_{1}$ is the Lipschitz constant of $\nabla f$ with respect to the 1-norm. Note that $L_{1}$ could be substantially larger than $L$, so the second part of the maximum in $A_{k}$ is likely to be the smaller part unless the $\epsilon_{i}$ are large. This regime is closer to the case of having an error in the gradient, as to obtain convergence the $\epsilon_{k}$ must decrease to zero. This result implies that a sufficient condition for the algorithm to obtain a linear convergence rate is that the errors $\epsilon_{k}$ converge to zero at a linear rate. Further, if the errors satisfy $\epsilon_{k}=O(\rho^{k})$ for some $\rho<(1-\mu_{1}/L)$, then the convergence rate of the method is the same as if we used an exact GS rule. On the other hand, if $\epsilon_{k}$ does not decrease to zero, we may end up repeatedly updating the same wrong coordinate and the algorithm will not converge (though we could switch to the randomized method if this is detected).

Proximal-Gradient Gauss-Southwell

One of the key motivations for the resurgence of interest in coordinate descent methods is their performance on problems of the form

With random coordinate selection, Richtárik and Takáč (2014) show that this method has a convergence rate of

similar to the unconstrained/smooth case.

However, the length of the step ($\|x^{k+1}-x^{k}\|$) could be arbitrarily small under this choice. In contrast, the GS-$r$ rule chooses the coordinate that maximizes the length of the step (Tseng and Yun, 2009; Dhillon et al., 2011),

This rule is effective for bound-constrained problems, but it ignores the change in the non-smooth term ($g_{i}(x_{i}^{k+1})-g_{i}(x_{k}^{k})$). Finally, the GS-$q$ rule maximizes progress assuming a quadratic upper bound on $f$ (Tseng and Yun, 2009),

While the least intuitive rule, the GS-$q$ rule seems to have the best theoretical properties. Further, if we use $L_{i}$ in place of $L$ in the GS-$q$ rule (which we call the GSL-$q$ strategy), then we obtain the GSL rule if the $g_{i}$ are not present. In contrast, using $L_{i}$ in place of $L$ in the GS-$r$ rule (which we call the GSL-$r$ strategy) does not yield the GSL rule as a special case.

In Appendix H, we show that using the GS-$q$ rule yields a convergence rate of

thus matching the convergence rate of randomized coordinate descent (but deterministically rather than in expectation). In contrast, in Appendix H we also give counter-examples showing that the above rate does not hold for the GS-$s$ or the GS-$r$ rule. Thus, any bound for the GS-$s$ or the GS-$r$ rule would be slower than the expected rate under random selection, while the GS-$q$ rule matches this bound. It is an open problem whether the GS-$q$ rule obtains the rate $(1-\mu_{1}/L)$ in general, but in the next section we discuss special cases where rates depending on $\mu_{1}$ can be obtained.

First, we note that if the $g_{i}$ are linear then the GS-$q$ rule obtains

since in this particular (smooth) case the algorithm and assumptions are identical to the setting of Section 4: the GS-q rule chooses the same coordinate to update as the GS rule applied to $F$, while $F$ has the same $L$ and $\mu_{1}$ as $f$ because the $g_{i}$ are linear.

Finally, consider the case of $g_{i}$ that are piecewise-linear. Under a suitable non-degeneracy assumption, coordinate descent methods achieve a particular “active set” property in a finite number of iterations (Wright, 2012; Nutini et al., 2017). Specifically, for values of $x_{i}^{*}$ that occur at non-smooth values of $g_{i}$, we will have $x_{i}^{k}=x_{i}^{*}$ for all sufficiently large $k$. At this point, none of the four GS-$*$ rules would select such coordinates again. Similarly, for values where $x_{i}^{*}$ occurs at smooth values of $g_{i}$, the iterates will eventually be confined to a region where the $g_{i}$ is smooth. Once this “active set” identification happens for piecewise-linear $g_{i}$, the iterates will be confined to a region where the selected $g_{i}$ are linear. At this point, linearity means that only one coordinate will be selected by the GS-$1$ rule and it will select the same coordinate as the GS-$q$ rule. Further, at this point the analysis of Nutini (2018, Appendix A.8) can be applied with $L$ instead $L_{1}$ for the GS-$q$ rule which leads to a rate of $(1-\mu_{1}/L)$ as in the smooth case.

Experiments

an instance of problem $h_{1}$. We set $A$ to be an $m$ by $n$ matrix with entries sampled from a $\mathcal{N}(0,1)$ distribution (with $m=1000$ and $n=1000$). We then added 1 to each entry (to induce a dependency between columns), multiplied each column by a sample from $\mathcal{N}(0,1)$ multiplied by ten (to induce different Lipschitz constants across the coordinates), and only kept each entry of $A$ non-zero with probability $10\log(n)/n$ (a sparsity level that allows the Gauss-Southwell rule to be applied with cost $O(\log^{3}(n))$. We set $\lambda=1$ and $b=Ax+e$, where the entries of $x$ and $e$ were drawn from a $\mathcal{N}(0,1)$ distribution. In this setting, we used a step-size of $1/L_{i}$ for each coordinate $i$, which corresponds to exact coordinate optimization.

We set the $a_{i}^{T}$ to be the rows of $A$ from the previous problem, and set $b=$ sign$(Ax)$, but randomly flipping each $b_{i}$ with probability $0.1$. In this setting, we compared using a step-size of $1/L_{i}$ to using exact coordinate optimization.

Over-determined dense least squares: Here we consider the problem

but, unlike the previous case, we do not set elements of $A$ to zero and we make $A$ have dimension $1000$ by $100$. Because the system is over-determined, it does not need an explicit strongly-convex regularizer to induce global strong-convexity. In this case, the density level means that the exact GS rule is not efficient. Hence, we use a balltree structure (Omohundro, 1989) to implement an efficient approximate GS rule based on the connection to the NNS problem discovered by Dhillon et al. (2011). On the other hand, we can compute the exact GSL rule for this problem as a NNS problem as discussed in Section 6.3.

We next consider an instance of problem $h_{2}$, performing label propagation for semi-supervised learning in the ‘two moons’ dataset (Zhou et al., 2004). We generate $500$ samples from this dataset, randomly label five points in the data, and connect each node to its five nearest neighbours. This high level of sparsity is typical of graph-based methods for semi-supervised learning, and allows the exact Gauss-Southwell rule to be implemented efficiently. We use the quadratic labeling criterion of Bengio et al. (2006), which allows exact coordinate optimization and is normally optimized with cyclic coordinate descent. We plot the performance under different selection rules in Figure 2. Here, we see that even cyclic coordinate descent outperforms randomized coordinate descent, but that the GS and GSL rules give even better performance. We note that the GS and GSL rules perform similarly on this problem since the Lipschitz constants do not vary much.

Discussion

It is clear that the GS rule is not practical for every problem where randomized methods are applicable. Nevertheless, we have shown that even approximate GS rules can obtain better convergence rate bounds than fully-randomized methods. We have given a similar justification for the use of exact coordinate optimization, and we note that our argument could also be used to justify the use of exact coordinate optimization within randomized coordinate descent methods (as used in our experiments). We have also proposed the improved GSL rule, and considered approximate/proximal variants. We expect our analysis also applies to block updates by using mixed norms $\|\cdot\|_{p,q}$, and could be used for accelerated/parallel methods (Fercoq and Richtárik, 2013), for primal-dual rates of dual coordinate ascent (Shalev-Shwartz and Zhang, 2013), for successive projection methods (Leventhal and Lewis, 2010), for boosting algorithms (Rätsch et al., 2001), and for scenarios without strong-convexity under general error bounds (Luo and Tseng, 1993).

Acknowledgements

We would like to thank the anonymous referees for their useful comments that significantly improved the paper. Julie Nutini is funded by an NSERC Canada Graduate Scholarship.

Appendix A Efficient calculation of GS rules for sparse problems

We first give additional details on how to calculate the GS rule efficiently for sparse instances of problems $h_{1}$ and $h_{2}$. We will consider the case where each $g_{i}$ is smooth, but the ideas can be extended to allow a non-smooth $g_{i}$. Further, note that the efficient calculation does not rely on convexity, so these strategies can also be used for non-convex problems.

where each $g_{i}$ and $f_{ij}$ are differentiable and $G=\{V,E\}$ is a graph where the number of vertices $|V|$ is the same as the number of variables $n$. If all nodes in the graph have a degree (number of neighbours) bounded above by some constant $d$, we can implement the GS rule in $O(d\log n)$ after an $O(n+|E|)$ time initialization by maintaining the following information about $x^{k}$:

A vector containing the values $\nabla_{i}g_{i}(x_{i}^{k})$.

A matrix containing the values $\nabla_{i}f_{ij}(x_{i}^{k},x_{j}^{k})$ in the first column and $\nabla_{j}f_{ij}(x_{i}^{k},x_{j}^{k})$ in the second column.

The elements of the gradient vector $\nabla h_{2}(x^{k})$ stored in a binary max heap data structure (see Cormen et al., 2001, Chapter 6).

Given the heap structure, we can compute the GS rule in $O(1)$ by simply reading the index value of the root node in the max heap. The costs for initializing these structures are:

$O(n)$ to compute $g_{i}(x_{i}^{0})$ for all $n$ nodes.

$O(|E|)$ to compute $\nabla_{ij}f_{ij}(x_{i}^{0},x_{j}^{0})$ for all $|E|$ edges.

$O(n+|E|)$ to sum the values in the above structures to compute $\nabla h(x^{0})$, and $O(n)$ to construct the initial max heap.

Thus, the one-time initialization cost is $O(n+|E|)$. The costs of updating the data structures after we update $x_{i_{k}}^{k}$ to $x_{i_{k}}^{k+1}$ for the selected coordinate $i_{k}$ are:

$O(1)$ to compute $g_{i_{k}}(x_{i_{k}}^{k+1})$.

$O(d)$ to compute $\nabla_{ij}f_{ij}(x_{i}^{k+1},x_{j}^{k+1})$ for $(i,j)\in E$ and $i=i_{k}$ or $j=i_{k}$ (only $d$ such values exist by assumption, and all other $\nabla_{ij}f_{ij}(x_{i},x_{j})$ are unchanged).

$O(d)$ to update up to $d$ elements of $\nabla h(x^{k+1})$ that differ from $\nabla h(x^{k})$ by using differences in changed values of $g_{i}$ and $f_{ij}$, followed by $O(d\log n)$ to perform $d$ updates of the heap at a cost of $O(\log n)$ for each update.

The most expensive part of the update is modifying the heap, and thus the total cost is $O(d\log n)$.For less-sparse problems where $n<d\log n$, using a heap is actually inefficient and we should simply store $\nabla h(x^{k})$ as a vector. The initialization cost is the same, but we can then perform the GS rule in $O(n)$ by simply searching through the vector for the maximum element.

where $g_{i}$ and $f$ are differentiable, and $A$ is an $m$ by $n$ matrix where we denote column $i$ by $a_{i}$ and row $j$ by $a_{j}^{T}$. Note that $f$ is a function from ${\rm I\!R}^{m}$ to ${\rm I\!R}$, and we assume $\nabla_{j}f$ only depends on $a_{j}^{T}x$. While this is a strong assumption (e.g., it rules out $f$ being the product function), this class includes a variety of notable problems like the least squares and logistic regression models from our experiments. If $A$ has $z$ non-zero elements, with a maximum of $c$ non-zero elements in each column and $r$ non-zero elements in each row, then with a pre-processing cost of $O(z)$ we can implement the GS rule in this setting in $O(cr\log n)$ by maintaining the following information about $x^{k}$:

A vector containing the values $\nabla_{i}g_{i}(x_{i}^{k})$.

A vector containing the product $Ax^{k}$.

A vector containing the values $\nabla f(Ax^{k})$.

A vector containing the product $A^{T}\nabla f(Ax^{k})$.

The elements of the gradient vector $\nabla h_{1}(x^{k})$ stored in a binary max heap data structure.

The heap structure again allows us to compute the GS rule in $O(1)$, and the costs of initializing these structures are:

$O(n)$ to compute $g_{i}(x_{i}^{0})$ for all $n$ variables.

$O(m)$ to compute $\nabla f(Ax^{0})$ (using that $\nabla_{j}f$ only depends on $a_{j}^{T}x^{0}$).

$O(z)$ to compute $A^{T}\nabla f(Ax^{0})$.

$O(n)$ to add the $\nabla_{i}g_{i}(x_{i}^{0})$ to the above product to obtain $\nabla h_{1}(x^{0})$ and construct the initial max heap.

As it is reasonable to assume that $z\geq m$ and $z\geq n$ (e.g., we have at least one non-zero in each row and column), the cost of the initialization is thus $O(z)$. The costs of updating the data structures after we update $x_{i_{k}}^{k}$ to $x_{i_{k}}^{k+1}$ for the selected coordinate $i_{k}$ are:

$O(1)$ to compute $g_{i_{k}}(x_{i_{k}}^{k+1})$.

$O(c)$ to update the product using $Ax^{k+1}=Ax^{k}+(x_{i_{k}}^{k+1}-x_{i_{k}}^{k})a_{i}$, since $a_{i}$ has at most $c$ non-zero values.

$O(c)$ to update up to $c$ elements of $\nabla f(Ax^{k+1})$ that have changed (again using that $\nabla_{j}f$ only depends on $a_{j}^{T}x^{k+1}$).

$O(cr)$ to perform up to $c$ updates of the form $A^{T}\nabla f(Ax^{k+1})=A^{T}\nabla f(Ax^{k})+(\nabla_{j}f(Ax^{k+1})-\nabla_{j}f(Ax^{k}))(a_{i})^{T}$, where each update costs $O(r)$ since each $a_{i}$ has at most $r$ non-zero values.

$O(cr\log n)$ to update the gradients in the heap.

The most expensive part is again the heap update, and thus the total cost is $O(cr\log n)$.

We can establish the relationship between $\mu$ and $\mu_{1}$ by using the known relationship between the $2$-norm and the $1$-norm,

In particular, if we assume that $f$ is $\mu$-strongly convex in the $2$-norm, then for all $x$ and $y$ we have

implying that $f$ is at least $\frac{\mu}{n}$-strongly convex in the $1$-norm. Similarly, if we assume that a given $f$ is $\mu_{1}$-strongly convex in the $1$-norm then for all $x$ and $y$ we have

implying that $f$ is at least $\mu_{1}$-strongly convex in the $2$-norm. Summarizing these two relationships, we have

Appendix C Analysis for separable quadratic case

We first establish an equivalent definition of strong-convexity in the $1$-norm, along the lines of Nesterov (2004, Theorem 2.1.9). Subsequently, we use this equivalent definition to derive $\mu_{1}$ for a separable quadratic function.

Assume that $f$ is $\mu_{1}$-strongly convex in the $1$-norm, so that for any $x,y\in{\rm I\!R}^{n}$ we have

Conversely, assume that for all $x$ and $y$ we have

and consider the function $g(\tau)=f(x+\tau(y-x))$ for $\tau\in{\rm I\!R}$. Then

Thus, $\mu_{1}$-strong convexity in the $1$-norm is equivalent to having

Consider a strongly convex quadratic function $f$ with a diagonal Hessian $H=\nabla^{2}f(x)=\hbox{diag}({\lambda_{1},\dots,\lambda_{n}})$, where $\lambda_{i}>0$ for all $i=1,\dots,n$. We show that in this case

From the previous section, $\mu_{1}$ is the minimum value such that (20) holds,

Using $\nabla f(x)=Hx+b$ for some $b$ and letting $z=y-x$, we get

where the last two lines use that the objective is invariant to scaling of $z$ and to the sign of $z$ (respectively), and where $e$ is a vector containing a one in every position. This is an equality-constrained strictly-convex quadratic program, so its solution is given as a stationary point $(z^{*},\eta^{*})$ of the Lagrangian,

Differentiating with respect to each $z_{i}$ for $i=1,\dots,n$ and equating to zero, we have for all $i$ that $2\lambda_{i}z_{i}^{*}-\eta^{*}=0$, or

Differentiating the Lagrangian with respect to $\eta$ and equating to zero we obtain $1-e^{T}z^{*}=0$, or equivalently

Combining this result for $\eta^{*}$ with equation (21), we have

This gives the minimizer, so we evaluate the objective at this point to obtain $\mu_{1}$,

Appendix D Gauss-Southwell with exact optimization

We can obtain a faster convergence for GS using exact coordinate optimization for sparse variants of problems $h_{1}$ and $h_{2}$, by observing that the convergence rate can be expressed in terms of the sequence of $(1-\mu_{1}/L_{i_{k}})$ values,

The worst case occurs when the product of the $(1-\mu_{1}/L_{i_{k}})$ values is as large as possible. However, using exact coordinate optimization guarantees that, after we have updated coordinate $i$, the GS rule will never select it again until one of its neighbours has been selected. Thus, we can obtain a tighter bound on the worst-case convergence rate using GS with exact coordinate optimization on iteration $k$, by solving the following combinatorial optimization problem defined on a weighted graph:

We are given a graph $G=(V,E)$ with $n$ nodes, a number $M_{i}$ associated with each node $i$, and an iteration number $k$. Choose a sequence $\{i_{t}\}_{t=1}^{k}$ that maximizes the sum of the $M_{i_{t}}$, subject to the following constraint: after each time node $i$ has been chosen, it cannot be chosen again until after a neighbour of node $i$ has been chosen.

We can use the $M_{i}$ chosen by this problem to obtain an upper-bound on the sequence of $\log(1-\mu_{1}/L_{i})$ values, and if the largest $M_{i}$ values are not close to each other in the graph, then this rate can be much faster than the rate obtained by alternating between the largest $M_{i}$ values. In the particular case of chain-structured graphs, a worst-case sequence can be constructed that spends all but $O(n)$ iterations in one of two solution modes: (i) alternate between two nodes $i$ and $j$ that are connected by an edge with the highest value of $\frac{M_{i}+M_{j}}{2}$, or (ii) alternate between three nodes $\{i,j,k\}$ with the highest value of $\frac{M_{i}+M_{j}+M_{k}}{3}$, where there is an edge from $i$ to $j$ and from $j$ to $k$, but not from $i$ to $k$. To show that these are the two solution modes, observe that the solution must eventually cycle because there are a finite number of nodes. If you have more than three nodes in the cycle, then you can always remove one node from the cycle to obtain a better average weight for the cycle without violating the constraint. We will fall into mode (i) if the average of $M_{i}$ and $M_{j}$ in this mode is larger than the average of $M_{i}$, $M_{j}$ and $M_{k}$ in the second mode. We can construct a solution to this problem that consists of a ‘burn-in’ period, where we choose the largest $M_{i}$, followed by repeatedly going through the better of the two solution modes up until the final three steps, where a ‘burn-out’ phase arranges to finish with several large $M_{i}$. By setting $M_{i}=\log(1-\mu_{1}/L_{i})$, this leads to a convergence rate of the form

where $\rho_{2}^{G}$ is the maximizer of $\sqrt{(1-\mu_{1}/L_{i})(1-\mu_{1}/L_{j})}$ among all consecutive nodes $i$ and $j$ in the chain, and $\rho_{3}^{G}$ is the maximizer of $\sqrt{(1-\mu_{1}/L_{i})(1-\mu_{1}/L_{j})(1-\mu_{1}/L_{k})}$ among consecutive nodes $i$, $j$, and $k$. The $O()$ notation gives the constant due to choosing higher $(1-\mu_{1}/L_{i})$ values during the burn-in and burn-out periods. The implication of this result is that, if the large $L_{i}$ values are more than two edges away from each other in the graph, then the convergence rate can be much faster.

Appendix E Gauss-Southwell-Lipschitz rule: convergence rate

The coordinate-descent method with a constant step-size of $L_{i_{k}}$ uses the iteration

Because $f$ is coordinate-wise $L_{i_{k}}$-Lipschitz continuous, we obtain the following bound on the progress made by each iteration:

By choosing the coordinate to update according to the Gauss-Southwell-Lipchitz (GSL) rule,

we obtain the tightest possible bound on (22). We define the following norm,

Under this notation, and using the GSL rule, (22) becomes

Measuring strong-convexity in the norm $\|\cdot\|_{L}$ we get

Minimizing both sides with respect to $y$ we get

By the logic Appendix B, to establish a relationship between different strong-convexity constants under different norms, it is sufficient to establish the relationships between the squared norms. In this section, we use this to establish the relationship between $\mu_{L}$ defined in (23) and both $\mu_{1}$ and $\mu$.

Assuming $c\geq\sqrt{L}$, where $L=\max_{i}\{L_{i}\}$, the expression is non-negative and we get

and assuming $\displaystyle c\geq\frac{1}{\sqrt{L_{min}}}$, where $L_{min}=\min_{i}\{L_{i}\}$, this expression is nonnegative and we get

The relationship between $\mu_{L}$ and $\mu_{1}$ is based on the squared norm, so in summary we have

Let $\vec{L}$ denote a vector with elements $\sqrt{L_{i}}$, and we note that

Note that we can also show that $\mu_{L}\leq\frac{\mu}{L_{min}}$, but this is less tight than the upper bound from the previous section because $\mu_{1}\leq\mu$.

Appendix G Approximate Gauss-Southwell with additive error

In the additive error regime, the approximate Gauss-Southwell rule chooses an $i_{k}$ satisfying

and we note that we can assume $\epsilon_{k}\leq\|\nabla f(x^{k})\|_{\infty}$ without loss of generality because we must always choose an $i$ with $|\nabla_{i_{k}}f(x^{k})|\geq 0$. Applying this to our bound on the iteration progress, we get

We first give a result that assumes $f$ is $L_{1}$-Lipschitz continuous in the $1$-norm. This implies an inequality that we prove next, followed by a convergence rate that depends on $L_{1}$. However, note that $L\leq L_{1}\leq Ln$, so this potentially introduces a dependency on $n$. We subsequently give a slightly less concise result that has a worse dependency on $\epsilon$ but does not rely on $L_{1}$.

We say that $\nabla f$ is $L_{1}$-Lipschitz continuous in the $1$-norm if we have for all $x$ and $y$ that

Similar to Nesterov (2004, Theorem 2.1.5), we now show that this implies

where we have used that $f(x^{k})\leq f(x^{k-1})$ for all $k$ and any choice of $i_{k-1}$ (this follows from the basic bound on the progress of coordinate descent methods).

We first show that $\nabla f$ being $L_{1}$-Lipschitz continuous in the $1$-norm implies that

for all $x$ and $y$. Consider the function $g(\tau)=f(x+\tau(y-x))$ with $\tau\in{\rm I\!R}$. Then

To subsequently show (26), fix $x\in{\rm I\!R}^{n}$ and consider the function

which is convex on ${\rm I\!R}^{n}$ and also has an $L_{1}$-Lipschitz continuous gradient in the $1$-norm, as

As the minimizer of $\phi$ is $x$ (i.e., $\phi^{\prime}(x)=0$), for any $y\in{\rm I\!R}^{n}$ we have

Substituting in the definition of $\phi$, we have

Using (27) in (25) and noting that $\epsilon_{k}\geq 0$, we obtain

Applying strong convexity (taken with respect to the $1$-norm), we get

G.3 Additive error bound in terms of L𝐿L

By our additive error inequality, we have

Further, from our basic progress bound that holds for any $i_{k}$ we have

Applying strong convexity and applying the inequality recursively we obtain

Although uglier than the expression depending on $L_{1}$, this expression will tend to be smaller unless $\epsilon_{k}$ is not small.

Appendix H Convergence Analysis of GS-s𝑠s, GS-r𝑟r, and GS-q𝑞q Rules

In this section, we consider problems of the form

where $f$ satisfies our usual assumptions, but the $g_{i}$ can be non-smooth. We first introduce some notation and state the convergence result, and then show that it holds. We then show that the rate cannot hold in general for the GS-$s$ and GS-$r$ rules.

To analyze this case, an important inequality we will use is that the $L$-Lipschitz-continuity of $\nabla_{i}f$ implies that for all $x$, $i$, and $d$, we have

Notice that the GS-$q$ rule is defined by

We use the notation $d_{i}^{k}=\mathop{\hbox{argmin}}_{d}V_{i}(x^{k},d)$ and we will use $d^{k}$ to denote the vector containing these values for all $i$. When using the GS-$q$ rule, the iteration is defined by

In this notation the GS-$r$ rule is given by

Under this notation, we can show that coordinate descent with the GS-$q$ rule satisfies the bound

We show this result by showing that the GS-$q$ rule makes at least as much progress as randomized selection.

H.2 GS-q𝑞q is at least as fast as random

Our argument in this section follows a similar approach to Richtárik and Takáč (2014). In particular, combining (28) and (29) we have the following upper bound on the iteration progress

From strong convexity of $f$, we have that $F$ is also $\mu$-strongly convex and that

for any $y\in{\rm I\!R}^{n}$ and any $\alpha\in$ (see Nesterov, 2004, Theorem 2.1.9). Using these gives us

Subtracting $F(x^{*})$ from both sides of this inequality gives us

H.3 Lack of progress of the GS-s𝑠s rule

We now show that the rate $(1-\mu/Ln)$ cannot hold for the GS-$s$ rule. We do this by constructing a problem where an iteration of the GS-$s$ method does not make sufficient progress. In particular, consider the bound-constrained problem

The parameter values for this problem are

where the $\lambda_{i}$ are the eigenvalues of $A^{T}A$, and $\mu$ and $\mu_{1}$ are the corresponding strong-convexity constants for the $2$-norm and $1$-norm, respectively.

The proximal operator of the indicator function is the projection onto the set $C$, which involves setting negative elements to zero. Thus, our iteration update is given by

For this problem, the GS-$s$ rule is given by

Based on the value of $\nabla f(x^{0})$, the GS-$s$ rule thus chooses to update coordinate 2, setting it to zero and obtaining

Thus, the GS-$s$ rule does not satisfy either bound. On the other hand, the GS-$r$ and GS-$q$ rules are given in this context by

and thus both these rules choose to update coordinate 1, setting it to zero to obtain $f(x^{1})\approx 5.2$ and a progress ratio of

H.4 Lack of progress of the GS-r𝑟r rule

We use the same $A$ as the previous section, so that $n$, $\mu$, $L$, and $\mu_{1}$ are the same. However, we now take

The proximal operator of the absolute value function is given by the soft-threshold function, and our coordinate update of variable $i_{k}$ is given by

where $d_{i}^{k}=\mathop{\rm prox}_{\lambda|\cdot|}[x_{i}^{k+\frac{1}{2}}]-x_{i}^{k}$ and in this case

Thus, the GS-$r$ rule chooses to update coordinate $1$. After this update the function value is

However, the bounds suggest faster progress ratios of

so the GS-$r$ rule does not satisfy either bound. In contrast, in this setting the GS-$q$ rule chooses to update coordinate $2$ and obtains $F(x^{1})\approx 2.2$, obtaining a progress ratio of

which satisfies both bounds by a substantial margin. Indeed, we used a genetic algorithm to search for a setting of the parameters of this problem (values of $x^{0}$, $\lambda$, $b$, and the diagonals of $A$) that would make the GS-$q$ rule not satisfy the bound depending on $\mu_{1}$, and it easily found counter-examples for the GS-$s$ and GS-$r$ rules but was not able to produce a counter example for the GS-$q$ rule.

Appendix I Runtime Experiments

References