Random Permutations Fix a Worst Case for Cyclic Coordinate Descent

Introduction

The basic (component-wise) coordinate descent framework for the smooth unconstrained optimization problem

When $f$ is a convex quadratic function, and when $\alpha_{k}$ in Algorithm 1 is chosen to minimize $f$ exactly along each coordinate direction, these variants are simply different variants of the Gauss-Seidel approach for solving the equivalent system of linear equations.

The coordinate descent approach is enjoying renewed popularity because of its usefulness in data analysis applications. Its convergence properties have come under renewed scrutiny. We refer to (Wright, 2015b) for a discussion of the state of the art as of 2015, but make a few additions and updates here, with an emphasis on results concerning linear convergence of the function values, by which we mean epoch-wise convergence of the form

where $\rho$ is typically much closer to $1$ than to , and $f^{*}$ is the optimal value of (1). For randomized methods, we consider a corresponding expression in expectation:

where the expectation is taken over all random variables encountered in the algorithm. When (3) holds, a reduction in function error by a factor of $\epsilon$ can be attained in approximately $|\log\epsilon|/(1-\rho)$ epochs. We sometimes refer to $1/(1-\rho)$ as the “complexity” of an algorithm for which (3) or (4) holds.

The standard Lipschitz constant $L$ is defined so that

(Here and throughout, we use $\|\cdot\|$ to denote $\|\cdot\|_{2}$.) For reasonable choices of the constants in (5), (6), and (7), the following bounds are satisfied:

The following property of Łojasiewicz (1963) is useful in proving linear convergence:

This property holds for $f$ strongly convex (with modulus of strong convexity $\mu$), and for the case in which $f$ grows quadratically with distance from a non-unique minimizing set, as in the “optimal strong convexity” condition of Liu & Wright (2015, (1.2)). It also holds generically for convex quadratic programs, even when the Hessians are singular. Further, condition (9) holds for the functional form considered by Luo & Tseng (1992, 1993), which is

without any conditions on $E$. (For a proof, see Appendix C.) In (Karimi et al., 2016), property (9) is called the “Polyak-Łojasiewicz condition.”

In this paper, our focus is on the case of $f$ convex quadratic, that is,

For this function, the values of $L_{i}$, $L$, $L_{\max}$, and $\mu$ are as follows:

where $\lambda_{\text{min,nz}}(\cdot)$ denotes the minimum nonzero eigenvalue. For such functions, the upper bound in (8) is achieved by $A=\mathbf{1}\mathbf{1}^{T}$ (where $\mathbf{1}=(1,1,\dotsc,1)^{T}$), for which $L_{i}=1$, $i=1,2,\dotsc,n$; $L_{\max}=1$; and $L=n$.

We have not included a linear term in (11), but note that there is no loss of generality in doing so. If we were to consider instead

(note that $x^{*}$ is the minimizer of this function), the main results of Sections 2 and 3 would continue to hold, except that in several theorems the initial iterate $x^{0}$ would be replaced by $x^{0}-x^{*}$, and $f(x)$ is replaced by $f(x)-f(x^{*})$.

2 Linear Convergence Results for CD Variants

Luo & Tseng (1992) prove linear convergence for a function of the form (10), where they require $E$ to have no zero columns. They obtain expressions for the constant $\rho$ in (3) for two variants of CD — a Gauss-Southwell variant and an “almost cyclic” rule — but these constants are difficult to characterize in terms of fundamental properties of $f$. In (Luo & Tseng, 1993), the same authors analyze a family of methods (including CD) for more general functions that satisfy a local error bound of the form $\|x-P(x)\|\leq\chi\|\nabla f(x)\|$ holds (where $P(x)$ is the projection of $x$ onto the solution set of (1) and $\chi$ is some constant). Again, their analysis is not clear about how the constant $\rho$ of (3) depends on the properties of $f$.

A family of linear convergence results is proved in Beck & Tetruashvili (2013, Theorem 3.9) for the case in which $f$ is strongly convex (immediately extendable to the case in which $f$ satisfies the condition (9)). For constant stepsizes $\alpha_{k}\equiv\alpha\leq 1/L_{\max}$, convergence of the form (3) holds with

For the random-permutations cyclic version RPCD, the convergence theory in (Beck & Tetruashvili, 2013) can be applied without modification to attain the bounds given above. As we discuss below, however, the practical performance of RPCD is sometimes much better than these bounds would suggest.

A different convergence rate is proved in Nesterov (2012, Theorem 5), namely,

for some constant $C$ depending on the initial point. This is an R-linear expression, obtained from Q-linear convergence of the modified function $f(x)-f(x^{*})+\sum_{i}L_{i}(x_{i}-x^{*}_{i})^{2}/2$, where $x^{*}$ is the (unique) solution of (1). It indicates a complexity of approximately $|\log\epsilon|(L_{\max}+\mu)/(2\mu)$.

An important benchmark in studying the convergence rates of coordinate descent is the steepest-descent (SD) method, which takes a step from $x^{k}$ along all coordinates simultaneously, in the direction $-\nabla f(x^{k})$. For some important types of functions, including empirical-risk-minimization functions that arise in data analysis, the computational cost of one steepest-descent step is comparable to the cost of one epoch of Algorithm 1 (see (Wright, 2015b)). Standard analysis of steepest descent shows that fixed-steplength variants applied to functions satisfying (9) have linear convergence of the form (3) (with one iteration of SD replacing one epoch of Algorithm 1) with $\rho=1-\mu/L$. This worst-case complexity is not improved qualitatively by using exact line searches.

In comparing convergence rates between CCD and SD (on the one hand) and RCD (on the other hand), we see that the former tend to depend on $L$ while the latter depends on $L_{\max}$. These bounds suggest that CCD may tend to track the performance of SD, while RCD could be significantly better if the ratio $L/L_{\max}$ is large, that is, toward the upper end of its range in (8). The phenomenon of large values of $L/L_{\max}$ is captured well by convex quadratic problems (11) in which the Hessian $A$ has a large contribution from $\mathbf{1}\mathbf{1}^{T}$. Such matrices were used in computations by one of the authors in 2015 (see (Wright, 2015a); reported briefly in (Wright, 2015b)). These tests showed that on such matrices, RCD was indeed much faster than CCD (and also SD). The performance of RPCD was as fast as that of RCD; it did not track CCD as the obvious worst-case analysis would suggest. Later work, reported in (Wright, 2015c), identified the matrix

3 Motivation and Outline

Our focus in this paper is to analyze the performance of RPCD for minimizing (11) with $A$ defined in (17). Our interest in RPCD is motivated by computational practice. Much has been written about randomized optimization algorithms (particularly stochastic gradient and coordinate descent) in recent years. The analysis usually applies to sampling-with-replacement versions, but the implementations almost always involve a sampling-without-replacement scheme. The reasons are clear: Convergence analysis is much more straightforward for sampling with replacement, while for sampling without replacement, implementations are more efficient, involving less data movement. Moreover, it has long been folklore in the machine learning community that sampling-without-replacement schemes perform better in practice. In this paper we take a step toward bringing the analysis into line with the practice, by giving a tight analysis of the sampling-without-replacement scheme RPCD, on a special but important function that captures perfectly the advantages of randomized schemes over a deterministic scheme.

In Section 2, we derive tools for analyzing epoch-wise convergence of CD variants on convex quadratic problems (11), focusing on the permutation-invariant matrix (17) and recalling results for the CCD and RCD variant in this case (obtained from (Sun & Ye, 2016) and (Nesterov, 2012)). Section 3 contains our results for RPCD applied to (11) with the permutation-invariant matrix (17), characterizing its convergence rate in terms of a two-parameter recurrence. The relationship of this two-parameter sequence to the expected function value at the end of each epoch is described in Theorem 3.5. Our main result, Theorem 3.7, gives bounds on these two parameters in terms of $\delta$ (the parameter that defines (17)) and epoch number. These bounds indicate that the convergence rate of RPCD matches that of RCD, and both are much faster than CCD on the problem defined by (11) and (17). We also note that a slightly tighter bound on the asymptotic behavior of the two-variable recurrence can be obtained from the spectral radius of the $2\times 2$ matrix governing this recurrence, in a regime in which $\delta$ is close to zero. We derive an estimate of this spectral radius in (52), using results from Appendix B. Theorem 3.9 explores the behavior of the randomized methods on the very first iteration, showing that a significant decrease can be expected just on this one iteration. (Similar results can be expected for the cyclic variant CCD, as we remark in comments following Theorem 3.9.)

Empirical verification of our analysis of RPCD, and computational comparisons with CCD and RCD, are presented in Section 4. The theoretical results are confirmed nicely in all cases. We conclude with some discussions in Section 5.

Convergence of CD Variants on Convex Quadratics

We consider the application of CCD and RPCD to the convex quadratic problem (11). This problem has solution $x^{*}=0$ with optimal objective $f^{*}=0$. We assume that the matrix $A$ is diagonally scaled so that

Under this assumption, the step of Algorithm 1 with exact line search will have the form

Some variants of CD methods applied (11) can be viewed as Gauss-Seidel methods applied to the system $Ax=0$. Cyclic CD corresponds to standard Gauss-Seidel, whereas RCD and RPCD are variants of randomized Gauss-Seidel.

Writing $A=L+D+L^{T}$, where $L$ is strictly lower triangular and $D$ is the diagonal, one epoch of the CCD method can be written as follows:

The average improvement in $f$ per epoch is obtained from the formula

To obtain a bound on this quantity, we denote the eigenvalues of $C$ by $\gamma_{i}$, $i=1,2,\dotsc,n$, and recall that the spectral radius $\rho(C)$ is $\max_{i=1,2,\dotsc,n}\,|\gamma_{i}|$. Since $A$ is positive definite, we have $\rho(C)<1$ (Golub & Van Loan, 2012, Theorem 11.2.3). We have from Gelfand’s formula (Gelfand, 1941) that

We can obtain a bound on $\rho_{\text{CCD}}(A,x^{0})$ in terms of $\rho(C)$ as follows:

We can describe each epoch of RPCD algebraically by using a permutation matrix $P_{l}$ to represent the permutation $\pi_{l}$ on epoch $l-1$. We split the matrix $P_{l}^{T}AP_{l}$ and define the operator $C_{l}$ as follows:

2 CD Variants Applied to Permutation-Invariant A𝐴A

In our search for the simplest instance of a matrix $A$ for which the superiority of randomization is observed, we arrived at the matrix (17). As mentioned above, the eigenvalues of $A$ are

The restriction in (17) ensures that $A$ has the following properties:

unit diagonals: $A_{ii}=1$, $i=1,2,\dotsc,n$;

invariant under symmetric permutations of the rows and columns, that is, $P^{T}AP=A$ for any $n\times n$ permutation matrix $P$;

$L/L_{\max}$ is close to its maximum value of $n$ when $\delta$ is small, opening a wide gap between the worst-case theoretical behaviors of CCD and RCD.

Figure 2 shows results for the CCD, RPCD, and RCD variants on the matrix $A$ from (17) with $n=100$ and two different values of $\delta$. Here, the vertical axis shows actual function values (not expected values) relative to $f(x^{0})$, for some particular $x^{0}$ whose elements are drawn i.i.d from $N(0,1)$. For both values of $\delta$, both randomized variants are much faster than CCD. For the larger value of $\delta$, RPCD has a clear advantage over RCD. Our analysis below supports these empirical observations.

We now derive expressions for the epoch iteration matrix $C$ of Section 2.1 for the specific case of the permutation-invariant matrix (17). This is needed for the analysis of RPCD on this matrix. By applying the splitting (20) to (17), we have

We have from (31b) and the properties of $E$ and $\bar{L}$ that

For the complementary case $i\geq j$, we have

3 Convergence Rates of CCD and RCD on the Permutation-Invariant A𝐴A

Here, we examine the theoretical convergence rate of CCD on the quadratic function with Hessian (17) by using the results of Sun & Ye (2016).

Recalling the rate (14) from Sun & Ye (2016, Proposition 3.1), and substituting the following quantities for (17):

(We use $\rho_{\text{CCD}}(\delta,x^{0})$ in place of $\rho_{\text{CCD}}(A,x^{0})$, to emphasize the dependence of $A$ in (17) on the parameter $\delta$.) By making the mild assumption that $\delta\leq 3/4$, this expression simplifies to

On the other hand, Sun and Ye show the following lower bound on $\rho_{\text{CCD}}(\delta,x^{0})$ (obtained by substituting from (33) into Theorem 3.1 of Sun & Ye (2016)):

By combining (34) and (35), we see that for small values of $\delta/n$, the average epoch-wise decrease in error is $\rho_{\text{CCD}}(\delta,x^{0})=1-c\delta/n^{2}$, for some moderate value of $c$. Classical numerical analysis for Gauss-Seidel derives similar dependency on $n^{2}$ for this case from the eigenvalues of $A$, $D$, and $L$; see (Samarskii & Nikolaev, 1989), Young & Rheinboldt (1971, p. 464), and Hackbusch (2016, Theorem 3.44). This dependency on $n$ is confirmed empirically, by running CCD for $A$ with the same $\delta$ but different $n$, as shown in Figure 3(a).

For RCD, we have by substituting the values in (33) into (15) that the expected per-epoch improvement in error is given by

This result suggests that complexity of RCD is $O(n^{2})$ times better than CCD for small $\delta$, and that its rate does not depend strongly on $n$. This independence of $n$ is confirmed empirically by Figure 3(b). The expression (16) suggests a slightly better complexity for RCD of roughly $|\log\epsilon|(1+\delta)/(2\delta)$ epochs, rather than $|\log\epsilon|/\delta$ epochs, corresponding to replacing $1-\delta$ in (36) by

A kind of lower bound on the per-iterate improvement of RCD on the problem (11), (17) can be found by setting $x_{i}=(-1)^{i}$, $i=1,2,\dotsc,n$ with $n$ even. It can be shown that the function values for this $x$ and the next RCD iterate $x^{+}$ are

Figure 3(c) shows that RPCD too has a convergence rate independent of $n$ on this matrix. (The performances of RPCD and RCD are quite similar on the problems graphed.) The convergence rate of CCD deteriorates with $n$, according to the predictions above.

Convergence of RPCD for the Permutation-Invariant A𝐴A

and note that $\bar{A}^{(0)}=A$ and (by comparison with (39)) that

We have the following recursive relationship between successive terms in the sequence of $\bar{A}^{(t)}$ matrices:

For any $P\in\Pi$, if $P$ shifts the $i$th position to the $j$th position, then $(PQP^{T})_{jj}=Q_{ii}$. Since the probability of taking any permutation from $\Pi$ is identical, we have that

(where ${\mathcal{P}}(\cdot)$ denotes probability). Therefore, each diagonal entry $B$ is the average over all diagonal entries of $Q$.

Consider permutations that shift the $i$th and the $j$th entries to the $k$th and the $l$th positions, respectively, that is,

Note that we always have that $i\neq j\Rightarrow k\neq l$ because permutations are bijections from and to $\{1,\dotsc,n\}$. Thus, there are $(n-2)!$ permutations in $\Pi$ with the property (44). Under the same reasoning as before, each off-diagonal entry of $B$ is the average of all off-diagonal entries of $Q$.

Finally, we obtain (43) by noting that $B_{ii}=\tau_{1}+\tau_{2}$, while $B_{ij}=\tau_{2}$ for $i\neq j$.

Note that for (46c) and (46d) we used the property $\text{trace}(AB)=\text{trace}(BA)$.

The following theorem reveals the relationship between successive matrices in the sequence $\bar{A}^{(0)},\bar{A}^{(1)},\dotsc$.

Consider solving (11) with the matrix $A$ defined in (17) using RPCD. For $\bar{A}^{(t)}$ defined in (41), with $\bar{A}^{(0)}=A$, we have

where $(\eta_{0},\nu_{0})=(\delta,1-\delta)$ and

and $d_{1},d_{2},m_{1},m_{2}$ are defined in (46).

We first prove (47) by induction. By (17), it holds at $t=0$. Now assume it holds for $t=k$, for some $\eta_{k}$ and $\nu_{k}$, then for $k+1$ we have from (41)

Because $\bar{A}^{(k)}$ is in the form (47), it is invariant to row and column permutations, that is, $P^{T}\bar{A}^{(k)}P=\bar{A}^{(k)}$ for all $P\in\Pi$. Hence,

Consider solving (11) with the matrix $A$ defined in (17) using RPCD. Then, using the notation of Theorem 3.3, we have

2 Convergence of the Two-Parameter Recurrence

It is evident from Theorems 3.3 and 3.5 (and Gelfand’s formula) that the asymptotic convergence of the expected value of $f$ is governed by $\rho(M)$, which, because of definitions (49) and (46), is a function of $\delta$ and $n$. In Figure 2(b) of Section 2.2 and Table 1 of Section 4, we see that this rate is significantly better than those obtained for RCD and CCD when $\delta$ is not too close to zero (that is, $\delta\geq.2$). In this section, we estimate the convergence rate of RPCD for $\delta$ close to zero, showing that in this regime, it is close to the rate of approximately $1-2\delta$ obtained by RCD (37), and much faster than the rate of CCD discussed in (34), (35), which is $1-c\delta/n^{2}$, for some modest value of $c$.

By appealing to Theorems 3.3 and 3.5, we obtain our main result.

Consider solving (11), (17) with $\delta\in[0,0.4]$ and $n\geq 10$ using RPCD. Then, using the notation of Theorem 3.3, we have that

Thus, we have the following bound on the convergence of the expected value of the function:

Since $(\eta_{0},\nu_{0})=(\delta,1-\delta)$, we have from (48) and using $\delta\in[0,0.4]$ that

The final claim follows directly from Theorem 3.5.

This result indicates a global linear rate of at worst $1-2\delta+4\delta^{2}$, similar to the rate (37) obtained for RCD (identical to $O(\delta)$) and much faster than the rate obtained for CCD in (34), (35).

By using slightly more refined estimates of the elements of $M$, which involve not strict upper bounds as in (A.11) but rather remainder terms containing higher powers of $\delta$ and/or $1/n$, we can obtain an estimate of $\rho(M)$. In Appendix B, we obtain the following estimates of $d_{1}$, $d_{2}$, $m_{1}$, and $m_{2}$:

By substituting these estimates into (49) and calculating the spectral radius $\rho(M)$ as the largest root of the characteristic quadratic $\det(M-\lambda I)$, we obtain

This asymptotic rate is identical to the rate for RCD (37) in the $1$, $\delta$, and $\delta^{2}$, terms, and is slightly better because of the presence of the $-2\delta/n$ term.

3 The First Iteration

We noted in the numerical experiments (Figures 3 and 4) that the decrease in $f$ over the first epoch of RPCD is rather dramatic. In fact, after just a single iteration, the function value was often of order $\delta$, for all three variants (CCD, RPCD, and RCD). The following result supports this observation.

Consider solving (11) with the matrix $A$ defined in (17) using RCD or RPCD with exact line search (19). Given any $x^{0}$, the expected function value after a single iteration satisfies

where $i$ denotes the coordinate chosen for updating at the first iteration.

Note that $i$ is chosen uniformly at random from $\{1,2,\dotsc,n\}$ for both RPCD and RCD. After one step of CD, we have

we have by taking expectation with respect to $i$ in (54) that the equality in (53) holds.

For CCD, we have from (54) with $i=1$ that

If $x^{0}$ is independent of $\delta$, we have that $f(x^{1})=O(\delta)$. However there is no guarantee that $f(x^{1})$ is substantially smaller than $f(x^{0})$. If $x^{0}$ is chosen “adversarially” in such a way that $|\mathbf{1}^{T}x^{0}|\ll\|x^{0}\|$, we may find that $f(x^{1})$ is not much smaller than $f(x^{0})$. For random choices of $x^{0}$, however, we would expect a significant decrease on the first iteration, similar to that observed for RPCD and RCD.

Computational Results

Conclusions

Recent work has shown that the problem (11) with Hessian matrix (17) is a case that reveals significant differences in performance between cyclic and randomized variants of coordinate descent. Here, we provide an analysis of the performance of random-permutations cyclic coordinate descent that sharply predicts the practical convergence behavior of this approach, showing an asymptotic convergence rate that at least matches (and is even slightly better than) that obtained by a random sampling-with-replacement scheme.

Acknowledgments

We thank two referees and the Editor-in-Chief for their comments on earlier drafts, which caused us to improve the presentation and sharpen the results of the paper.

References

Appendix A Estimating Terms in the Recurrence Matrix M𝑀M

Here we first find upper and (in some cases) lower bounds for the following quantities, for the matrix $A$ given in (17) and the corresponding value of $C$ defined in (29) and (31):

We then use these quantities to obtain bounds on $d_{1}$, $d_{2}$, $m_{1}$, and $m_{2}$ from (46). We assume throughout that $n\geq 10$ and $\delta\in[0,0.4]$.

For $\mathbf{1}^{T}C\mathbf{1}$, we have from (29) and (31) that

where $u=\bar{L}^{T}\mathbf{1}$ and $v=E^{T}\mathbf{1}$ have the following components:

(from (31a)). For $\delta\in[0,0.4]$, we have

We next seek an upper bound for $\|C^{T}\mathbf{1}\|_{2}^{2}$. We have from (32) that

We now use (32) to compute bounds on the other quantities in (A.1). We have

Noting that $[-2n+2(n-1)\delta]<0$ and $[n^{2}-2n(n-1)\delta+(n-1)^{2}\delta^{2}]>0$ for the values of $\delta$ and $n$ of interest, and using $2n\delta^{n}(1-\delta^{2})\leq(2\delta^{8})n\delta^{2}\leq.01n\delta^{2}$, we continue as follows:

Thus, dividing by $n(n-1)$, and using $\delta\in[0,0.4]$ to deduce that $(1-\delta^{2})^{-1}\leq 1+1.5\delta^{2}$, we obtain

For the corresponding lower bound, we pick up from (A.4) and again use $[-2n+2(n-1)\delta]<0$ and $[n^{2}-2n(n-1)\delta+(n-1)^{2}\delta^{2}]>0$, together with $[n^{2}-2n(n-1)\delta+(n-1)^{2}\delta^{2}]\leq n^{2}(1+\delta^{2})$ and $\delta^{2n}(1+\delta^{2})\leq.001\delta^{2}$, to obtain the following:

Note that this lower bound is strictly positive in the regime $\delta\in[0,0.4]$ and $n\geq 10$.

For $\|C\|_{F}^{2}$, we obtain from (32) that

where in (A.8) we used $2\delta^{n}\leq 2(\delta^{8})\delta^{2}\leq.01\delta^{2}$. It therefore follows that

It follows, using again $\delta\in[0,0.4]$ and $n\geq 10$, that

where we used (A.6) for the final inequality. It follows that

From the formulas (46) together with (A.2), (A.5), (A.6), (A.3), (A.9), and (A.10), and using $n\geq 10$, we have the following:

From (A.11c), (A.11d), (A.2) and (A.3), we have

For the two terms $d_{1}$ and $d_{2}$, we first need better approximations of $\|C\mathbf{1}\|_{2}^{2}$ and $\|C\|_{F}^{2}$. From (A.4), we proceed with

For $\|C\|_{F}^{2}$, we obtain from (A.7) that

Appendix C Condition (9) for g​(E​x)𝑔𝐸𝑥g(Ex) with g𝑔g Strongly Convex

Meanwhile we have by convexity of $f$ that

Dividing both sides by $(f(x)-f^{*})^{1/2}$ we obtain