Breaking Locality Accelerates Block Gauss-Seidel

Introduction

The randomized Gauss-Seidel method is a commonly used iterative algorithm to compute the solution of an $n\times n$ linear system $Ax=b$ by updating a single coordinate at a time in a randomized order. While this approach is known to converge linearly to the true solution when $A$ is positive definite (see e.g. ), in practice it is often more efficient to update a small block of coordinates at a time due to the effects of cache locality.

In extending randomized Gauss-Seidel to the block setting, a natural question that arises is how one should sample the next block. At one extreme a fixed partition of the coordinates is chosen ahead of time. The algorithm is restricted to randomly selecting blocks from this fixed partitioning, thus favoring data locality. At the other extreme we break locality by sampling a new set of random coordinates to form a block at every iteration.

Our first contribution in this paper is to show that, when compared to the random coordinate selection model, the fixed partition model can perform very poorly in terms of iteration complexity to reach a pre-specified error. Specifically, we present a family of instances (similar to the matrices recently studied by Lee and Wright ) where non-accelerated Gauss-Seidel with random coordinate selection performs arbitrarily faster than both non-accelerated and even accelerated Gauss-Seidel, using any fixed partition. Our result thus shows the importance of the sampling strategy and that acceleration cannot make up for a poor choice of sampling distribution.

In the process of deriving our results, we also develop a general proof framework for randomized accelerated methods based on Wilson et al. which avoids the use of estimate sequences in favor of an explicit Lyapunov function. Using our proof framework we are able to recover recent results on accelerated coordinate descent. Furthermore, our proof framework allows us to immediately transfer our results on Gauss-Seidel over to the randomized accelerated Kaczmarz algorithm, extending a recent result by Liu and Wright on updating a single constraint at a time to the block case.

Finally, we empirically demonstrate that despite its theoretical nuances, accelerated Gauss-Seidel using random coordinate selection can provide significant speedups in practical applications over Gauss-Seidel with fixed partition sampling, as well as the classical conjugate-gradient (CG) algorithm. As an example, for a kernel ridge regression (KRR) task in machine learning on the augmented CIFAR-10 dataset ($n=250,000$), acceleration with random coordinate sampling performs up to $1.5\times$ faster than acceleration with a fixed partitioning to reach an error tolerance of $10^{-2}$, with the gap substantially widening for smaller error tolerances. Furthermore, it performs over $3.5\times$ faster than conjugate-gradient on the same task.

Background

We assume that we are given an $n\times n$ matrix $A$ which is positive definite, and an $n$ dimensional response vector $b$. We also fix an integer $p$ which denotes a block size. Under the assumption of $A$ being positive definite, the function $f(x)=\frac{1}{2}x^{\mathsf{T}}Ax-x^{\mathsf{T}}b$ is strongly convex and smooth. Recent analysis of Gauss-Seidel proceeds by noting the connection between Gauss-Seidel and (block) coordinate descent on $f$. This is the point of view we will take in this paper.

We first describe the sketching framework of and show how it yields rates on Gauss-Seidel when blocks are chosen via a fixed partition or randomly at every iteration. While we will only focus on the special case when the sketch matrix represents column sampling, the sketching framework allows us to provide a unified analysis of both cases.

Under the assumptions stated above, show that for every $k\geq 0$, the sequence (1) satisfies

2 Accelerated rates for fixed partition Gauss-Seidel

Based on the interpretation of Gauss-Seidel as block coordinate descent on the function $f$, we can use Theorem 1 of Nesterov and Stich to recover a procedure and a rate for accelerating (1) in the fixed partition case; the specific details are discussed in Section A.4.2 of the appendix. We will refer to this procedure as ACDM.

The convergence guarantee of the ACDM procedure is that for all $k\geq 0$,

Results

We now present the main results of the paper. All proofs are deferred to the appendix.

Our first result is to construct instances where Gauss-Seidel with fixed partition sampling runs arbitrarily slower than random coordinate sampling, even if acceleration is used.

Consider the family of $n\times n$ positive definite matrices $\mathscr{A}$ given by $\mathscr{A}=\{A_{\alpha,\beta}:\alpha>0,\alpha+\beta>0\}$ with $A_{\alpha,\beta}$ defined as $A_{\alpha,\beta}=\alpha I+\frac{\beta}{n}\mathbf{1}_{n}\mathbf{1}_{n}^{\mathsf{T}}$. The family $\mathscr{A}$ exhibits a crucial property that $\Pi^{\mathsf{T}}A_{\alpha,\beta}\Pi=A_{\alpha,\beta}$ for every $n\times n$ permutation matrix $\Pi$. Lee and Wright recently exploited this invariance to illustrate the behavior of cyclic versus randomized permutations in coordinate descent.

We explore the behavior of Gauss-Seidel as the matrices $A_{\alpha,\beta}$ become ill-conditioned. To do this, we consider a particular parameterization which holds the minimum eigenvalue equal to one and sends the maximum eigenvalue to infinity via the sub-family $\{A_{1,\beta}\}_{\beta>0}$. Our first proposition characterizes the behavior of Gauss-Seidel with fixed partitions on this sub-family.

Next, we perform a similar calculation under the random column sampling model.

Our next proposition states that the rate of Gauss-Seidel from (2) is tight order-wise in that for any instance there always exists a starting point which saturates the bound.

2 A Lyapunov analysis of accelerated Gauss-Seidel and Kaczmarz

Motivated by our findings, our goal is to understand the behavior of accelerated Gauss-Seidel under random coordinate sampling. In order to do this, we establish a general framework from which the behavior of accelerated Gauss-Seidel with random coordinate sampling follows immediately, along with rates for accelerated randomized Kaczmarz and the accelerated coordinate descent methods of and .

For conciseness, we describe a simpler version of our framework which is still able to capture both the Gauss-Seidel and Kaczmarz results, deferring the general version to the full version of the paper. Our general result requires a bit more notation, but follows the same line of reasoning.

Consider the following sequence $\{(x_{k},y_{k},z_{k})\}_{k\geq 0}$ defined by the recurrence

where $H_{0},H_{1},...$ are independent realizations of $H$ and $\tau$ is a parameter to be chosen. Following , we construct a candidate Lyapunov function $V_{k}$ for the sequence (10) defined as

The following theorem demonstrates that $V_{k}$ is indeed a Lyapunov function for $(x_{k},y_{k},z_{k})$.

holds for a.e. $H$, where $\Phi(x;H)=x-H\nabla f(x)$. Set $\tau$ in (10) as $\tau=\sqrt{\mu/\nu}$, with

We now proceed to specialize Theorem 3.4 to both the Gauss-Seidel and Kaczmarz settings.

Note that in the setting of Theorem 3.5, by the definition of $\nu$ and $\mu$, it is always the case that $\nu\leq 1/\mu$. Therefore, the iteration complexity of acceleration is at least as good as the iteration complexity without acceleration.

We conclude our discussion of Gauss-Seidel by describing the analogue of Proposition 3.3 for Algorithm 1, which shows that our analysis in Theorem 3.5 is tight order-wise. The following proposition applies to ACDM as well; we show in the full version of the paper how ACDM can be viewed as a special case of Algorithm 1.

2.2 Accelerated Kaczmarz

The argument for Theorem 3.5 can be slightly modified to yield a result for randomized accelerated Kaczmarz in the sketching framework, for the case of a consistent overdetermined linear system.

Specifically, suppose we are given an $m\times n$ matrix $A$ which has full column rank, and $b\in\mathcal{R}(A)$. Our goal is to recover the unique $x_{*}$ satisfying $Ax_{*}=b$. To do this, we apply a similar line of reasoning as . We set $f(x)=\frac{1}{2}\lVert x-x_{*}\rVert^{2}_{2}$ and $H=P_{A^{\mathsf{T}}S}$, where $S$ again is our random sketching matrix. At first, it appears our choice of $f$ is problematic since we do not have access to $f$ and $\nabla f$, but a quick calculation shows that $H\nabla f(x)=(S^{\mathsf{T}}A)^{{\dagger}}S^{\mathsf{T}}(Ax-b)$. Hence, with $r_{k}=Ax_{k}-b$, the sequence (10) simplifies to

The remainder of the argument proceeds nearly identically, and leads to the following theorem.

Specialized to the setting of where each row of $A$ has unit norm and is sampled uniformly at every iteration, it can be shown (Section A.5.1) that $\nu\leq m$ and $\mu=\frac{1}{m}\lambda_{\min}(A^{\mathsf{T}}A)$. Hence, the above theorem states that the iteration complexity to reach $\varepsilon$ error is $O\left(\frac{m}{\sqrt{\lambda_{\min}(A^{\mathsf{T}}A)}}\log(1/\varepsilon)\right)$, which matches Theorem 5.1 of order-wise. However, Theorem 3.7 applies in general for any sketching matrix.

3 Specializing accelerated Gauss-Seidel to random coordinate sampling

Let $A$ be an $n\times n$ positive definite matrix and let $p$ satisfy $1<p<n$. We have that

We can now combine Theorem 3.5 with (5) to derive the following upper bound on the iteration complexity of accelerated Gauss-Seidel with random coordinates as

We conclude our results by illustrating our bounds on a simple example. Consider the sub-family $\{A_{\delta}\}_{\delta>0}\subseteq\mathscr{A}$, with

Related Work

We split the related work into two broad categories of interest: (a) work related to coordinate descent (CD) methods on convex functions and (b) randomized solvers designed for solving consistent linear systems.

When $A$ is positive definite, Gauss-Seidel can be interpreted as an instance of coordinate descent on a strongly convex quadratic function. We therefore review related work on both non-accelerated and accelerated coordinate descent, focusing on the randomized setting instead of the more classical cyclic order or Gauss-Southwell rule for selecting the next coordinate. See for a discussion on non-random selection rules, for a comparison of random selection versus Gauss-Southwell, and for efficient implementations of Gauss-Southwell.

Nesterov’s original paper in first considered randomized CD on convex functions, assuming a partitioning of coordinates fixed ahead of time. The analysis included both non-accelerated and accelerated variants for convex functions. This work sparked a resurgence of interest in CD methods for large problems. Most relevant to our paper are extensions to the block setting , handling arbitrary sampling distributions , and second order updates for quadratic functions .

For accelerated CD, Lee and Sidford generalize the analysis of Nesterov . While the analysis of was limited to selecting a single coordinate at a time, several follow on works generalize to block and non-smooth settings. More recently, both Allen-Zhu et al. and Nesterov and Stich independently improve the results of by using a different non-uniform sampling distribution. One of the most notable aspects of the analysis in is a departure from the (probabilistic) estimate sequence framework of Nesterov. Instead, the authors construct a valid Lyapunov function for coordinate descent, although they do not explicitly mention this. In our work, we make this Lyapunov point of view explicit. The constants in our acceleration updates arise from a particular discretization and Lyapunov function outlined from Wilson et al. . Using this framework makes our proof particularly transparent, and allows us to recover results for strongly convex functions from and as a special case.

From the numerical analysis side both the Gauss-Seidel and Kaczmarz algorithm are classical methods. Strohmer and Vershynin were the first to prove a linear rate of convergence for randomized Kaczmarz, and Leventhal and Lewis provide a similar kind of analysis for randomized Gauss-Seidel. Both of these were in the single constraint/coordinate setting. The block setting was later analyzed by Needell and Tropp . More recently, Gower and Richtárik provide a unified analysis for both randomized block Gauss-Seidel and Kaczmarz in the sketching framework. We adopt this framework in this paper. Finally, Liu and Wright provide an accelerated analysis of randomized Kaczmarz once again in the single constraint setting and we extend this to the block setting.

Experiments

In this section we experimentally validate our theoretical results on how our accelerated algorithms can improve convergence rates. Our experiments use a combination of synthetic matrices and matrices from large scale machine learning tasks.

Setup. We run all our experiments on a 4 socket Intel Xeon CPU E7-8870 machine with 18 cores per socket and 1TB of DRAM. We implement all our algorithms in Python using numpy, and use the Intel MKL library with 72 OpenMP threads for numerical operations. We report errors as relative errors, i.e. $\lVert x_{k}-x_{*}\rVert_{A}^{2}/\lVert x_{*}\rVert_{A}^{2}$. Finally, we use the best values of $\mu$ and $\nu$ found by tuning each experiment.

We implement fixed partitioning by creating random blocks of coordinates at the beginning of the experiment and cache the corresponding matrix blocks to improve performance. For random coordinate sampling, we select a new block of coordinates at each iteration.

For our fixed partition experiments, we restrict our attention to uniform sampling. While Gower and Richtárik propose a non-uniform scheme based on $\mathbf{Tr}(S^{\mathsf{T}}AS)$, for translation-invariant kernels this reduces to uniform sampling. Furthermore, as the kernel block Lipschitz constants were also roughly the same, other non-uniform schemes also reduce to nearly uniform sampling.

Our first set of experiments numerically verify the separation between fixed partitioning sampling versus random coordinate sampling.

Figure 2 shows the progress per iteration on solving $A_{1,\beta}x=b$, with the $A_{1,\beta}$ defined in Section 3.1. Here we set $n=5000$, $p=500$, $\beta=1000$, and $b\sim N(0,I)$. Figure 2 verifies our analytical findings in Section 3.1, that the fixed partition scheme is substantially worse than uniform sampling on this instance. It also shows that in this case, acceleration provides little benefit in the case of random coordinate sampling. This is because both $\mu$ and $1/\nu$ are order-wise $p/n$, and hence the rate for accelerated and non-accelerated coordinate descent coincide. However we note that this only applies for matrices where $\mu$ is as large as it can be (i.e. $p/n$), that is instances for which Gauss-Seidel is already converging at the optimal rate (see , Lemma 4.2).

2 Kernel ridge regression

We next evaluate how fixed partitioning and random coordinate sampling affects the performance of Gauss-Seidel on large scale machine learning tasks. We use the popular image classification dataset CIFAR-10 and evaluate a kernel ridge regression (KRR) task with a Gaussian kernel. Specifically, given a labeled dataset $\{(x_{i},y_{i})\}_{i=1}^{n}$, we solve the linear system $(K+\lambda I)\alpha=Y$ with $K_{ij}=\exp(-\gamma\lVert x_{i}-x_{j}\rVert_{2}^{2})$, where $\lambda,\gamma>0$ are tunable parameters (see e.g. for background on KRR). The key property of KRR is that the kernel matrix $K$ is positive semi-definite, and hence Algorithm 1 applies.

For the CIFAR-10 dataset, we augment the datasetSimilar to https://github.com/akrizhevsky/cuda-convnet2. to include five reflections, translations per-image and then apply standard pre-processing steps used in image classification . We finally apply a Gaussian kernel on our pre-processed images and the resulting kernel matrix has $n=250000$ coordinates.

Results from running 500 iterations of random coordinate sampling and fixed partitioning algorithms are shown in Figure 4. Comparing convergence across iterations, similar to previous section, we see that un-accelerated Gauss-Seidel with random coordinate sampling is better than accelerated Gauss-Seidel with fixed partitioning. However we also see that using acceleration with random sampling can further improve the convergence rates, especially to achieve errors of $10^{-3}$ or lower.

We also compare the convergence with respect to running time in Figure 4. Fixed partitioning has better performance in practice random access is expensive in multi-core systems. However, we see that this speedup in implementation comes at a substantial cost in terms of convergence rate. For example in the case of CIFAR-10, using fixed partitions leads to an error of $1.2\times 10^{-2}$ after around 7000 seconds. In comparison we see that random coordinate sampling achieves a similar error in around 4500 seconds and is thus $1.5\times$ faster. We also note that this speedup increases for lower error tolerances.

3 Comparing Gauss-Seidel to Conjugate-Gradient

We also compared Gauss-Seidel with random coordinate sampling to the classical conjugate-gradient (CG) algorithm. CG is an important baseline to compare with, as it is the de-facto standard iterative algorithm for solving linear systems in the numerical analysis community. While we report the results of CG without preconditioning, we remark that the performance using a standard banded preconditioner was not any better. However, for KRR specifically, there have been recent efforts to develop better preconditioners, and we leave a more thorough comparison for future work. The results of our experiment are shown in Figure 4. We note that Gauss-Seidel both with and without acceleration outperform CG. As an example, we note that to reach error $10^{-1}$ on CIFAR-10, CG takes roughly 7000 seconds, compared to less than 2000 seconds for accelerated Gauss-Seidel, which is a $3.5\times$ improvement.

4 Kernel ridge regression on smaller datasets

In addition to using the large CIFAR-10 augmented dataset, we also tested our algorithms on the smaller MNISThttp://yann.lecun.com/exdb/mnist/ dataset. To generate a kernel matrix, we applied the Gaussian kernel on the raw MNIST pixels to generate a matrix $K$ with $n=60000$ rows and columns.

Results from running 500 iterations of random coordinate sampling and fixed partitioning algorithms are shown in Figure 5. We plot the convergence rates both across time and across iterations. Comparing convergence across iterations we see that random coordinate sampling is essential to achieve errors of $10^{-4}$ or lower. In terms of running time, similar to the CIFAR-10 experiment, we see that the benefits in fixed partitioning of accessing coordinates faster comes at a cost in terms of convergence rate, especially to achieve errors of $10^{-4}$ or lower.

5 Effect of block size

We next analyze the importance of the block size $p$ for the accelerated Gauss-Seidel method. As the values of $\mu$ and $\nu$ change for each setting of $p$, we use a smaller MNIST matrix for this experiment. We apply a random feature transformation to generate an $n\times d$ matrix $F$ with $d=5000$ features. We then use $A=F^{\mathsf{T}}F$ and $b=F^{\mathsf{T}}Y$ as inputs to the algorithm. Figure 2 shows the wall clock time to converge to $10^{-5}$ error as we vary the block size from $p=50$ to $p=1000$.

Increasing the block-size improves the amount of progress that is made per iteration but the time taken per iteration increases as $O(p^{3})$ (Line 5, Algorithm 1). However, using efficient BLAS-3 primitives usually affords a speedup from systems techniques like cache blocking. We see the effects of this in Figure 2 where using $p=500$ performs better than using $p=50$. We also see that these benefits reduce for much larger block sizes and thus $p=1000$ is slower.

6 Computing the μ𝜇\mu and ν𝜈\nu constants

In our last experiment, we explicitly compute the $\mu$ and $\nu$ constants from Theorem 3.5 for a few $16\times 16$ positive definite matrices constructed as follows.

Linearly spaced eigenvalues. We first draw $Q$ uniformly at random from $n\times n$ orthogonal matrices. We then construct $A_{i}=Q\Sigma_{i}Q^{\mathsf{T}}$ for $i=1,2,3$, where $\Sigma_{1}$ is diag(linspace(1, 10, 16)), $\Sigma_{2}$ is diag(linspace(1, 100, 16)), and $\Sigma_{3}$ is diag(linspace(1, 1000, 16)).

Tridiagonal matrix. We let $A$ be a tridiagonal matrix with the diagonal value equal to one, and the off diagonal value equal to $(\delta-a)/(2\cos(\pi n/(n+1)))$ for $\delta=1/10$. The matrix has a minimum eigenvalue of $\delta$.

Figure 6 shows the results of our computation for the linearly spaced eigenvalues ensemble, the random Wishart ensemble and the other deterministic structured matrices. Alongside with the actual $\nu$ values, we plot the bound given for each instance by Lemma 3.8. From the figures we see that our bound is quite close to the computed value of $\nu$ for circulant matrices and for random matrices with linearly spaced eigenvalues with small $\kappa$. We plan to extend our analysis to derive a tighter bound in the future.

Conclusion

In this paper, we extended the accelerated block Gauss-Seidel algorithm beyond fixed partition sampling. Our analysis introduced a new data-dependent parameter $\nu$ which governs the speed-up of acceleration. Specializing our theory to random coordinate sampling, we derived an upper bound on $\nu$ which shows that well conditioned blocks are a sufficient condition to ensure speedup. Experimentally, we showed that random coordinate sampling is readily accelerated beyond what our bound suggests.

The most obvious question remains to derive a sharper bound on the $\nu$ constant from Theorem 3.5. Another interesting question is whether or not the iteration complexity of random coordinate sampling is always bounded above by the iteration complexity with fixed coordinate sampling.

We also plan to study an implementation of accelerated Gauss-Seidel in a distributed setting . The main challenges here are in determining how to sample coordinates without significant communication overheads, and to efficiently estimate $\mu$ and $\nu$. To do this, we wish to explore other sampling schemes such as shuffling the coordinates at the end of every epoch .

Acknowledgements

We thank Ross Boczar for assisting us with Mathematica support for non-commutative algebras, Orianna DeMasi for providing useful feedback on earlier drafts of this manuscript, and the anonymous reviewers for their helpful feedback. ACW is supported by an NSF Graduate Research Fellowship. BR is generously supported by ONR awards N00014-11-1-0723 and N00014-13-1-0129, NSF award CCF-1359814, the DARPA Fundamental Limits of Learning (Fun LoL) Program, a Sloan Research Fellowship, and a Google Research Award. This research is supported in part by DHS Award HSHQDC-16-3-00083, NSF CISE Expeditions Award CCF-1139158, DOE Award SN10040 DE-SC0012463, and DARPA XData Award FA8750-12-2-0331, and gifts from Amazon Web Services, Google, IBM, SAP, The Thomas and Stacey Siebel Foundation, Apple Inc., Arimo, Blue Goji, Bosch, Cisco, Cray, Cloudera, Ericsson, Facebook, Fujitsu, HP, Huawei, Intel, Microsoft, Mitre, Pivotal, Samsung, Schlumberger, Splunk, State Farm and VMware.

References

Appendix A.1 Preliminaries

In what follows, unless stated otherwise, whenever we discuss a partition of $[n]$ we assume that the partition is given by $\bigcup_{i=1}^{n/p}J_{i}$, where

This is without loss of generality because for any arbitrary equal sized partition of $[n]$, there exists a permutation matrix $\Pi$ such that all our results apply by the change of variables $A\leftarrow\Pi^{\mathsf{T}}A\Pi$ and $b\leftarrow\Pi^{\mathsf{T}}b$.

Appendix A.2 Proofs for Separation Results (Section 3.1)

Recall the family of $n\times n$ positive definite matrices $\mathscr{A}$ defined in (17) as

We first gather some elementary formulas. By the matrix inversion lemma,

The fact that the right hand side is independent of $S$ is the key property which makes our calculations possible. Indeed, we have that

With these formulas in hand, our next proposition gathers calculations for the case when $S$ represents uniformly choosing $p$ columns without replacement.

First, we have the following elementary expectation calculations,

To compute $G_{\alpha,\beta}$, we simply plug (24) and (25) into (21). After simplification,

From this formula for $G_{\alpha,\beta}$, (22) follows immediately.

Next, we note for any $r,q$, using the properties that $S^{\mathsf{T}}S=I$, $\mathbf{1}_{n}^{\mathsf{T}}S\mathbf{1}_{p}=p$, and $\mathbf{1}_{p}^{\mathsf{T}}\mathbf{1}_{p}=p$, we have that

Taking expectations of both sides of the above equation and using the formulas in (24), (25), (26), and (27),

We now set $r=\alpha^{-1}$, $q=-\frac{\beta/n}{\alpha(\alpha+\beta p/n)}$, and $\gamma,\eta$ from (28) to reach the desired formula for (23). ∎

Consider the family of $n\times n$ positive definite matrices $\{A_{\alpha,\beta}\}$ from (18), and let $n$, $p$, and $S$ be described as in the preceding paragraph. We have that

Once again, the expectation calculations are

A.2.2 Proof of Proposition 3.3

Unrolling this recursion yields for all $k\geq 0$,

Choose $A^{1/2}e_{0}=v$, where $v$ is an eigenvector of $I-A^{1/2}GA^{1/2}$ with eigenvalue $\lambda_{\max}(I-A^{1/2}GA^{1/2})=1-\lambda_{\min}(GA)=1-\mu$. Now by Jensen’s inequality,

Appendix A.3 Proofs for Convergence Results (Section 3.2)

It is easily verified that $(x,y,z)=(x_{*},x_{*},x_{*})$ is a fixed point of the aforementioned dynamical system. Our goal for now is to describe conditions on $f$, $\mu$, and $\tau$ such that the sequence of updates (31a), (31b), and (31c) converges to this fixed point. As described in Wilson et al. , our main strategy for proving convergence will be to introduce the following Lyapunov function

Furthermore, suppose that $\nu>0$ satisfies

Then as long as we set $\tau>0$ such that $\tau$ satisfies for almost every $\omega\in\Omega$,

we have that $V_{k}$ defined in (32) satisfies for all $k\geq 0$,

First, recall the following two point equality valid for any vectors $a,b,c\in V$ in a real inner product space $V$,

Above, (40a) follows from $\mu$-strong convexity, (40b) and (40e) both use the definition of the sequence (31), (40c) follows from $L$-Lipschitz gradients, (40d) uses the two-point inequality (37), and the last inequality follows from the assumption of $\tau\leq\sqrt{\frac{\mu}{L}}$. The claim (36) now follows by re-arrangement. ∎

Next, we describe how to recover Theorem 3.5 from Theorem A.3.1. We do this by applying Theorem A.3.1 to the function $f(x)=\frac{1}{2}x^{\mathsf{T}}Ax-x^{\mathsf{T}}b$.

It remains to check the gradient inequality (33) and compute the strong convexity and Lipschitz parameters. These computations fall directly from the calculations made in Theorem 1 of , but we replicate them here for completeness.

To check the gradient inequality (33), because $f$ is a quadratic function, its second order Taylor expansion is exact. Hence for almost every $\omega\in\Omega$,

Hence the inequality (33) holds with equality.

The strong convexity bound now follows since

Hence, we can upper bound $V_{0}$ as follows

On the other hand, we have that $\frac{1}{2}\lVert y_{k}-x_{*}\rVert^{2}_{A}\leq V_{k}$. Putting the inequalities together,

where the first inequality holds by Jensen’s inequality. The claimed inequality (14) now follows.

A.3.2 Proof of Proposition 3.6

We first state and prove an elementary linear algebra fact which we will use below in our calculations.

Let $A,B,C,D$ be $n\times n$ diagonal matrices, and define $M=\begin{bmatrix}A&B\\ C&D\end{bmatrix}$. The eigenvalues of $M$ are given by the union of the eigenvalues of the $2\times 2$ matrices

where $A_{i},B_{i},C_{i},D_{i}$ denote the $i$-th diagonal entry of $A,B,C,D$ respectively.

Now we proceed with the proof of Proposition 3.6. Define $e_{k}=\begin{bmatrix}y_{k}-x_{*}\\ z_{k}-x_{*}\end{bmatrix}$. It is easy to see from the definition of Algorithm 1 that $\{e_{k}\}$ satisfies the recurrence

Define $P=\begin{bmatrix}A&0\\ 0&\mu G^{-1}\end{bmatrix}$. By taking and iterating expectations,

Denote the matrix $Q=A^{1/2}G^{1/2}$. Unrolling the recurrence above yields that

Write the SVD of $Q$ as $Q=U\Sigma V^{\mathsf{T}}$. Both $U$ and $V$ are $n\times n$ orthonormal matrices. It is easy to see that $R^{k}$ is given by

Suppose we choose $P^{1/2}e_{0}$ to be a right singular vector of $R^{k}$ corresponding to the maximum singular value $\sigma_{\max}(R^{k})$. Then we have that

where $\rho(\cdot)$ denotes the spectral radius. The first inequality is Jensen’s inequality, and the second inequality uses the fact that the spectral radius is bounded above by any matrix norm. The eigenvalues of $R^{k}$ are the $k$-th power of the eigenvalues of $R$ which, using the similarity transform (42) along with Proposition A.3.2, are given by the eigenvalues of the $2\times 2$ matrices $R_{i}$ defined as

where the first inequality holds since $\mu G^{-1}\preccurlyeq A$ and the second inequality holds since $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for non-negative $a,b$.

Appendix A.4 Recovering the ACDM Result from Nesterov and Stich [15]

We next show how to recover Theorem 1 of Nesterov and Stich using Theorem A.3.1, in the case of $\alpha=1$. A nearly identical argument can also be used to recover the result of Allen-Zhu et al. under the strongly convex setting in the case of $\beta=0$. Our argument proceeds in two steps. First, we prove a convergence result for a simplified accelerated coordinate descent method which we introduce in Algorithm 2. Then, we describe how a minor tweak to ACDM shows the equivalence between ACDM and Algorithm 2.

Now consider the following accelerated randomized coordinate descent algorithm in Algorithm 2.

Theorem A.3.1 is readily applied to Algorithm 2 to give a convergence guarantee which matches the bound of Theorem 1 of Nesterov and Stich. We sketch the argument below.

We next compute the $\nu$ constant defined in (34). We do this by checking the sufficient condition that $H_{i}G^{-1}H_{i}\preccurlyeq\nu H_{i}$ for $i=1,...,m$. Doing so yields that $\nu=1$, since

To complete the argument, we set $\mu$ as the strong convexity constant and $L$ as the Lipschitz gradient constant of $f$ with respect to the $\lVert\cdot\rVert_{G^{-1}}$ norm. It is straightforward to check that

Above, (a) follows by the convexity of the maximum eigenvalue, (b) holds since $S_{i}^{\mathsf{T}}S_{i}=I$, (c) uses the fact that for any matrix $Q$ satisfying $Q^{\mathsf{T}}Q=I$ and $M$ positive semi-definite, we have $(QMQ^{\mathsf{T}})^{1/2}=QM^{1/2}Q^{\mathsf{T}}$, and (d) follows since $\lambda_{\max}(S_{i}MS_{i}^{\mathsf{T}})=\lambda_{\max}(M)$ for any $p\times p$ symmetric matrix $M$. Using the fact that $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for any non-negative $a,b$, the inequality $\sqrt{L}\leq\sum_{i=1}^{m}\sqrt{L_{i}}$ immediately follows. To conclude the proof, it remains to calculate the requirement on $\tau$ via (35). Since $\frac{\gamma_{i}}{\sqrt{\Gamma_{i}}}=\frac{p_{i}}{\sqrt{L_{i}}}=\frac{1}{\sum_{i=1}^{m}\sqrt{L_{i}}}$, we have that $\frac{\gamma_{i}}{\sqrt{\Gamma_{i}}}\leq\frac{1}{\sqrt{L}}$, and hence the requirement is that $\tau\leq\frac{\sqrt{\mu}}{\sum_{i=1}^{m}\sqrt{L_{i}}}$.

A.4.2 Relating Algorithm 2 to ACDM

For completeness, we replicate the description of the ACDM algorithm from Nesterov and Stich in Algorithm 3. We make one minor tweak in the initialization of the $A_{k},B_{k}$ sequence which greatly simplifies the exposition of what follows.

We first write the sequence produced by Algorithm 3 as

Since $\beta_{k}B_{k+1}=\mu a_{k+1}$, the $z_{k+1}$ update simplifies to

Hence as long as $\mu A_{0}=B_{0}$ (which is satisfied by our modification), we have that $\mu A_{k+1}=B_{k+1}$ for all $k\geq 0$. With this identity, we have that $\alpha_{k}=\beta_{k}$ for all $k\geq 0$. Therefore, (44) simplifies to

We now calculate the value of $\beta_{k}$. At every iteration, we have that

Combining these identities, we have shown that (43a), (43b), and (43c) simplifies to

This sequence directly coincides with the sequence generated by Algorithm 2 after a simple relabeling.

A.4.3 Accelerated Gauss-Seidel for fixed partitions from ACDM

We now describe Algorithm 4, which is the specialization of ACDM (Algorithm 3) to accelerated Gauss-Seidel in the fixed partition setting.

Appendix A.5 A Result for Randomized Block Kaczmarz

We first describe the randomized accelerated block Kaczmarz algorithm in Algorithm 5. Our main convergence result concerning Algorithm 5 is presented in Theorem A.5.1.

Hence the gradient inequality (33) holds with equality. ∎

We first state a proposition which will be useful in our analysis of $\nu$.

We will prove that for every $1\leq i\leq s$,

from which the claim immediately follows. By Schur complements, (48) holds iff

Since the eigenvalues of a Kronecker product are given by the Cartesian product of the individual eigenvalues, (48) holds. ∎

where (a) follows from Proposition A.5.2. Hence, $\nu\neq m$. On the other hand,

Appendix A.6 Proofs for Random Coordinate Sampling (Section 3.3)

Our primary goal in this section is to provide a proof of Lemma 3.8. Along the way, we prove a few other results which are of independent interest. We first provide a proof of the lower bound claim in Lemma 3.8.

Since trace commutes with expectation and respects the positive semi-definite ordering, taking trace of both sides of (49) yields that

Next, the upper bound relies on the following lemma, which generalizes Lemma 2 of .

Our proof follows the strategy in the proof of Theorem 3.2 from . First, write $P_{B}=B(B^{\mathsf{T}}B)^{{\dagger}}B^{\mathsf{T}}$. Since $\mathcal{R}(B^{\mathsf{T}})=\mathcal{R}(B^{\mathsf{T}}B)$, we have by generalized Schur complements (see e.g. Theorem 1.20 from ) and the fact that expectation preserves the semi-definite order,

We are now in a position to prove the upper bound of Lemma 3.8. We apply Lemma A.6.2 to $M=A^{1/2}SS^{\mathsf{T}}A^{1/2}$ to conclude, using the fact that $\mathcal{R}(M)=\mathcal{R}(MM^{\mathsf{T}})$, that

Again, since conjugation by $A^{1/2}$ preserves semi-definite ordering, we have that

Using the fact that for positive definite matrices $X,Y$ we have $X\preccurlyeq Y$ iff $Y^{-1}\preccurlyeq X^{-1}$, (53) is equivalent to

Conjugating both sides by $P_{A^{1/2}S}$ and taking expectations,

Next, letting $J\subseteq 2^{[n]}$ denote the index set associated to $S$, for every $S$ we have

Plugging this calculation back into (54) yields the desired upper bound of Lemma 3.8.