Straggler Mitigation in Distributed Optimization Through Data Encoding

Introduction

Solving large-scale optimization problems has become feasible through distributed implementations. However, the efficiency can be significantly hampered by slow processing nodes, network delays or node failures. In this paper we develop an optimization framework based on encoding the dataset, which mitigates the effect of straggler nodes in the distributed computing system. Our approach can be readily adapted to the existing distributed computing infrastructure and software frameworks, since the node computations are oblivious to the data encoding.

In this paper, we focus on problems of the form

instead. In doing so, we proceed with the computation in each iteration without waiting for the stragglers, with the idea that the inserted redundancy will compensate for the lost data. The goal is to design the matrix $S$ such that, when the nodes obliviously solve the problem (2) without waiting for the slowest $(m-k)$ nodes (where $k$ is a design parameter) the achieved solution approximates the original solution $w^{*}=\operatorname*{arg\,min}_{w}f(w)$ sufficiently closely. Since in large-scale machine learning and data analysis tasks one is typically not interested in the exact optimum, but rather a “sufficiently" good solution that achieves a good generalization error, such an approximation could be acceptable in many scenarios. Note also that the use of such a technique does not preclude the use of other, non-coding straggler-mitigation strategies (e.g., and references therein), which can still be implemented on top of the redundancy embedded in the system, to potentially further improve performance.

Focusing on gradient descent and L-BFGS algorithms, we show that under a spectral condition on $S$, one can achieve an approximation of the solution of (1), by solving (2), without waiting for the stragglers. We show that with sufficient redundancy embedded, and with updates from a sufficiently large, yet strict subset of the nodes in each iteration, it is possible to deterministically achieve linear convergence to a neighborhood of the solution, as opposed to convergence in expectation (see Fig. 4). Further, one can adjust the approximation guarantee by increasing the redundancy and number of node updates waited for in each iteration. Another potential advantage of this strategy is privacy, since the nodes do not have access to raw data itself, but can still perform the optimization task over the jumbled data to achieve an approximate solution.

Although in this paper we focus on quadratic objectives and two specific algorithms, in principle our approach can be generalized to more general, potentially non-smooth objectives and constrained optimization problems, as we discuss in Section 4 ( adding a regularization term is also a simple generalization).

Our main contributions are as follows. (i) We demonstrate that gradient descent (with constant step size) and L-BFGS (with line search) applied in a coding-oblivious manner on the encoded problem, achieves (universal) sample path linear convergence to an approximate solution of the original problem, using only a fraction of the nodes at each iteration. (ii) We present three classes of coding matrices; namely, equiangular tight frames (ETF), fast transforms, and random matrices, and discuss their properties. (iii) We provide experimental results demonstrating the advantage of the approach over uncoded ($S=I$) and data replication strategies, for ridge regression using synthetic data on an AWS cluster, as well as matrix factorization for the Movielens 1-M recommendation task.

Use of data replication to aid with the straggler problem has been proposed and studied in , and references therein. Additionally, use of coding in distributed computing has been explored in . However, these works exclusively focused on using coding at the computation level, i.e., certain linear computational steps are performed in a coded manner, and explicit encoding/decoding operations are performed at each step. Specifically, used MDS-coded distributed matrix multiplication and focused on breaking up large dot products into shorter dot products, and perform redundant copies of the short dot products to provide resilience against stragglers. considers a gradient descent method on an architecture where each data sample is replicated across nodes, and designs a code such that the exact gradient can be recovered as long as fewer than a certain number of nodes fail. However, in order to recover the exact gradient under any potential set of stragglers, the required redundancy factor is on the order of the number of straggling nodes, which could mean a large amount of overhead for a large-scale system. In contrast, we show that one can converge to an approximate solution with a redundancy factor independent of network size or problem dimensions (e.g., $2$ as in Section 5).

Our technique is also closely related to randomized linear algebra and sketching techniques , used for dimensionality reduction of large convex optimization problems. The main difference between this literature and the proposed coding technique is that the former focuses on reducing the problem dimensions to lighten the computational load, whereas coding increases the dimensionality of the problem to provide robustness. As a result of the increased dimensions, coding can provide a much closer approximation to the original solution compared to sketching techniques.

Encoded Optimization Framework

Figure 1 shows a typical data-distributed computational model in large-scale optimization (left), as well as our proposed encoded model (right). Our computing network consists of $m$ machines, where machine $i$ stores $\left(\widetilde{X}_{i},\widetilde{y}_{i}\right)=\left(S_{i}X,S_{i}y\right)$ and $S=\left[S_{1}^{\top}\;\;S_{2}^{\top}\;\;\dots\;\;S_{m}^{\top}\right]^{\top}$. The optimization process is oblivious to the encoding, i.e., once the data is stored at the nodes, the optimization algorithm proceeds exactly as if the nodes contained uncoded, raw data $(X,y)$.

In each iteration $t$, the central server broadcasts the current estimate $w_{t}$, and each worker machine computes and sends to the server the gradient terms corresponding to its own partition $g_{i}(w_{t}):=\widetilde{X}_{i}^{\top}(\widetilde{X}_{i}w_{t}-\widetilde{y}_{i})$.

Note that this framework of distributed optimization is typically communication-bound, where communication over a few slow links constitute a significant portion of the overall computation time. We consider a strategy where at each iteration $t$, the server only uses the gradient updates from the first $k$ nodes to respond in that iteration, thereby preventing such slow links and straggler nodes from stalling the overall computation:

where $A_{t}\subseteq[m]$, $|A_{t}|=k$ are the indices of the first $k$ nodes to respond at iteration $t$, $\eta:=\frac{k}{m}$ and $\widetilde{X}_{A}=\left[S_{i}X\right]_{i\in A_{t}}$. (Similarly, $S_{A}=\left[S_{i}\right]_{i\in A_{t}}$.) Given the gradient approximation, the central server then computes a descent direction $d_{t}$ through the history of gradients and parameter estimates. For the remaining nodes $i\not\in A_{t}$, the server can either send an interrupt signal, or simply drop their updates upon arrival, depending on the implementation.

Next, the central server chooses a step size $\alpha_{t}$, which can be chosen as constant, decaying, or through exact line search Note that exact line search is not more expensive than backtracking line search for a quadratic loss, since it only requires a single matrix-vector multiplication. by having the workers compute $\widetilde{X}d_{t}$ that is needed to compute the step size. We again assume the central server only hears from the fastest $k$ nodes, denoted by $D_{t}\subseteq[m]$, where $D_{t}\neq A_{t}$ in general, to compute

where $\widetilde{X}_{D}=\left[S_{i}X\right]_{i\in D_{t}}$, and $0<\nu<1$ is a back-off factor of choice.

Our goal is to especially focus on the case $k<m$, and design an encoding matrix $S$ such that, for any sequence of sets $\left\{A_{t}\right\}$, $\left\{D_{t}\right\}$, $f(w_{t})$ universally converges to a neighborhood of $f(w^{*})$. Note that in general, this scheme with $k<m$ is not guaranteed to converge for traditionally batch methods like L-BFGS. Additionally, although the algorithm only works with the encoded function $\widetilde{f}$, our goal is to provide a convergence guarantee in terms of the original function $f$.

Algorithms and Convergence Analysis

Let the smallest and largest eigenvalues of $X^{\top}X$ be denoted by $\mu>0$ and $M>0$, respectively.

Let $\eta$ with $\frac{1}{\beta}<\eta\leq 1$ be given. In order to prove convergence,we will consider a family of matrices $\left\{S^{(\beta)}\right\}$ where $\beta$ is the aspect ratio (redundancy factor), such that for any $\epsilon>0$, and any $A\subseteq[m]$ with $\left|A\right|=\eta m$,

for sufficiently large $\beta\geq 1$, where $S_{A}=\left[S_{i}\right]_{i\in A}$ is the submatrix associated with subset $A$ (we drop dependence on $\beta$ for brevity). Note that this is similar to the restricted isometry property (RIP) used in compressed sensing , except that (4) is only required for submatrices of the form $S_{A}$. Although this condition is needed to prove worst-case convergence results, in practice the proposed encoding scheme can work well even when it is not exactly satisfied, as long as the bulk of the eigenvalues of $S_{A}^{\top}S_{A}$ lie within a small interval $\left[1-\epsilon,1+\epsilon\right]$. We will discuss several specific constructions and their relation to property (4) in Section 4.

We consider gradient descent with constant step size, i.e.,

The following theorem characterizes the convergence of the encoded problem under this algorithm.

Let $f_{t}=f(w_{t})$, where $w_{t}$ is computed using gradient descent with updates from a set of (fastest) workers $A_{t}$, with constant step size $\alpha_{t}\equiv\alpha=\frac{2\zeta}{M(1+\epsilon)}$ for some $0<\zeta\leq 1$, for all $t$. If $S$ satisfies (4) with $\epsilon>0$, then for all sequences of $\{A_{t}\}$ with cardinality $\left|A_{t}\right|=k$,

where $\kappa=\frac{1+\epsilon}{1-\epsilon}$, and $\gamma_{1}=\left(1-\frac{4\mu\zeta(1-\zeta)}{M\left(1+\epsilon\right)}\right)$, and $f_{0}=f(w_{0})$ is the initial objective value.

The proof is provided in Appendix B, which relies on the fact that the solution to the effective “instantaneous" problem corresponding to the subset $A_{t}$ lies in the set $\{w:f(w)\leq\kappa^{2}f(w^{*})\}$, and therefore each gradient descent step attracts the estimate towards a point in this set, which must eventually converge to this set. Note that in order to guarantee linear convergence, we need $\kappa\gamma_{1}<1$, which can be ensured by property (4).

Theorem 1 shows that gradient descent over the encoded problem, based on updates from only $k<m$ nodes, results in deterministically linear convergence to a neighborhood of the true solution $w^{*}$, for sufficiently large $k$, as opposed to convergence in expectation. Note that by property (4), by controlling the redundancy factor $\beta$ and the number of nodes $k$ waited for in each iteration, one can control the approximation guarantee. For $k=m$ and $S$ designed properly (see Section 4), then $\kappa=1$ and the optimum value of the original function $f\left(w^{*}\right)$ is reached.

Although L-BFGS is originally a batch method, requiring updates from all nodes, its stochastic variants have also been proposed recently . The key modification to ensure convergence is that the Hessian estimate must be computed via gradient components that are common in two consecutive iterations, i.e., from the nodes in $A_{t}\cap A_{t-1}$. We adapt this technique to our scenario. For $t>0$, define $u_{t}:=w_{t}-w_{t-1}$, and

Then once the gradient terms $\left\{g_{t}\right\}_{i\in A_{t}}$ are collected, the descent direction is computed by $d_{t}=-B_{t}\widetilde{g}_{t}$, where $B_{t}$ is the inverse Hessian estimate for iteration $t$, which is computed by

For our convergence result for L-BFGS, we need another assumption on the matrix $S$, in addition to (4). Defining $\breve{S}_{t}=\left[S_{i}\right]_{i\in A_{t}\cap A_{t-1}}$ for $t>0$, we assume that for some $\delta>0$,

for all $t>0$. Note that this requires that one should wait for sufficiently many nodes to finish so that the overlap set $A_{t}\cap A_{t-1}$ has more than a fraction $\frac{1}{\beta}$ of all nodes, and thus the matrix $\breve{S}_{t}$ can be full rank. This is satisfied if $\eta\geq\frac{1}{2}+\frac{1}{2\beta}$ in the worst-case, and under the assumption that node delays are i.i.d., it is satisfied in expectation if $\eta\geq\frac{1}{\sqrt{\beta}}$. However, this condition is only required for a worst-case analysis, and the algorithm may perform well in practice even when this condition is not satisfied. The following lemma shows the stability of the Hessian estimate.

If (5) is satisfied, then there exist constants $c_{1},c_{2}>0$ such that for all $t$, the inverse Hessian estimate $B_{t}$ satisfies $c_{1}I\preceq B_{t}\preceq c_{2}I$.

The proof, provided in Appendix A, is based on the well-known trace-determinant method. Using Lemma 1, we can show the following result.

Let $f_{t}=f(w_{t})$, where $w_{t}$ is computed using L-BFGS as described above, with gradient updates from machines $A_{t}$, and line search updates from machines $D_{t}$. If $S$ satisfies (4) and (5), for all sequences of $\{A_{t}\},\{D_{t}\}$ with $\left|A_{t}\right|=\left|D_{t}\right|=k$,

where $\kappa=\frac{1+\epsilon}{1-\epsilon}$, and $\gamma_{2}=\left(1-\frac{4\mu c_{1}c_{2}}{M\left(c_{1}+c_{2}\right)^{2}}\right)$, and $f_{0}=f(w_{0})$ is the initial objective value.

The proof is provided in Appendix B. Similar to Theorem 1, the proof is based on the observation that the solution of the effective problem at time $t$ lies in a bounded set around the true solution $w^{*}$. As in gradient descent, coding enables linear convergence deterministically, unlike the stochastic and multi-batch variants of L-BFGS .

Although we focus on quadratic cost functions and two specific algorithms, our approach can potentially be generalized for objectives of the form $\left\|Xw-y\right\|^{2}+h(w)$ for a simple convex function $h$, e.g., LASSO; or constrained optimization $\min_{w\in\mathcal{C}}\left\|Xw-y\right\|^{2}$ (see ); as well as other first-order algorithms used for such problems, e.g., FISTA . In the next section we demonstrate that the codes we consider have desirable properties that readily extend to such scenarios.

Code Design

We consider three classes of coding matrices: tight frames, fast transforms, and random matrices.

Therefore, the solution to the encoded problem satisfies the optimality condition for the original problem as well:

and if $f$ is also strongly convex, then $\widetilde{w}^{*}=w^{*}$ is the unique solution. Note that since the computation is coding-oblivious, this is not true in general for an arbitrary full rank matrix, and this is, in addition to property (4), a desired property of the encoding matrix. In fact, this equivalency extends beyond smooth unconstrained optimization, in that

for any convex constraint set $\mathcal{C}$, as well as

for any non-smooth convex objective term $h(x)$, where $\partial h$ is the subdifferential of $h$. This means that tight frames can be promising encoding matrix candidates for non-smooth and constrained optimization too. In , it was shown that when $\left\{A_{t}\right\}$ is static, equiangular tight frames allow for a close approximation of the solution for constrained problems.

A tight frame is equiangular if $\left|\langle\phi_{i},\phi_{j}\rangle\right|$ is constant across all pairs $(i,j)$ with $i\neq j$.

Let $F=\left\{\phi_{i}\right\}_{i=1}^{n\beta}$ be a tight frame. Then $\omega(F)\geq\sqrt{\frac{\beta-1}{2n\beta-1}}$. Moreover, equality is satisfied if and only if $F$ is an equiangular tight frame.

Therefore, an ETF minimizes the correlation between its individual elements, making each submatrix $S_{A}^{\top}S_{A}$ as close to orthogonal as possible, which is promising in light of property (4). We specifically evaluate Paley and Hadamard ETFs (not to be confused with Hadamard matrix, which is discussed next) in our experiments. We also discuss Steiner ETFs in Appendix D, which enable efficient implementation.

Another computationally efficient method for encoding is to use fast transforms: Fast Fourier Transform (FFT), if $S$ is chosen as a subsampled DFT matrix, and the Fast Walsh-Hadamard Transform (FWHT), if $S$ is chosen as a subsampled real Hadamard matrix. In particular, one can insert rows of zeroes at random locations into the data pair $(X,y)$, and then take the FFT or FWHT of each column of the augmented matrix. This is equivalent to a randomized Fourier or Hadamard ensemble, which is known to satisfy the RIP with high probability .

A natural choice of encoding is using i.i.d. random matrices. Although such random matrices do not have the computational advantages of fast transforms or the optimality-preservation property of tight frames, their eigenvalue behavior can be characterized analytically. In particular, using the existing results on the eigenvalue scaling of large i.i.d. Gaussian matrices and union bound, it can be shown that

as $n\to\infty$, where $\sigma_{i}$ denotes the $i$th singular value. Hence, for sufficiently large redundancy and problem dimension, i.i.d. random matrices are good candidates for encoding as well. However, for finite $\beta$, even if $k=m$, in general for this encoding scheme the optimum of the original problem is not recovered exactly.

Using the analytical bounds (6)–(7) on i.i.d. Gaussian matrices, one can see that such matrices satisfy (4) with $\epsilon=O\left(\frac{1}{\sqrt{\beta\eta}}\right)$, independent of problem dimensions or number of nodes $m$. Although we do not have tight eigenvalue bounds for subsampled ETFs, numerical evidence (Figure 3) suggests that they may satisfy (4) with smaller $\epsilon$ than random matrices, and thus we believe that the required redundancy in practice is even smaller for ETFs.

Note that our theoretical results focus on the extreme eigenvalues due to a worst-case analysis; in practice, most of the energy of the gradient will be on the eigen-space associated with the bulk of the eigenvalues, which the following proposition suggests can be mostly 1 (also see Figure 3), which means even if (4) is not satisfied, the gradient (and the solution) can be approximated closely for a modest redundancy, such as $\beta=2$. The following result is a consequence of the Cauchy interlacing theorem, and the definition of tight frames.

If the rows of $S$ are chosen to form an ETF with redundancy $\beta$, then for $\eta\geq 1-\frac{1}{\beta}$, $\frac{1}{\beta}S_{A}^{\top}S_{A}$ has $n(1-\beta\eta)$ eigenvalues equal to 1.

Numerical Results

We generate the elements of matrix $X$ i.i.d. $\sim N(0,1)$, the elements of $y$ i.i.d. $\sim N(0,p)$, for dimensions $(n,p)=(4096,6000)$, and solve the problem $\min_{w}\frac{1}{2\beta n}\left\|\widetilde{X}w-\widetilde{y}\right\|^{2}+\frac{\lambda}{2}\|w\|^{2}$, for regularization parameter $\lambda=0.05$. We evaluate column-subsampled Hadamard matrix with redundancy $\beta=2$ (encoded using FWHT for fast encoding), data replication with $\beta=2$, and uncoded schemes. We implement distributed L-BFGS as described in Section 3 on an Amazon EC2 cluster using the mpi4py Python package, over $m=32$ m1.small worker node instances, and a single c3.8xlarge central server instance. We assume the central server encodes and sends the data variables to the worker nodes (see Appendix D for a discussion of how to implement this more efficiently).

Figure 4 shows the result of our experiments, which are aggregated over 20 trials. As baselines, we consider the uncoded scheme, as well as a replication scheme, where each uncoded partition is replicated $\beta=2$ times across nodes, and the server uses the faster copy in each iteration. It can be seen from the right figure that one can speed up computation by reducing $\eta$ from 1 to, for instance, 0.375, resulting in more than $40\%$ reduction in the runtime. Note that in this case, uncoded L-BFGS fails to converge, whereas the Hadamard-coded case stably converges. We also observe that the data replication scheme converges on average, but in the worst case, the convergence is much less smooth, since the performance may deteriorate if both copies of a partition are delayed.

We choose $\mu=3$, $p=15$, and $\lambda=10$, which achieves a test RMSE 0.861, close to the current best test RMSE on this dataset using matrix factorizationhttp://www.mymedialite.net/examples/datasets.html.

Problem (8) is often solved using alternating minimization, minimizing first over all $\left(x_{i},u_{i}\right)$, and then all $\left(y_{j},v_{j}\right)$, in repetition. Each such step further decomposes by row and column, made smaller by the sparsity of $R$. To solve for $\left(x_{i},u_{i}\right)$, we first extract $I_{i}=\{j\mid r_{ij}\text{ is observed}\}$, and solve the resulting sequence of regularized least squares problems in the variables $w_{i}=[x_{i}^{\top},u_{i}]^{\top}$ distributedly using coded L-BFGS; and repeat for $w=[y_{j}^{\top},v_{j}]^{\top}$, for all $j$. As in the first experiment, distributed coded L-BFGS is solved by having the master node encoding the data locally, and distributing the encoded data to the worker nodes (Appendix D discusses how to implement this step more efficiently). The overhead associated with this initial step is included in the overall runtime in Figure 6.

The Movielens experiment is run on a single 32-core machine with 256 GB RAM. In order to simulate network latency, an artificial delay of $\Delta\sim\text{exp}(\text{10 ms})$ is imposed each time the worker completes a task. Small problem instances ($n<500$) are solved locally at the central server, using the built-in function numpy.linalg.solve. Additionally, parallelization is only done for the ridge regression instances, in order to isolate speedup gains in the L-BFGS distribution. To reduce overhead, we create a bank of encoding matrices $\left\{S_{n}\right\}$ for Paley ETF and Hadamard ETF, for $n=100,200,\ldots,3500$, and then given a problem instance, subsample the columns of the appropriate matrix $S_{n}$ to match the dimensions. Overall, we observe that encoding overhead is amortized by the speed-up of the distributed optimization.

Figure 6 gives the final performance of our distributed L-BFGS for various encoding schemes, for each of the 5 epochs, which shows that coded schemes are most robust for small $k$. A full table of results is given in Appendix C.

Acknowledgments

This work was supported in part by NSF grants 1314937 and 1423271.

References

Appendix A Lemmas

In the proofs, we will ignore the normalization constants on the objective functions for brevity. Let $f^{A}_{t}:=\|\widetilde{X}_{A}w_{t}-\widetilde{y}_{A}\|^{2}$, and $f_{A}(w):=\|\widetilde{X}_{A}w-\widetilde{y}_{A}\|^{2}$ (we set $A\equiv A_{t}$). Let $\widetilde{w}_{t}^{*}$ denote the solution to the effective “instantaneous" problem at iteration $t$, i.e., $\widetilde{w}_{t}^{*}=\operatorname*{arg\,min}_{w}\|\widetilde{X}_{A}w-\widetilde{y}_{A}\|^{2}$.

Stronger versions of the following lemma has been proved in , but we include a weakened version of this result here for completeness.

Define $e=\widetilde{w}_{t}^{*}-w^{*}$ and note that

where (a) follows by expanding and re-arranging $\left\|\widetilde{X}_{A}\widetilde{w}_{t}^{*}-\widetilde{y}_{A}\right\|^{2}\leq\left\|\widetilde{X}_{A}w^{*}-\widetilde{y}_{A}\right\|^{2}$, which is since $\widetilde{w}_{t}^{*}$ is the minimizer of this function; (b) follows by the fact that $\nabla f(w^{*})=X^{\top}(Xw^{*}-y)=0$ by optimality of $w^{*}$ for $f$; (c) follows by Cauchy-Schwarz inequality; and (d) follows by the definition of matrix norm.

Since this is true for any $c>0$, we choose $c=\frac{\lambda_{\max}+\lambda_{\min}}{2}$, which gives

Plugging this back in (9), we get $f(\widetilde{w}^{*})\leq\kappa^{2}f(w^{*})$, which completes the proof. ∎

and similarly $\widetilde{f}^{A}(w)\leq\lambda_{\max}f(w)$, we have

which can be re-arranged into the linear recursive inequality

where $\kappa=\frac{\lambda_{\max}}{\lambda_{\min}}$. By considering such inequalities for $0\leq\tau\leq t$, multiplying each by $\left(\kappa\gamma\right)^{t-\tau}$ and summing, we get

$\widetilde{f}^{A}(w)$ is $\lambda_{\min}\mu$-strongly convex.

It is sufficient to show that the minimum eigenvalue of $\widetilde{X}_{A}^{\top}\widetilde{X}_{A}$ is bounded away from zero. This can easily be shown by the fact that

where since $W^{\top}MW$ is still a positive definite matrix it has the $2\times 2$ eigen-decomposition $Q^{\top}DQ$. Defining $q_{i}=Qr_{i}$ for $i=1,2$, note that the quantity we are interested in can be equivalently represented as

where $q_{2}=Dq_{1}$. Further note that for any unit vector $v$,

and since $\|Wv\|=1$, the condition number of $D$ (the ratio of the two non-zero elements of $D$) cannot be larger than that of $M$, which is $\kappa$ (since otherwise once could find unit vectors $u_{1}=Wv_{1}$ and $u_{2}=Wv_{2}$ such that $\frac{u_{1}^{\top}Mu_{1}}{u_{2}^{\top}Mu_{2}}>\kappa$, which is a contradiction). Representing $q_{1}=\left[\cos\phi\;\;\sin\phi\right]^{\top}$ for some angle $\phi$, $\frac{q_{2}}{\|q_{2}\|}$ can then be written as $q_{2}=\left[\frac{d_{1}}{\sqrt{d_{1}^{2}+d_{2}^{2}}}\cos\phi\;\;\frac{d_{2}}{\sqrt{d_{1}^{2}+d_{2}^{2}}}\sin\phi\right]^{\top}$. Note that minimizing the inner product $\frac{q_{1}^{\top}q_{2}}{\|q_{2}\|}$ is equivalent to maximizing the function

over $\phi$. By setting the derivative to zero, we find that the maximizing $\phi$ is given by $\cos^{-1}\frac{d_{2}-d_{1}}{d_{2}+d_{1}}$. Therefore

which implies $\textbf{tr}\left(B_{t}\right)\leq(1+\epsilon)M\left(\widetilde{\sigma}+d\right)$. It can also be shown (similar to ) that

which implies $\det\left(B_{t}\right)\geq\det\left(B_{t}^{(0)}\right)\left(\frac{\epsilon\mu}{(1+\epsilon)M\left(\widetilde{\sigma}+d\right)}\right)^{\widetilde{\sigma}}$. Since $B_{t}\geq 0$, its trace is bounded above, and its determinant is bounded away from zero, there must exist $0<c_{1}\leq c_{2}$ such that

Appendix B Proofs of Theorem 1 and Theorem 2

Throughout the section, we will consider a particular iteration $t$, and denote

where $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ denote the minimum and maximum eigenvalues of a matrix.

We will also denote with $\widetilde{w}_{t}^{*}$ the solution to the effective “instantaneous" problem at iteration $t$, i.e., $\widetilde{w}_{t}^{*}=\operatorname*{arg\,min}_{w}\left\|\widetilde{X}_{A}w-\widetilde{y}_{A}\right\|^{2}$, where $A\equiv A_{t}$, and we ignore the normalization constants on the objective functions for. Finally, we define $\widetilde{f}^{A}(w):=\|\widetilde{X}_{A}w-\widetilde{y}_{A}\|^{2}$, and $\widetilde{f}^{A}_{t}:=\|\widetilde{X}_{A}w_{t}-\widetilde{y}_{A}\|^{2}$.

Using convexity, and the choices that $d_{t}=-\widetilde{g}_{t}$ and $\alpha_{t}=\alpha$, we have

where (a) follows by the fact that $S_{A}^{\top}S_{A}\preceq\lambda_{\max}I$; (b) follows since $X^{\top}X\preceq MI$; and (c) follows by strong convexity. Re-arranging this inequality, and using the definition of $\gamma_{1}$, we get

Then, Lemma 3 with $\bar{w}=\widetilde{w}_{t}$ implies

Finally, Lemma 2 implies $f(\widetilde{w}_{t}^{*})\leq\kappa^{2}f(w^{*})$, which concludes the proof.

B.2 Proof of Theorem 2

Using convexity and the closed-form expression for the step size, we have

where (a) follows by defining $z=\frac{Xd_{t}}{\|Xd_{t}\|}$; (b) follows by the fact that the term in parenthesis is an increasing function of the quadratic form $z^{\top}\breve{S}^{\top}_{t}\breve{S}_{t}z$ and by Assumption 1; (c) follows by the assumption that $X^{\top}X\preceq MI$; (d) follows by the definition of $d_{t}$; (e) follows by Lemmas 5 and 1; (f) follows by strong convexity of $\widetilde{f}$ (by Lemma 4), which implies $\|\widetilde{g}_{t}\|^{2}\geq 2\mu\left(\widetilde{f}\left(\theta_{t}\right)-\widetilde{f}\left(\widetilde{w}_{t}^{*}\right)\right)$; (g) follows by choosing $\nu=\frac{\lambda_{\min}}{\lambda_{\max}}$; and (h) follows using the definition of $\gamma_{2}$.

and hence applying first Lemma 3 with $\bar{w}=\widetilde{w}_{t}$, and then Lemma 2, we get the desired result.

Appendix C Full results of Movielens 1-M experiment

Tables 1 and 2 give the test and train RMSE for the Movielens 1-M recommendation task, with a random 80/20 train/test split.

Appendix D Efficient encoding using Steiner ETF

We first describe Steiner ETF, based on the construction proposed in .

Note that each of the $v$ rows have exactly $v-1$ non-zero elements. We construct Steiner ETF $S$ as a $v^{2}\times\frac{v(v-1)}{2}$ matrix by replacing each 1 in a row with a distinct column of $H$, and normalizing by $\sqrt{v-1}$. For instance, for the above example, we have

In general, this procedure results in a matrix $S$ with redundancy factor $\beta=2\frac{v}{v-1}$. In full generality, Steiner ETFs can be constructed for larger redundancy levels; we refer the reader to for a full discussion of these constructions.

Steiner ETF allows for a distributed and efficient implementation for encoding a given matrix $X$. Note that the encoding matrix $S$ consists of $v$ blocks, each corresponding to a row of $V$. Consider the following partition for $S$:

Multiplication can then be implemented by simply taking a Hadamard transform of the corresponding rows of $X$, whose indices are given by $B_{1,i}\cup B_{2,i}$. In the case of dimension mismatch, one can append zero rows to $X$, or remove some rows from $S$ to make the multiplication well-defined.

Therefore, one can partition the blocks $\left\{S_{i}\right\}$ across worker nodes, and worker node $k$ can read the corresponding rows of $X$ from the pool of data, given by $\bigcup_{iinI_{k}}B_{1,i}\cup B_{2,i}$, where $I_{k}$ is the blocks assigned to worker $k$, and then apply Fast Hadamard Transform for each block. Note that processing of blocks can be further parallelized within each node, using multiple cores.

In practice, we have observed that the performance of Steiner ETF significantly improves if the rows of $SX$ are shuffled after encoding. This can be implemented by having the nodes exchange rows of $SX$ after encoding; however, this incurs a significant communication cost. A more practical approach could be one where each block is assigned to multiple nodes at random, i.e., each block is encoded by multiple nodes. The worker nodes can then drop a subset of their encoded rows such that each row is retained by exactly one node, based on some pre-defined row-allocation rule. Note that this has the same effect as having the nodes randomly exchange encoded rows with each other.

One might raise the point that many large-scale datasets are sparse, and this sparsity is lost after encoding, significantly increasing the memory usage. To address this issue, first consider the case where each node has access to the entire dataset $(X,y)$, for the sake of argument. Then a worker node would compute the gradient corresponding to block $i$ using the order of operations represented by the following parenthesization:

since this would only require matrix-vector multiplications. Now, note that for a given block $S_{i}$, there are only $v-1$ non-zero rows, out of $v(v-1)/2$. Therefore, in order to compute $g_{i}(w_{t})$ in the above order of operations, one would only need to store the rows of $X$ and $y$ corresponding to these non-zero rows, and then apply the transformation corresponding to $S_{i}$ whenever needed. If each node is assigned $\frac{v}{m}$ blocks, this would require storing at most $\frac{v(v-1)}{m}$ sparse rows, which means that the memory usage would only increase by a constant factor (on the order of redundancy factor).