# RidgeSketch: A Fast Sketching Based Solver for Large Scale Ridge Regression

Nidham Gazagnadou, Mark Ibrahim, Robert M Gower

► To cite this version:

Nidham Gazagnadou, Mark Ibrahim, Robert M Gower. RidgeSketch: A Fast Sketching Based Solver for Large Scale Ridge Regression. SIAM Journal on Matrix Analysis and Applications, 2022, 43, pp.1440 - 1468. 10.1137/21m1422963. hal-04182655

HAL Id: hal-04182655

https://telecom-paris.hal.science/hal-04182655v1

Submitted on 17 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
RIDGESKETCH: A FAST SKETCHING BASED SOLVER FOR
LARGE SCALE RIDGE REGRESSION*

NIDHAM GAZAGNADOU†, MARK IBRAHIM‡, AND ROBERT M. GOWER§

**Abstract.** We propose new variants of the sketch-and-project method for solving large scale ridge regression problems. First, we propose a new momentum alternative and provide a theorem showing it can speed up the convergence of sketch-and-project, through a fast *sublinear* convergence rate. We carefully delimit under what settings this new sublinear rate is faster than the previously known linear rate of convergence of sketch-and-project without momentum. Second, we consider combining the sketch-and-project method with new modern sketching methods such as Count sketch, SubCount sketch (a new method we propose), and subsampled Hadamard transforms. We show experimentally that when combined with the sketch-and-project method, the (Sub)Count sketch is very effective on sparse data and the standard Subsample sketch is effective on dense data. Indeed, we show that these sketching methods, combined with our new momentum scheme, result in methods that are competitive even when compared to the conjugate gradient method on real large scale data. On the contrary, we show the subsampled Hadamard transform does not perform well in this setting, despite the use of fast Hadamard transforms, and nor do recently proposed acceleration schemes work well in practice. To support all of our experimental findings, and invite the community to validate and extend our results, with this paper we are also releasing an open source software package: *RidgeSketch*. We designed this object-oriented package in Python for testing sketch-and-project methods and benchmarking ridge regression solvers. *RidgeSketch* is highly modular, and new sketching methods can easily be added as subclasses. We provide code snippets of our package in the appendix.

**Key words.** randomized numerical linear algebra, sketch-and-project methods, randomized sketching, randomized coordinate descent, momentum, acceleration

MSC codes. 65K10, 68T99, 90C06, 68W20, 65Y20, 62J07

DOI. 10.1137/21M1422963

**1. Introduction.** Consider the regression problem given by

$$ (1.1) \qquad \min_{w \in \mathbb{R}^d} \frac{1}{2} \|\mathbf{X}w - y\|_2^2 + \frac{\lambda}{2} \|w\|_2^2, $$

where $\mathbf{X} \in \mathbb{R}^{n \times d}$ is the feature matrix, $y \in \mathbb{R}^n$ the targets, and $\lambda > 0$ the regularization parameter.

The need to solve linear regression such as (1.1) occurs throughout scientific computing, where it is known as ridge regression in the statistics community [40, 48], data assimilation in weather forecasting [23], and linear least squares with Tikhonov regularization among numerical analysts [14].

Here we focus on applications where the number of rows and columns of $\mathbf{X} \in \mathbb{R}^{n \times d}$ are both large scale. With datasets now being collected automatically and electronically, there has been a surge of new stochastic incremental methods that can gracefully scale with the dimensions of the data in (1.1). Yet, it is still unclear if the new stochastic methods are capable of outperforming the classic conjugate gradient

<footnote>*Received by the editors May 26, 2021; accepted for publication (in revised form) by M. Tygert June 2, 2022; published electronically August 22, 2022.</footnote>

<footnote>https://doi.org/10.1137/21M1422963</footnote>

<footnote>†LTCI, Télécom Paris, Institut Polytechnique de Paris, Palaiseau, France (nidham.gazagnadou@gmail.com).</footnote>

<footnote>‡Facebook AI Research (FAIR), New York, NY 10003 USA (marksibrahim@fb.com).</footnote>

<footnote>§LTCI, Télécom Paris, Institut Polytechnique de Paris, Palaiseau, France, and Center for Computational Mathematics, Flatiron Institute, New York, NY 10010 USA (rgower@flatironinstitute.org).</footnote>
(CG) methods [19], even though the CG method was developed in the 1950s. We address this question by designing and implementing new stochastic methods based on randomized sketching [36, 37] together with the iterative projection methods known as the sketch-and-project methods [16]. Our main objective can be stated simply as follows:

**Main objective:** Develop variants of the sketch-and-project method that are competitive in practice, and provide new theoretical support for these variants.

We do this by exploring the use of new sketching transforms such as count sketch [7] and proposing a new momentum variant of the sketch-and-project method. To demonstrate the benefits of our new momentum and sketching variant, we show how our resulting method is competitive with the CG methods and supersedes previous momentum based variants [27]. We will also give several negative results, such as showing how, in our setting, recently developed accelerated variants [45, 15] are not viable in practice, and nor are the fast Johnson-Lindenstrauss (JL) transforms [1].

In section 2 we outline some of the background and state our main contributions. In section 3, we describe the sketch-and-project method for solving ridge regression and the different sketching methods we explored in section 4. Then, we give a convergence theory for our method without and with momentum respectively in sections 5 and 6, which we specialize in section 7 for single column sketches. To illustrate these theoretical results, we provide large scale numerical experiments highlighting our theoretical findings in section 8 and showing that our momentum sketch-and-project method competes against CG and direct solvers. Finally, in section 9 we present the accelerated version of our method and show its inefficiency in practice since it requires additional overhead costs and spectral information of **A** to compute acceleration parameters.

## 2. Background and contributions.

*Iterative sketching in linear systems.* When the dimensions $n$ and $d$ are large, direct methods for solving (1.1) can be infeasible, and iterative methods are favored. In particular, the Krylov methods including the CG algorithm [19] are the industrial standard so long as one can afford full matrix vector products and the system matrix fits in memory. On the other hand, if a single matrix vector product is considerably expensive, or the data matrix is too large to fit in memory, then iterative methods exploiting only a few rows or columns of **A** are effective. This includes, for example, randomized Kaczmarz methods [22, 43, 33], its greedy or deterministic variants [35, 12, 11, 18, 5], and the coordinate descent (CD) method [25, 28, 51]. In [16], the authors united Kaczmarz, CD, and a host of other randomized iterative methods under the sketch-and-project framework.

*Contributions.* We revisit the sketch-and-project method and examine how to make specialized variants for solving ridge regression that are competitive as compared to CGs. To do so, we are also releasing a high-quality and modular Python package called RidgeSketch<fnref n="1" /> for efficiently implementing and testing sketch-and-project methods. More information about the code and how to contribute by adding sketches or datasets is given in Appendix B.

<footnote n="1">¹Our fully documented and tested package is available at https://github.com/facebookresearch/RidgeSketch.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Areas of superiority of sublinear convergence of CD with momentum (in blue) and linear convergence of CD without momentum (in red), where the upper and lower bounds are respectively the left- and the right-hand side of (7.11).</figcaption></figure>

*Momentum*. Recently, a variant of sketch-and-project with momentum was proposed in [27]. The authors of [27] show a theoretical advantage in terms of convergence in expectation (of the first moment only), but do not present any advantage in terms of convergence in L2 or high probability. There has even more recent work [30] analyzing momentum together with sketch-and-project, and extensions to solving the linear feasibility problem [31, 32]. Yet none of these works show any theoretical benefit of using momentum over using no momentum at all. This lack of a theoretical benefit of using momentum is echoed throughout stochastic optimization: there are almost no theoretical benefits in using momentum outside of the strongly convex full batch setting [38].

*Contributions*. Here we establish the first theoretical advantage of using momentum together with the sketch-and-project method. We show that when using sketch-and-project with momentum, the last iterate enjoys a fast sublinear convergence. Without momentum, this result does not hold. Instead, without momentum, it has only been shown to hold for the average of the iterates. As depicted in Figure 1, we show that the fast sublinear rate of our new momentum variant gives a tighter complexity bound than the previously best known linear rate of convergence [17] when the desired precision is moderate and the condition number of the underlining problem is moderate to large. Our new convergence theory for momentum also suggests completely iteration dependent schedules for setting the momentum parameter. We perform extensive numerical tests showing the superiority of this new scheduling as compared to using a constant momentum.

*Acceleration*. Lee and Sidford [24] show how the Kaczmarz method could be accelerated through its connection to CD, and also compare the resulting accelerated convergence rate to the rate of the CG algorithm. However, they provide no numerical experiments, leaving it unclear if this form of acceleration can afford any practical advantage. Later, Liu and Wright [26] developed an accelerated Kaczmarz method and showed that it can be faster than CG on densely generated artificial data. But they did not provide examples of this on real data or an affordable rule for setting the acceleration parameters. More recently, it was shown that the entire family of the sketch-and-project methods could be accelerated [15]. Yet, experiments of the authors rely on a grid search that would defeat any gains in using acceleration.

*Contributions*. We investigate the possibility of developing a practical setting for the two acceleration parameters proposed in [15] that would result in a robust
performance gain over standard sketch-and-project. Unlike Nesterov's acceleration for gradient descent in the convex setting, there are no default parameter settings that work consistently across a significant class of problems. We show through a careful grid search that finding good parameters is like "looking for a needle in the haystack" and that identifying any practical settings for these parameters is virtually impossible. We thus recommend that, to advance the use of acceleration for linear systems, one would need to study a smaller class of problems and certain spectral bounds to derive a working rule for setting the parameters.

*Johnson-Lindenstrauss sketches.* In [21], the authors show how high dimension data can be projected onto a low dimension subspace using a Gaussian matrix in such a way that it approximately preserves the pairwise distance between points. This result is the celebrated JL lemma. Since then, many more random transforms have been shown to satisfy this property, such as Count sketch [9], subsampled Fourier [1], and Hadamard transforms [6, 44, 20]. These JL transforms have been used to speed up LAPACK solvers [4], solving linear regression [49] and Newton based methods [37].

*Contributions.* We propose new combinations of the sketch-and-project method with JL sketches such as the subsampled randomized Hadamard transform (SRHT) and a new subsampled count (SubCount) sketch. Our new package RidgeSketch is also set up in a way that is easily extensible, where new sketching methods can easily be added as a new instance of a *sketching class* (see Appendix B). Our results for the subsampled Hadamard sketch are negative. We show, despite the favorable theoretical complexity of using Hadamard sketches, that the overhead costs make them a completely impractical choice for sketch-and-project methods. SubCount sketch, on the other hand, when combined with sketch-and-project, results in an efficient method.

*Alternative iterative sketching based methods.* A method closely related to the sketch-and-project method is the iterative Hessian sketch [37], which makes use of iterative sketching to solve constrained quadratics. In the unconstrained setting such as (1.1), iterative Hessian sketch is efficient when the dimension *d* is significantly smaller than number of data points *n*. This rules out one of our main applications: kernel ridge regression where *n* = *d*.

**2.1. Linear system formulation.** Since the optimization problem (1.1) is differentiable and without constraints, its solution satisfies the stationarity conditions given by

$$ (2.1) \qquad (\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})\mathbf{w} = \mathbf{X}^\top\mathbf{y}. $$

We can also rewrite the above linear system in its *dual form* given by

$$ (2.2) \qquad \mathbf{w} = \mathbf{X}^\top\alpha, \quad \text{where} \quad (\mathbf{X}\mathbf{X}^\top + \lambda\mathbf{I})\alpha = \mathbf{y}. $$

The equivalence between solving (2.1) and (2.2) is well known and proven in the appendix in Lemma A.1 for completion. The choice of solving (2.1) or (2.2) will depend on the dimensions of the data. On the one hand, the system in (2.1) involves a $d \times d$ matrix, and thus the *primal form* (2.1) is preferred when $d \le n$. On the other hand, the system in (2.2) involves an $n \times n$ matrix, and thus solving the *dual form* (2.2) is preferred when $n < d$.

In either case, the bottleneck cost is the solution of a linear system where the system matrix is symmetric and positive definite. To simplify notation we introduce $\mathbf{A} := \mathbf{B} + \lambda\mathbf{I}$ and let

$$ (2.3) \qquad \mathbf{A}\mathbf{w} = \mathbf{b}, $$
where $\mathbf{B} = \mathbf{X}^{\top}\mathbf{X}$ and $b = \mathbf{X}^{\top}y$, for the primal form (2.1), or $\mathbf{B} = \mathbf{XX}^{\top}$ and $b = y$ for the dual one (2.2). Note that in the latter case, a multiplication by $\mathbf{X}^{\top}$ is required to recover the weights vector. Let $m$ be the dimensions of $\mathbf{A} \in \mathbb{R}^{m \times m}$; thus $m$ equals the smallest dimension of the design matrix $\mathbf{X} \in \mathbb{R}^{n \times d}$. Indeed, $m = d$ if we choose to solve the primal version (2.1) and $m = n$ if we choose to solve the dual one (2.2).

Thus solving the ridge regression problem boils down to finding the solution of the linear system (2.3). If the dimension of the squared matrix $\mathbf{A}$, denoted $m$, is not too large, one can solve this problem using a direct solver (for instance, through an SVD or an Cholesky decomposition). But when $m$ is large, direct methods become intractable as their computational cost grows with $O(m^3)$.

**2.2. Using a kernel.** We also consider kernel ridge regression, where the feature matrix is the result of applying a feature map. This leads to particular considerations since the resulting feature matrix may have an infinite number of columns.

The idea behind kernel ridge regression is that, instead of learning using the original input (or feature) vectors $\mathbf{X} := [x_1, \dots, x_n]$, we can learn using a high dimensional feature map of the inputs $\phi: \mathbb{R}^d \to \mathbb{R}^r$ where $r > d$ or even an infinite dimensional space. For instance, $\phi(x)$ could encode a high dimensional polynomial. By replacing each $x_i$ with $\phi(x_i)$ in (1.1) we arrive at

$$ (2.4) \qquad \min_{w \in \mathbb{R}^r} \frac{1}{2} \sum_{i=1}^n (\phi(x_i)w - y_i)^2 + \frac{\lambda}{2} \|w\|_2^2 . $$

When $r$ is large, or even infinite, solving (2.4) directly can be difficult or intractable. Fortunately, the dual formulation of (2.4) is always an $n$-dimensional problem independently of the dimension $r$.<fnref n="2" /> The dual formulation of (2.4) is given by

$$ (2.5) \qquad \alpha^* = \arg \min_{\alpha \in \mathbb{R}^n} \frac{1}{2} \|\mathbf{K}\alpha - y\|_2^2 + \frac{\lambda}{2} \alpha^\top \mathbf{K} \alpha, $$

where $\mathbf{K} = (\mathbf{K}(x_i, x_j))_{ij} := \langle \phi(x_i), \phi(x_j) \rangle_{ij}$ is the kernel matrix. This is equivalent to solving in $\alpha$ the linear system

$$ (2.6) \qquad \mathbf{K}(\mathbf{K} + \lambda\mathbf{I})\alpha = \mathbf{K}y. $$

With $\alpha^*$, the solution to the above, we can then predict the output of a new input vector $x$ using

$$ \text{predict}(x) = \sum_{i=1}^{n} \alpha_{i}^{*} \langle \phi(x_{i}), \phi(x) \rangle = \sum_{i=1}^{n} \alpha_{i}^{*} \mathbf{K}(x_{i}, x) = \mathbf{K}\alpha^{*}. $$

Consequently to solve (2.5) and make predictions, we only need access to the kernel matrix. Fortunately, there are several feature maps for which the kernel matrix is easily computable, including the one we use in our experiments, which is the Gaussian kernel, otherwise known as the radial basis function (RBF)

$$ (2.7) \qquad K(x, x') = \exp\left(\frac{-\|x - x'\|^2}{2\sigma^2}\right), $$

where $\sigma > 0$ is the kernel parameter.

<footnote n="2">²This is commonly known as the kernel trick; see Chapter 16 in [42].</footnote>
Ultimately, despite the addition of a kernel, the resulting problem (2.6) is still a linear system of the form $Aw = b$ where

$$A = K^{\top}(K + \lambda I) \quad \text{and} \quad b = K^{\top}y.$$

The only marked difference now is that $A$ tends to be dense, and because of this, matrix vector products are particularly expensive.

**3. The sketch-and-project method.** Sketch-and-project is an archetypal algorithm that unifies a variety of randomized iterative methods including both randomized Kaczmarz and CD [16], and all their block and importance sampling variants.

At each iteration, the sketch-and-project methods randomly compressed the linear system using what is known as a sketching matrix.

**DEFINITION 3.1.** Let $\tau \in \mathbb{N}$ and let $\mathcal{D}$ be a distribution over matrices in $\mathbb{R}^{m \times \tau}$. We refer to $\tau$ as the sketch size and to $S \in \mathbb{R}^{m \times \tau}$ drawn from the distribution $\mathcal{D}$ as a sketching matrix.

We can use a sketching matrix $S$ to reduce the number of rows of the linear system (2.3) to $\tau$ rows as follows:

$$(3.1) \qquad S^{\top}Aw = S^{\top}b.$$

If the sketch size $\tau$ is sufficiently large and the sketching matrix is appropriately chosen, then we can guarantee with high probability that the solution to the sketched linear system (3.1) is close to the solution $w^*$ of the original system (2.3) (see [29]). But this one-shot sketching approach poses several challenges: (1) it may be hard to determine how large $\tau$ should be, (2) the sketched linear system now has multiple solutions, and (3) with some low probability the solution to (3.1) could be far from $w^*$. To address these issues, we use an iterative projection scheme.

<caption>**Algorithm 3.1** The Sketch-and-Project method.</caption>

1: **parameters:** distribution over random $m \times \tau$ matrices in $\mathcal{D}$, tolerance $\epsilon > 0$
2: Set $w^0 = 0 \in \mathbb{R}^m$
3: Set $r^0 = Aw^0 - b = -b \in \mathbb{R}^m$
4: $k = 0$
5: **while** $\|r^k\|_2 / \|r^0\|_2 \le \epsilon$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;Sample an independent copy $S_k \sim \mathcal{D}$
7: &nbsp;&nbsp;&nbsp;&nbsp;$r_S^k = S_k^\top r^k$
8: &nbsp;&nbsp;&nbsp;&nbsp;$\delta_k = \text{least\_norm\_solution}(S_k^\top AS_k, r_S^k)$
9: &nbsp;&nbsp;&nbsp;&nbsp;$w^{k+1} = w^k - S_k \delta_k$
10: &nbsp;&nbsp;&nbsp;&nbsp;$r^{k+1} = r^k - AS_k \delta_k$
11: &nbsp;&nbsp;&nbsp;&nbsp;$k = k + 1$
12: **end while**
13: **return** $w^t$

$\triangleright$ weights initialization
$\triangleright$ residual initialization

$\triangleright$ compute sketched residual
$\triangleright$ solve sketched system
$\triangleright$ update the iterates
$\triangleright$ update the residual

$\triangleright$ return weights vector

Let $W$ be a symmetric positive definite matrix of order $n$ (which will typically be chosen as $A$); here we project with respect to the $W$-norm$^3$ given by $\|\cdot\|_W = \sqrt{\langle \cdot, W\cdot \rangle}$.

<footnote>$^3$Using this norm for symmetric positive definite matrices has shown to result in algorithms with a fast convergence rate [16].</footnote>
At the kth iteration of the sketch-and-project algorithm, a sketching matrix $S_k$ is drawn from $\mathcal{D}$ and the current iterate $w^k$ is projected onto the solution space of the sketched system $S_k^\top A x = S_k^\top b$ with respect to the $\mathbf{W}$-norm, that is,

$$ (3.2) \qquad w^{k+1} = \arg \min_{w \in \mathbb{R}^m} \|w - w^k\|_{\mathbf{W}}^2 \quad \text{subject to} \quad S_k^\top A w = S_k^\top b. $$

The closed form solution to (3.2) is given by

$$ (3.3) \qquad w^{k+1} = w^k - \mathbf{W}^{-1} \mathbf{A} S_k (S_k^\top \mathbf{A} \mathbf{W}^{-1} \mathbf{A} S_k)^\dagger S_k^\top (\mathbf{A} w^k - b), $$

where $\dagger$ denotes the pseudoinverse. When $\mathbf{A}$ is known to be positive definite, as in our case, using $\mathbf{W} = \mathbf{A}$ often results in an overall faster convergence of (3.3) as shown in [16]. Using $\mathbf{W} = \mathbf{A}$ in (3.3) gives the updates

$$ (3.4) \qquad w^{k+1} = w^k - S_k (S_k^\top \mathbf{A} S_k)^\dagger S_k^\top (\mathbf{A} w^k - b). $$

We refer to (3.4) as the RidgeSketch method since it is specialized for solving ridge regression. Here we give the details on how to efficiently implement the RidgeSketch update (3.4); see Algorithm 3.1.

One practical detail we have added to the pseudocode in Algorithm 3.1 is a stopping criteria. For any iterative algorithm, it is important to know when to stop. We can do this by monitoring the residual $r^k := \mathbf{A}w^k - b$. From an initial residual $r^0 \in \mathbb{R}^m$, when the relative residual $\|r^k\| / \|r^0\|$ is below a given tolerance, we stop. We also need this residual for computing the update (3.4). We can efficiently update the residual from one iteration to the next since

$$ (3.5) \qquad r^{k+1} = \mathbf{A} w^{k+1} - b = \mathbf{A} (w^k - \mathbf{S}_k \delta_k) - b = r^k - \mathbf{A} \mathbf{S}_k \delta_k, $$

where $\delta_k := (S_k^\top \mathbf{A} S_k)^\dagger S_k^\top (\mathbf{A} w^k - b) = (S_k^\top \mathbf{A} S_k)^\dagger S_k^\top r^k$. Thus we can update the residual at the cost of $O(m\tau)$, that is, multiplying a $m \times \tau$ matrix $\mathbf{A} S_k$ with the $\tau$ dimensional vector $\delta_k$. Note that $\delta_k$ can be efficiently computed as the least-norm solution of the following linear system in $x$

$$ (3.6) \qquad S_k^\top \mathbf{A} S_k x = S_k^\top r^k. $$

**4. Sketching methods and matrices.** Here we introduce several sketching matrices that can be used in Algorithm 3.1. The ideal sketch is one that reduces the dimension of the linear system (2.3) as much as possible, while preserving as much information as possible, and that can be efficiently implemented. As we discuss throughout this section, there is no sketch that has all three of these qualities, and ultimately, one must make a trade-off between them.

**4.1. Classical sketches.** One of the most classical and simple sketching methods is the Gaussian sketch.

**DEFINITION 4.1.** A Gaussian sketch is a random matrix $S \in \mathbb{R}^{m \times \tau}$ where each element is sampled independent and identically distributed (i.i.d.) for the standard Gaussian distribution.

As argued in [36], the resulting sketched matrix $S^\top A$ can be a good approximation of the full matrix and is easy to control with probability bounds. Though simple to implement, the cost of forming $S^\top A$ is $O(\tau m^2)$, which is expensive.

A much cheaper option is to use a *subsampling sketch*.
DEFINITION 4.2. A subsampling sketch is based on a randomly sampled subset $C \subset \{1, \dots, m\}$ with $|C| = \tau$ elements drawn uniformly on average from all such subsets. Let $\mathbf{I}_C \in \mathbb{R}^{m \times \tau}$ denote the concatenation of the columns of the identity matrix $\mathbf{I} \in \mathbb{R}^{m \times m}$ whose columns are indexed in $C$. We define the subsampling sketch distribution as

$$ (4.1) \qquad \mathbb{P}[\mathbf{S} = \mathbf{I}_C] = \frac{1}{\binom{m}{\tau}} \quad \text{for all } C \subset \{1, \dots, m\}, |C| = \tau. $$

Subsampling sketches are very cheap to compute; indeed, we need not even compute $\mathbf{S}^\top \mathbf{A}$ since it is simply equivalent to fetching the rows of $\mathbf{A}$ indexed by a random subset. This can be done in Python without generating any copies of the data by slicing the selected rows. Slicing is a very well optimized operation in NumPy [46] and SciPy [47] compressed sparse row (CSR) sparse arrays, which makes it one the fastest sketching methods.

Though subsampling sketches are cheap and fast, the sketched matrix $\mathbf{S}^\top \mathbf{A}$ can be a poor approximation of $\mathbf{A}$, since it is always possible that some vital part of $\mathbf{A}$ is “left out” in the rows that were not sampled. Still, the subsampling sketch will prove to work well within the iterative sketch-and-project scheme.

Next we consider a sketch that makes use of subsampling, addition, and subtraction of rows.

## 4.2. Count and SubCount sketch.
In order to avoid losing too much information by just subsampling rows, one can also sum and subtract groups of rows. This is the idea behind count sketch which stems from the streaming data literature [7, 10] and got popularized as a matrix sketching tool by [8]. Count sketch selects rows of $\mathbf{A}$, flips their sign with probability 1/2, and adds it to a random row, sampled uniformly, of the output matrix $\mathbf{S}^\top \mathbf{A}$.

To decrease the overall cost of Count sketch, we also combined it with a subsampling step. We call the resulting method SubCount sketch. SubCount sketch has two parameters, the subsampling size $s \in \{1, \dots, m\}$ and the sum size $k \in \{1, \dots, s\}$ that must be such that it divides $s$. SubCount sketch can be broken down into three steps: subsampling $s$ random rows of the input, then randomly flipping their sign and finally summing $k$ contiguous rows together. This can also be illustrated in terms of matrix multiplications as follows.

DEFINITION 4.3. A SubCount sketch is a matrix $\mathbf{S} \in \mathbb{R}^{m \times \tau}$ such that $\mathbf{S}^\top = \Sigma \mathbf{D} \mathbf{I}_C$, where

* $\mathbf{I}_C \in \mathbb{R}^{s \times m}$ is a subsampling matrix based on a set $C \subset \{1, \dots, m\}$ chosen uniformly at random from all sets with $s$ elements,

* $\mathbf{D} \in \mathbb{R}^{s \times s}$ is a diagonal matrix with elements sampled uniformly from $\{-1, 1\}$,

* $\Sigma \in \mathbb{R}^{k \times s}$ is a sum matrix that sums every $k$ contiguous rows together. It reads

$$ (4.2) \qquad \Sigma = \begin{bmatrix} \underbrace{1 \cdots 1}_{k} & 0 \cdots 0 & \cdots & 0 \cdots 0 \\ 0 \cdots 0 & \underbrace{1 \cdots 1}_{k} & \cdots & 0 \cdots 0 \\ 0 \cdots 0 & 0 \cdots 0 & \cdots & \underbrace{1 \cdots 1}_{k} \end{bmatrix}. $$

See Figure 2 for a depiction of the first two steps of SubCount sketch. The resulting sketch size is given by $\tau = s/k \in \mathbb{N}^*$. In our RidgeSketch package, since
<figure><img src="image_3.png" /><figcaption>FIG. 2. SubCount sketch example.</figcaption></figure>

the sketch size $\tau$ is a parameter selected by the user, we adjusted $s$ and $k$ such that $\tau = s/k$ as follows: If $10\tau \le m$, we set arbitrarily $k$ to 10; else, we fix $k = \lfloor m/\tau \rfloor$ and then compute $s = k\tau$. When we have no subsampling, that is, $s = m$, then we simply refer to SubCount sketch and count sketch.<fnref n="4" /> The advantage of count sketch is that it enables the computation of $S^\top A$ with a cost of $O(nnz(A))$, where $nnz(A)$ is the number of nonzero values of $A$. Thus, this sketching method is fast for large sparse matrices, especially for the CSR format for which accessing rows is efficient. Moreover, because it linearly combines groups of rows of the input matrix, it avoids the pitfall of leaving out meaningful rows of $A$, unlike subsampling. This will later be confirmed numerically in section 8.

**4.3. Subsampled randomized Hadamard transform.** Here we consider the SRHT [50, 44], which we refer to as the Hadamard sketch.

**DEFINITION 4.4.** A Hadamard sketch or SRHT is a matrix $S \in \mathbb{R}^{m \times \tau}$ such that

$$ (4.3) \qquad S^\top = \frac{1}{\sqrt{\tau m}} I_C H_m D, $$

where we assume there exists $q \in \mathbb{N}^*$ such that $m = 2^q$, and

* $D \in \mathbb{R}^{m \times m}$ is a diagonal matrix with elements sampled uniformly from $\{-1, 1\}$,

* $H_m \in \mathbb{R}^{m \times m}$ is the Hadamard matrix of order $m$ defined recursively through (4.4)

$$ H_1 = [1], \quad H_2 = \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}, \quad \dots, \quad H_{2^q} = \begin{bmatrix} H_{2^{q-1}} & H_{2^{q-1}} \\ H_{2^{q-1}} & -H_{2^{q-1}} \end{bmatrix} \in \mathbb{R}^{m \times m}, $$

* $I_C \in \mathbb{R}^{\tau \times m}$ is a subsampling matrix based on a set $C \subset \{1, \dots, m\}$ chosen uniformly at random from all sets with $\tau$ elements

The Hadamard sketch has recently become popular in several applications in numerical linear algebra because it satisfies the JL lemma [6] and it can be computed

<footnote n="4">⁴The standard definition of count sketch also shuffles the columns of the summing matrix (4.2). We did not use this shuffling since we found that it had little to no effect on performance in our setting.</footnote>
using $O(m^2 \log \tau)$ operations using the fast Walsh-Hadamard transform [13] or the recent trimmed Walsh-Hadamard transform [2]. Let $\mathbf{H}$ be the Hadamard matrix defined in (4.4) of order $m$. Such methods use a butterfly structure, like the Cooley-Tukey FFT algorithm, to compute $\mathbf{H}y$ in $O(m \log m)$ time for any $y \in \mathbb{R}^m$, but they require $m$ to be a power of 2.

This last detail is often glossed over in theory, but we find in practice this poses a challenge when $m$ is large. One can remove some rows of $\mathbf{A}$ to meet this assumption, with the risk of losing vital information. Instead, when $m$ is not a power of 2, we need to pad $\mathbf{A}$ with zero rows until the augmented matrix has $m'$ rows, where $m'$ is a power of 2. That is, we need $m' = 2^{\lceil \log_2 m \rceil}$ rows. In the worst-case scenario, when $m = 2^q + 1$ for some $q \in \mathbb{N}^*$, then $m' = 2^{q+1}$, effectively doubling the number of rows.

**5. Convergence theory.** Here we present theoretical convergence guarantees of the RidgeSketch method (3.4) for sketch matrices $\mathbf{S} \in \mathbb{R}^{m \times \tau}$ drawn from a fixed distribution $\mathcal{D}$. Later on, in section 6 we provide a new convergence theory for a new momentum variant of the sketch-and-project method. To contextualize our contribution, we will first present the previously known convergence theory of the sketch-and-project method.

**5.1. Convergence of the iterates.** The sketch-and-project method (3.4) enjoys a linear convergence in $L^2$ given as follows.

**THEOREM 5.1** (convergence of the iterates [17, 16]). Let $w^*$ be a solution of (2.3) and let $w^0 \in \text{Range}(\mathbf{A})$. Let $t \in \mathbb{N}$ and consider the iterates given in (3.3). It follows that

$$ (5.1) \qquad \mathbb{E} \left[ \|w^t - w^*\|_{\mathbf{A}}^2 \right] \le (1-\rho)^t \|w^0 - w^*\|_{\mathbf{A}}^2, $$

where

$$ (5.2) \qquad \rho = \lambda_{\min}^{+}(\mathbf{A}^{1/2}\mathbb{E}[\mathbf{S}(\mathbf{S}^{\top}\mathbf{A}\mathbf{S})^{\dagger}\mathbf{S}]^{1/2}) . $$

Consequently, we also have for the residual that

$$ (5.3) \qquad \mathbb{E} \left[ \| \mathbf{A} w^t - b \|^2 \right] \le (1-\rho)^t \lambda_{\max}(\mathbf{A}) \| w^0 - w^* \|_{\mathbf{A}}^2. $$

*Proof.* The proof of (5.1) is given in Theorem 1.1 in [17]. The convergence of the residual in (5.3) follows as a consequence since

$$ (5.4) \qquad \| \mathbf{A} w^t - b \|^2 = \| w^t - w^* \|_{\mathbf{A}^2}^2 \le \lambda_{\max}(\mathbf{A}) \| w^t - w^* \|_{\mathbf{A}}^2. \quad \square $$

This linear convergence in $L^2$ is often thought of as a gold standard for convergence of stochastic sequences, since it implies convergence in high probability and of all the moments of the sequence. Furthermore, the error decays at an exponential rate determined by $\rho$. Yet the downside of (5.1) and (5.3) is that the rate of convergence $\rho$ can be very small. Next we present a sublinear rate of convergence in $L^2$ that has an improved rate of convergence.

**5.2. Convergence of the residuals with step size $0 < \gamma < 1$.** The convergence proofs we present next rely on viewing the sketch-and-project method as an instance of SGD (stochastic gradient descent) [39, 52]. To establish this SGD viewpoint, we first reformulate the problem of solving (2.3) as the following minimization problem:

$$ (5.5) \qquad \min_{w \in \mathbb{R}^d} f(w) := \frac{1}{2} \| \mathbf{A}w - b \|_{\mathbb{E}[\mathbf{H}\mathbf{S}]}^2 = \mathbb{E} \left[ \frac{1}{2} \| \mathbf{A}w - b \|_{\mathbf{H}\mathbf{S}}^2 \right] =: \mathbb{E}[f_{\mathbf{S}}(w)], $$
where

$$ (5.6) \qquad \mathbf{H}_{\mathbf{S}} := \mathbf{S}(\mathbf{S}^{\top}\mathbf{A}\mathbf{W}^{-1}\mathbf{A}\mathbf{S})^{\dagger}\mathbf{S}^{\top}. $$

Solving our linear systems is now equivalent to solving the stochastic minimization problem (5.5), for which the classic method is SGD. It just so happens that SGD with a step size of $\gamma = 1$ is exactly the sketch-and-project iteration (3.3). Indeed, if we compute the gradients with respect to the $\|w\|_{\mathbf{W}} = \sqrt{\langle \mathbf{W}w, w \rangle}$ norm, then SGD is given by

$$ (5.7) \qquad w^{k+1} = w^k - \gamma \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k), $$

where the gradient relative to the weighted inner product $\langle \cdot, \cdot \rangle_{\mathbf{W}}$ is given by

$$ (5.8) \qquad \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k) := \mathbf{W}^{-1} \mathbf{A} \mathbf{H}_{\mathbf{S}} (\mathbf{A} w^k - b), $$

where $\mathbf{S} \sim \mathcal{D}$ is sampled i.i.d. at each iteration and $\gamma > 0$ is the step size. Since $f_{\mathbf{S}}(w)$ is an unbiased estimate of $f(w)$ we have that

$$ (5.9) \qquad \mathbb{E}[\nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k)] = \nabla_{\mathbf{W}} f(w^k). $$

Using the interpretation as SGD, we can provide the convergence of $\mathbb{E}[f(w^k)]$ to zero. We can then relate the convergence of $f(w^k)$ to the convergence of the residual using the following lemma.

LEMMA 5.2.

$$ (5.10) \qquad \lambda_{\min}(\mathbb{E}[\mathbf{H}_{\mathbf{S}}]) \|\mathbf{A}w - b\|^2 \le 2f(w). $$

*Proof.* This follows from

$$ \lambda_{\min}(\mathbb{E}[\mathbf{H}_{\mathbf{S}}]) \|\mathbf{A}w - b\|^2 = \|\mathbf{A}w - b\|^2_{\lambda_{\min}(\mathbb{E}[\mathbf{H}_{\mathbf{S}}])\mathbf{I}_m} \le \|\mathbf{A}w - b\|^2_{\mathbb{E}[\mathbf{H}_{\mathbf{S}}]} = 2f(w). \quad \square $$

First we need the following property of the gradient taken from Lemma 3.1 in [39].

LEMMA 5.3 (gradient norm—function identity). It follows that

$$ (5.11) \qquad \|\nabla_{\mathbf{W}} f_{\mathbf{S}}(w)\|^2_{\mathbf{W}} = 2f_{\mathbf{S}}(w). $$

*Proof.* For completeness we give the proof. By straightforward computation we have that

$$ \begin{align*} \|\nabla_{\mathbf{W}} f_{\mathbf{S}}(w)\|^2_{\mathbf{W}} & \stackrel{(5.8)}{=} \langle \mathbf{W}^{-1} \mathbf{A} \mathbf{H}_{\mathbf{S}} (\mathbf{A}w - b), \mathbf{W}^{-1} \mathbf{A} \mathbf{H}_{\mathbf{S}} (\mathbf{A}w - b) \rangle_{\mathbf{W}} \\ & = (\mathbf{A}w - b)^{\top} \mathbf{H}_{\mathbf{S}} \mathbf{A} \mathbf{W}^{-1} \mathbf{A} \mathbf{H}_{\mathbf{S}} (\mathbf{A}w - b) \\ & \stackrel{(5.6)}{=} (\mathbf{A}w - b)^{\top} \mathbf{H}_{\mathbf{S}} (\mathbf{A}w - b) \\ & \stackrel{(5.5)}{=} 2f_{\mathbf{S}}(w), \end{align*} $$

where in the third equality we expanded $\mathbf{H}_{\mathbf{S}}$ and applied the identity $\mathbf{M}^{\dagger} = \mathbf{M}^{\dagger}\mathbf{M}\mathbf{M}^{\dagger}$ with $\mathbf{M} = \mathbf{S}^{\top}\mathbf{A}\mathbf{W}^{-1}\mathbf{A}\mathbf{S}$. $\square$

Furthermore, the functions $f_{\mathbf{S}}$ are convex.

LEMMA 5.4. The function $f_{\mathbf{S}}$ is a convex quadratic. Consequently

$$ (5.12) \qquad f_{\mathbf{S}}(y) \ge f_{\mathbf{S}}(x) + \langle \nabla_{\mathbf{W}} f_{\mathbf{S}}(x), y - x \rangle_{\mathbf{W}}. $$
Proof. Letting $w \in \mathbb{R}^d$ and $S \sim D$, the gradient relative to the Euclidean inner product $\langle \cdot, \cdot \rangle$ is

$$ \nabla f_S(x) = A H_S (Ax - b) , $$

and thus the Hessian is

$$ \nabla^2 f_S(x) = A H_S A , $$

which is a semidefinite positive matrix. This implies that $f_S$ is convex. As a consequence (5.12) holds since

$$ \begin{aligned} f_S(y) &\ge f_S(x) + \langle \nabla f_S(x), y - x \rangle \\ &= f_S(x) + \langle \nabla_W f_S(x), y - x \rangle_W \end{aligned} $$

and given that $\langle \nabla_W f_S(x), y - x \rangle_W = \langle WW^{-1} \nabla f_S(x), y - x \rangle = \langle \nabla f_S(x), y - x \rangle$. $\square$

Next we establish the sublinear convergence of the *average of the iterates* of the sketch-and-project method. This result is a direct consequence of Theorem 4.10 in [39] and has already been proven in Theorem 3 in [27]. We present the complete statement and proof since it is a warm-up for our forthcoming results, and since the proof is substantially simpler than the one presented in [27].

**THEOREM 5.5** (convergence of the residuals). Let $\gamma \in (0, 1)$. Let $w^*$ be a solution of (2.3), let $t \in \mathbb{N}^*$, and let $(w^k)_{0 \le k \le t}$ be the iterates given by (5.7). It follows that

$$ (5.13) \qquad f(\bar{w}^t) \le \frac{1}{t} \sum_{k=0}^{t} f(w^k) \le \frac{1}{t} \frac{\|w^0 - w^*\|_{\mathbf{W}}^2}{2\gamma(1-\gamma)}, $$

where $\bar{w}^t := \frac{1}{t} \sum_{k=0}^{t} w^k$. Furthermore

$$ (5.14) \qquad \| \mathbf{A} \bar{w}^t - b \|^2 \le \frac{1}{t} \frac{1}{\lambda_{\min}(\mathbb{E}[\mathbf{H}_S])} \frac{\| w^0 - w^* \|_{\mathbf{W}}^2}{2\gamma(1-\gamma)}. $$

*Proof.*

$$ \begin{align*} \|w^{k+1} - w^*\|_{\mathbf{W}}^2 &= \|w^k - w^*\|_{\mathbf{W}}^2 - 2\gamma \langle \nabla_{\mathbf{W}} f_S(w^k), w^k - w^* \rangle_{\mathbf{W}} \\ &\quad + \gamma^2 \|\nabla_{\mathbf{W}} f_S(w^k)\|_{\mathbf{W}}^2 \\ &\stackrel{(5.11)}{=} \|w^k - w^*\|_{\mathbf{W}}^2 - 2\gamma \langle \nabla_{\mathbf{W}} f_S(w^k), w^k - w^* \rangle_{\mathbf{W}} + 2\gamma^2 f_S(w^k) \\ &\stackrel{(5.12)}{\le} \|w^k - w^*\|_{\mathbf{W}}^2 - 2\gamma(1-\gamma) f_S(w^k) . \end{align*} $$

Rearranging, and dividing through by $\gamma(1-\gamma) > 0$ we have that

$$ f_S(w^k) \le \frac{1}{2\gamma(1-\gamma)} \left( \|w^k - w^*\|_{\mathbf{W}}^2 - \|w^{k+1} - w^*\|_{\mathbf{W}}^2 \right) . $$

Summing up over both sides for $k = 0, \dots, t$ and using telescopic cancellation we have that

$$ (5.15) \qquad \sum_{k=0}^{t} f_S(w^k) \le \frac{1}{2\gamma(1-\gamma)} \left( \|w^0 - w^*\|_{\mathbf{W}}^2 - \|w^{t+1} - w^*\|_{\mathbf{W}}^2 \right) \le \frac{\|w^0 - w^*\|_{\mathbf{W}}^2}{2\gamma(1-\gamma)} . $$
By applying the Jensen's inequality to $f_S$, which is convex since it is a quadratic form,

$$ (5.16) \qquad f_S \left( \frac{1}{t} \sum_{k=0}^{t} w^k \right) \le \frac{1}{t} \sum_{k=0}^{t} f_S(w^k). $$

Now, by denoting $\bar{w}^t := \frac{1}{t} \sum_{k=0}^{t} w^k$, the convergence result (5.13) follows by dividing (5.15) by $t$ and using (5.16). Finally, the convergence of the residual in (5.14) follows from (5.10). □

The weakness of Theorem 5.5 is that it describes how the average of the iterates $\bar{w}^t$ converge instead of the last one $w^t$. This type of convergence is problematic since $\bar{w}^t$ gives as much importance, a $1/t$ weight, to the initial point $w^0$ as to the last one, $w^t$. Since $w^0$ is often chosen arbitrarily, the average of the iterates can converge only as fast as $w^0$ is forgotten. That is, the average cannot converge faster than $1/t$. This is apparent in experiments, where using averaging from the start results in a slow convergence. In practice, averaging only the last few iterates works substantially better, but it is not supported in theory. To resolve this issue, we will replace this equal averaging with a weighted average that gives more weight to recent iterates and forgets the initial conditions exponentially fast.

**6. Momentum.** A common variant of SGD is to add momentum. Since the sketch-and-project method can be interpreted as SGD (5.7), we can add momentum. Let $\gamma_k \in [0, 1]$ and $\beta_k \in [0, 1]$ be respectively the step size and the momentum parameter. The heavy ball formulation of momentum is given by

$$ (6.1) \qquad w^{k+1} = w^k - \gamma_k \nabla_W f_S(w^k) + \beta_k (w^k - w^{k-1}). $$

Note that we have now allowed for a step size $\gamma_k$ that is iteration dependent.

This same heavy ball formulation was considered in [27], where the authors propose a precise analysis and show no benefit using momentum. But, all of their analysis assumes that the momentum parameter is constant. It turns out that by allowing $\beta_k$ to be iteration dependent, we can do better. But first, we need the *iterative averaging viewpoint* of momentum.

**6.1. Iterate averaging viewpoint.** Recently, a new *iterate averaging* parametrization of the momentum method was proposed in [41]. This iterative averaging parametrization is given by

$$ (6.2) \qquad z^k = z^{k-1} - \eta_k \nabla_W f_S(w^k), $$

$$ (6.3) \qquad w^{k+1} = \left(1 - \frac{1}{\zeta_{k+1} + 1}\right) w^k + \frac{1}{\zeta_{k+1} + 1} z^k, $$

where we have introduced two new parameter sequences $\eta_k$ and $\zeta_k$ that map back to the $\gamma_k$ and $\beta_k$ parameters via

$$ (6.4) \qquad \gamma_k = \frac{\eta_k}{\zeta_{k+1} + 1} \quad \text{and} \quad \beta_k = \frac{\zeta_k}{\zeta_{k+1} + 1}. $$

Next we show that (6.3) produces the same $w_k$ iterates as (6.1).

**PROPOSITION 6.1** (equivalent formulations). Let $t \in \mathbb{N}^*$. The steps given in (6.1) and in (6.2)-(6.3) generate the same iterates $(w_k)_{0 \le k \le t}$.
*Proof.* First, let us expand the *iterate averaging* parametrization

$$ 
\begin{align} 
w^{k+1} & \stackrel{(6.3)}{=} \left(1 - \frac{1}{\zeta_{k+1} + 1}\right) w^k + \frac{1}{\zeta_{k+1} + 1} z^k \nonumber \\ 
& \stackrel{(6.2)}{=} w^k - \frac{1}{\zeta_{k+1} + 1} w^k + \frac{1}{\zeta_{k+1} + 1} \left(z^{k-1} - \eta \nabla_{\mathbf{w}} f_{\mathbf{s}}(w^k)\right) \nonumber \\ 
& = w^k - \frac{\eta}{\zeta_{k+1} + 1} \nabla_{\mathbf{w}} f_{\mathbf{s}}(w^k) + \frac{1}{\zeta_{k+1} + 1} \left(z^{k-1} - w^k\right) \nonumber \\ 
& \stackrel{(6.5)}{=} w^k - \gamma_k \nabla_{\mathbf{w}} f_{\mathbf{s}}(w^k) + \frac{1}{\zeta_{k+1} + 1} \left(z^{k-1} - w^k\right), 
\end{align} 
$$

where in the last step we used the left-hand-side relation in (6.4). To conclude the proof, we need only show that

$$ \frac{1}{\zeta_{k+1} + 1} (z^{k-1} - w^k) = \beta_k(w_k - w_{k-1}). $$

Indeed, this follows $w^k = (1 - \frac{1}{\zeta_{k+1}})w^{k-1} + \frac{1}{\zeta_{k+1}}z^{k-1}$, which can be rearranged in $z^{k-1} - w^k = \zeta_k(w^k - w^{k-1})$; consequently

$$ \frac{1}{\zeta_{k+1} + 1} (z^{k-1} - w^k) = \frac{\zeta_k}{\zeta_{k+1} + 1} (w^k - w^{k-1}) \stackrel{(6.4)}{=} \beta_k(w_k - w_{k-1}), $$

which together with (6.5) concludes the proof. □

Next we show how to leverage the iterative averaging viewpoint (6.2)–(6.3) to prove the convergence of the last iterates of sketch-and-project with momentum (6.1).

## 6.2. Convergence theorem.

**THEOREM 6.2.** Consider the iterates (6.2)–(6.3) with $w^0 = w^{-1}$. Let $\eta_k$ be a sequence of parameters with $0 < \eta_k < 1$. If

$$ (6.6) \qquad \zeta_0 = 0 \quad \text{and} \quad \zeta_k = \frac{1}{\eta_k} \sum_{t=0}^{k-1} \eta_t (1 - \eta_t) \quad \text{for all } k \ge 1, $$

then

$$ (6.7) \qquad \| \mathbf{A} w^k - b \|^2 \le \frac{1}{\lambda_{\min}(\mathbb{E}[\mathbf{H}_{\mathbf{s}}])} \frac{\| w^0 - w^* \|^2_{\mathbf{W}}}{\sum_{t=0}^k \eta_t (1 - \eta_t)} . $$

*Proof.* This proof is based on Theorem 3.1 in [41] for SGD applied to convex and smooth functions. Consider the Lyapunov function

$$ (6.8) \qquad L_k = \mathbb{E} \left[ \|z^k - w^*\|^2_{\mathbf{W}} \right] + 2\eta_k \zeta_k \mathbb{E} [f(w^k)] . $$

First note that

$$ 
\begin{align*} 
\|z^k - w^*\|^2_{\mathbf{W}} & \stackrel{(6.2)}{=} \|z^{k-1} - \eta_k \nabla_{\mathbf{w}} f_{\mathbf{s}}(w^k) - w^*\|^2_{\mathbf{W}} \\ 
& = \|z^{k-1} - w^*\|^2_{\mathbf{W}} - 2\eta_k \langle z^{k-1} - w^*, \nabla_{\mathbf{w}} f_{\mathbf{s}}(w^k) \rangle_{\mathbf{W}} \\ 
& \quad + \eta_k^2 \|\nabla_{\mathbf{w}} f_{\mathbf{s}}(w^k)\|^2_{\mathbf{W}} 
\end{align*} 
$$
$$
\begin{align*}
& \stackrel{(6.3)}{=} \|z^{k-1} - w^*\|_{\mathbf{W}}^2 - 2\eta_k \langle w^k - w^* + \zeta_k(w^k - w^{k-1}, \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k)) \rangle_{\mathbf{W}} \\
& \qquad + \eta_k^2 \|\nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k)\|_{\mathbf{W}}^2 \\
&= \|z^{k-1} - w^*\|_{\mathbf{W}}^2 - 2\eta_k \langle w^k - w^*, \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k) \rangle_{\mathbf{W}} \\
& \qquad - 2\eta_k \zeta_k \langle w^k - w^{k-1}, \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k) \rangle_{\mathbf{W}} + \eta_k^2 \|\nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k)\|_{\mathbf{W}}^2 \\
& \stackrel{(5.11)}{=} \|z^{k-1} - w^*\|_{\mathbf{W}}^2 + 2\eta_k^2 f_{\mathbf{S}}(w^k) - 2\eta_k \langle w^k - w^*, \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k) \rangle_{\mathbf{W}} \\
& \qquad - 2\eta_k \zeta_k \langle w^k - w^{k-1}, \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k) \rangle_{\mathbf{W}} .
\end{align*}
$$

The above holds to equality. Now we introduce the first inequality by calling upon
(5.12) so that

$$
\begin{equation} \tag{6.9}
\begin{aligned}
\|z^k - w^*\|_{\mathbf{W}}^2 &\stackrel{(5.12)}{\le} \|z^{k-1} - w^*\|_{\mathbf{W}}^2 + 2\eta_k^2 f_{\mathbf{S}}(w^k) - 2\eta_k f_{\mathbf{S}}(w^k) \\
&\quad - 2\eta_k \zeta_k (f_{\mathbf{S}}(w^k) - f_{\mathbf{S}}(w^{k-1})) \\
&= \|z^{k-1} - w^*\|_{\mathbf{W}}^2 - 2\eta_k(1 + \zeta_k - \eta_k)f_{\mathbf{S}}(w^k) + 2\eta_k \zeta_k f_{\mathbf{S}}(w^{k-1})\,.
\end{aligned}
\end{equation}
$$

The restriction on the parameters in (6.6) was designed so that

$$
(6.10) \qquad \eta_k(1 + \zeta_k - \eta_k) = \eta_{k+1} \zeta_{k+1} .
$$

Indeed, from (6.6) we have that

$$
\zeta_{k+1} = \frac{1}{\eta_{k+1}} \sum_{t=0}^{k} \eta_t (1 - \eta_t) = \frac{\eta_k}{\eta_{k+1}} \frac{1}{\eta_k} \sum_{t=0}^{k} \eta_t (1 - \eta_t) = \frac{\eta_k}{\eta_{k+1}} (\zeta_k + 1 - \eta_k) .
$$

Using (6.10) in (6.9) and taking expectation gives

$$
\mathbb{E} \left[ \|z^k - w^*\|_{\mathbf{W}}^2 \right] + 2\eta_{k+1}\zeta_{k+1}\mathbb{E} \left[ f(w^k) \right] \le \mathbb{E} \left[ \|z^{k-1} - w^*\|_{\mathbf{W}}^2 \right] + 2\eta_k\zeta_k\mathbb{E} \left[ f(w^{k-1}) \right] .
$$

Summing up both sides from $k=0, \dots, t$ and using telescopic cancellation gives

$$
\mathbb{E} \left[ \|z^t - w^*\|_{\mathbf{W}}^2 \right] + 2\eta_{t+1}\zeta_{t+1}\mathbb{E} \left[ f(w^t) \right] \le \mathbb{E} \left[ \|z^{t-1} - w^*\|_{\mathbf{W}}^2 \right] + 2\eta_0\zeta_0\mathbb{E} \left[ f(w^{t-1}) \right] .
$$

Using that $z^{-1} = w^0$ from (6.3), $\zeta_0 = 0$, and rearranging gives

$$
\mathbb{E} \left[ \|z^t - w^*\|_{\mathbf{W}}^2 \right] + 2\eta_{t+1}\zeta_{t+1} \mathbb{E} \left[ f(w^t) \right] \le \mathbb{E} \left[ \|w^0 - w^*\|_{\mathbf{W}}^2 \right] .
$$

Substituting the definition of $\zeta_{t+1}$ from (6.6) gives

$$
2\mathbb{E}[f(w^t)] \le \frac{1}{\eta_{t+1}\zeta_{t+1}} \mathbb{E}[\|w^0 - w^*\|_{\mathbf{W}}^2] = \frac{1}{\sum_{k=0}^{t} \eta_k (1 - \eta_k)} \mathbb{E}[\|w^0 - w^*\|_{\mathbf{W}}^2].
$$

The final step is a result of using (5.10) to lower bound $\mathbb{E}[f(w^t)]$.

<figure><img src="image_4.png" /></figure>

Theorem 6.2 shows that the convergences of the iterates $w^t$ depend on a sequence of parameters $\eta_k$. Next we give a corollary that shows that by simply choosing $\eta_k \equiv \eta < 1$ the iterates $w^t$ enjoy a fast sublinear convergence.
COROLLARY 6.3. Consider the setting of Theorem 6.2. Let $\eta_k \equiv \eta < 1$ and thus

$$ (6.11) \qquad \zeta_0 = 0 \quad \text{and} \quad \zeta_k = k(1-\eta) \quad \text{for all } k \ge 1, $$

and

$$ (6.12) \qquad \gamma_k = \frac{\eta}{(k+1)(1-\eta)+1} \quad \text{and} \quad \beta_k = 1 - \frac{2-\eta}{(k+1)(1-\eta)+1}, $$

then

$$ \|Aw^k - b\|^2 \le \frac{1}{\lambda_{\min}(\mathbb{E}[\mathbf{H}_S])} \frac{1}{k} \frac{\|w^0 - w^*\|_{\mathbf{W}}^2}{\eta(1-\eta)}. $$

Consequently, for $\eta = \frac{1}{2}$ and for a given tolerance $\epsilon > 0$, the iteration complexity of minimizing the residual is given by

$$ (6.13) \qquad t \ge \frac{4 \|w^0 - w^*\|_{\mathbf{W}}^2}{\lambda_{\min}(\mathbb{E}[\mathbf{H}_S])} \frac{1}{\epsilon} $$

to reach a desired precision

$$ \mathbb{E}[\|Aw^t - b\|^2] < \epsilon. $$

Corollary 6.3 shows that the last iterate of sketch-and-project with momentum converges at the same rate as the average of the iterates of sketch-and-project (see Theorem 5.5). This is the first tangible theoretical advantage of using momentum in this setting.

In section 7, through a more specialized setting for ridge regression, we show that the complexity given in (6.13) can be even tighter than the previously known linear convergence in Theorem 5.1. But first, we present a practical pseudocode implementation of sketch-and-project with momentum.

**6.3. Implementation.** In Algorithm 6.1 we have the pseudocode of sketch-and-project with momentum (6.1) for solving ridge regression (1.1) where $\mathbf{W} = \mathbf{A}$.

The residual $r^k := \|Aw^k - b\|^2$ is required for computing the update like in Algorithm 3.1. Fortunately we can efficiently compute the residual at step $k+1$ by storing the residuals at two steps $k-1$ and $k$ since

$$ (6.14) \qquad \begin{aligned} r^{k+1} &= \mathbf{A}w^{k+1} - b \\ &\stackrel{(6.1)}{=} \mathbf{A}((1 + \beta_k)w^k - \beta_k w^{k-1} - \gamma_k \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k)) - b \\ &\stackrel{(5.8)}{=} \mathbf{A}((1 + \beta_k)w^k - \beta_k w^{k-1} - \gamma_k \mathbf{S}_k \delta_k) - b \\ &\stackrel{(5.8)}{=} (1 + \beta_k)r^k - \beta_k r^{k-1} - \gamma_k \mathbf{A} \mathbf{S}_k \delta_k, \end{aligned} $$

where

$$ \delta_k = (\mathbf{S}_k^\top \mathbf{A} \mathbf{S}_k)^\dagger \mathbf{S}_k^\top r^k = \text{least\_norm\_solution} (\mathbf{S}_k^\top \mathbf{A} \mathbf{S}_k, r_{\mathbf{S}}^k), $$

and where we used that, since $\mathbf{W} = \mathbf{A}$, we have

$$ \nabla_{\mathbf{W}} f_{\mathbf{S}}(w^k) \stackrel{(5.8)}{=} \mathbf{W}^{-1} \mathbf{A} \mathbf{S} (\mathbf{S}^\top \mathbf{A} \mathbf{W}^{-1} \mathbf{A} \mathbf{S})^\dagger \mathbf{S}^\top (\mathbf{A} w^k - b) = \mathbf{S}_k (\mathbf{S}_k^\top \mathbf{A} \mathbf{S}_k)^\dagger \mathbf{S}_k^\top r^k. $$

Thus, for the momentum version of our algorithm we can keep the residual update at the cost of $O(m\tau)$ by just storing the residual $r^{k-1}$ at the previous time step.
**Algorithm 6.1** The sketch-and-project method with momentum.

1: **parameters:** distribution over random $m \times \tau$ matrices in $\mathcal{D}$, tolerance $\epsilon > 0$, momentum parameters $\eta_k, \zeta_k \in [0, 1]$

2: Set $w^{-1} = w^0 = 0 \in \mathbb{R}^m$ ▷ weights initialization

3: Set $r^{-1} = r^0 = \mathbf{A}w^0 - b = -b \in \mathbb{R}^m$ ▷ residual initialization

4: $k = 0$

5: **while** $\|r^k\|_2 / \|r^0\|_2 \le \epsilon$ **do**

6: $\gamma_k = \frac{\eta_k}{1+\zeta_{k+1}}$ and $\beta_k = \frac{\zeta_k}{1+\zeta_{k+1}}$

7: Sample an independent copy $\mathbf{S}_k \sim \mathcal{D}$

8: $r_{\mathbf{S}}^k = \mathbf{S}_k^\top r^k$ ▷ compute sketched residual

9: $\delta_k = \text{least\_norm\_solution} (\mathbf{S}_k^\top \mathbf{A} \mathbf{S}_k, r_{\mathbf{S}}^k)$ ▷ solve sketched system

10: $w^{k+1} = (1 + \beta_k)w^k - \beta_k w^{k-1} - \gamma_k \mathbf{S}_k \delta_k$ ▷ update the iterates

11: $r^{k+1} = (1 + \beta_k)r^k - \beta_k r^{k-1} - \gamma_k \mathbf{A} \mathbf{S}_k \delta_k$ ▷ update the residual

12: $k = k + 1$

13: **end while**

14: **return** $w^t$ ▷ return weights vector

**7. Specialized convergence theory for single column sketches.** Here we take a closer look at the rates of convergence given by Theorem 5.1 and Corollary 6.3 by considering a specialized setting of ridge regression ($\mathbf{W} = \mathbf{A}$) and single column sketches. That is, in this section we use a discrete distribution for $\mathbf{S}$ given by

$$ (7.1) \qquad \mathbb{P}[\mathbf{S} = s_i] = p_i \quad \text{for } i = 1, \dots, m, $$

where $s_1, \dots, s_m \in \mathbb{R}^m$ are a fixed collection of vectors and $\sum_{i=1}^q p_i = 1$.

To better understand the complexity (6.13) and (5.4), we first need to find a lower bound for $\lambda_{\min}(\mathbb{E}[\mathbf{H}_{\mathbf{S}}])$ where

$$ (7.2) \qquad \mathbb{E}[\mathbf{H}_{\mathbf{S}}] = \mathbb{E}[\mathbf{S}(\mathbf{S}^{\top}\mathbf{A}\mathbf{S})^{\dagger}\mathbf{S}^{\top}] . $$

By using a special choice for the probabilities in (7.1), we are able to give a convenient lower bound for $\mathbb{E}[\mathbf{H}_{\mathbf{S}}]$ in the following lemma.

**LEMMA 7.1.** Let $\mathbf{S}$ have a discrete distribution according to

$$ (7.3) \qquad \mathbb{P}[\mathbf{S} = s_i] = \frac{s_i^\top \mathbf{A} s_i}{\sum_{j=1}^q s_j^\top \mathbf{A} s_j} \quad \text{for } i = 1, \dots, q, $$

where $s_1, \dots, s_q \in \mathbb{R}^m$ are unit column vectors. Let $\mathbf{F} := [s_1, \dots, s_q] \in \mathbb{R}^{m \times q}$. It follows that

$$ (7.4) \qquad \lambda_{\min}(\mathbb{E}[\mathbf{H}_{\mathbf{S}}]) = \frac{\lambda_{\min}(\mathbf{F}\mathbf{F}^{\top})}{\sum_{j=1}^{m} s_j^{\top} \mathbf{A} s_j} . $$

*Proof.* This convenient probability distribution (7.3) was already considered in [16] in section 5.2. But there, the authors used these probabilities to study a different spectral quantity; thus for completion we adapt their proof to our setting. First note that

$$ \mathbb{E}[\mathbf{H}_{\mathbf{S}}] \stackrel{(7.2)}{=} \sum_{i=1}^{q} \frac{p_i}{s_i^{\top} \mathbf{A} s_i} s_i s_i^{\top} \stackrel{(7.3)}{=} \frac{\sum_{i=1}^{q} s_i s_i^{\top}}{\sum_{j=1}^{q} s_j^{\top} \mathbf{A} s_j} . $$
Consequently, the smallest eigenvalue is given by

$$ \lambda_{\min}(\mathbb{E}[\mathbf{H}_S]) = \lambda_{\min}\left(\frac{\sum_{i=1}^{q} s_i s_i^{\top}}{\sum_{j=1}^{q} s_j^{\top} \mathbf{A} s_j}\right) = \frac{\lambda_{\min}(\mathbf{F}\mathbf{F}^{\top})}{\sum_{j=1}^{m} s_j^{\top} \mathbf{A} s_j} . \quad \square $$

### 7.1. Comparing the complexity of CD with and without momentum.

Here we compare the fast sublinear convergence of RidgeSketch with momentum given in Theorem 6.2 to the linear convergence of RidgeSketch without momentum given in Theorem 5.1. Though linear convergence is generally preferred, we will show here that our new sublinear rate of convergence can be faster. To illustrate this, we will focus on the special case of CD.

The CD method is the result of applying RidgeSketch when the sketching matrices are unit coordinate vectors. That is, when $s_i = e_i \in \mathbb{R}^m$ is the $i$th column of the identity matrix and

$$ (7.5) \qquad \mathbb{P}[s_i = e_i] = \frac{\mathbf{A}_{ii}}{\text{tr}(\mathbf{A})} \quad \text{for } i = 1, \dots, m. $$

This particular nonuniform sampling in (7.5) was first given in [25]. With this sketch, the RidgeSketch method with momentum (6.1) becomes CD with momentum which is given by

$$ (7.6) \qquad w^{k+1} = w^k - \gamma_k \frac{\mathbf{A}_i w^k - b_i}{\mathbf{A}_{ii}} e_i + \beta_k (w^k - w^{k-1}) . $$

The CD method with constant momentum $\beta_k \equiv \beta$ is known to converge linearly [25, 27]. But the linear convergence of CD with momentum is always slower than CD without momentum; see Theorem 1 in [27]. Here we present the first convergence rate of CD with momentum that can be faster than CD without momentum. But first, we need the following corollary.

**COROLLARY 7.2.** Let $\epsilon > 0$. If we set the momentum parameters $\gamma_k$ and $\beta_k$ according to (6.12) with $\eta \equiv 0.5$, then the iterates of CD with momentum (7.6) satisfy

$$ (7.7) \qquad t \ge 4 \frac{\text{tr}(\mathbf{A})}{\epsilon} \implies \frac{\mathbb{E}[\|\mathbf{A}w^t - b\|^2]}{\|w^0 - w^*\|^2_{\mathbf{A}}} < \epsilon. $$

Alternatively, if we use no momentum (setting $\gamma_k = 1$ and $\beta_k = 0$), then the iterates (7.6) satisfy

$$ (7.8) \qquad t \ge \frac{\text{tr}(\mathbf{A})}{\lambda_{\min}(\mathbf{A})} \log \left( \frac{\lambda_{\max}(\mathbf{A})}{\epsilon} \right) \implies \frac{\mathbb{E}[\|\mathbf{A}w^t - b\|^2]}{\|w^0 - w^*\|^2_{\mathbf{A}}} < \epsilon. $$

*Proof.* Consider the sampling given by (7.5). From (7.4), since $\mathbf{F} = \mathbf{I}_m$, we have that

$$ (7.9) \qquad \lambda_{\min}(\mathbb{E}[\mathbf{H}_S]) \ge \frac{1}{\text{tr}(\mathbf{A})}. $$

Consequently, using (7.9) together with the complexity bound (6.13) and setting $\mathbf{W} = \mathbf{A}$ we have (7.7).
The linear complexity (7.8) follows from a special case of Theorem 5.1. Indeed, by using (7.9) we have that

$$ \rho = \frac{\lambda_{\min}(\mathbf{A})}{\mathrm{tr}(\mathbf{A})}. $$

The resulting complexity (7.8) follows by standard manipulations of logarithm. □

This linear convergence (7.8) is generally preferred because of the resulting logarithmic dependency of $\epsilon$. But, as we show next, the sublinear complexity given in (7.7) can be tighter when $\epsilon$ is not too small.

**COROLLARY 7.3 (domain of superiority of the sublinear over the linear convergence).** Consider the setting of Corollary 7.2. Let $\hat{\epsilon} := \frac{\epsilon}{\lambda_{\max}(\mathbf{A})}$ be the scaled precision and let $\kappa := \frac{\lambda_{\max}(\mathbf{A})}{\lambda_{\min}(\mathbf{A})}$ denote the condition number. If

$$ (7.10) \qquad \hat{\epsilon}(1 - \hat{\epsilon}) \geq \frac{4}{\kappa}, $$

then the complexity bound of momentum (7.7) is tighter than the bound in (7.8). Furthermore, if $\kappa \geq 16$, the solutions to (7.10) in $\hat{\epsilon}$ are given by

$$ (7.11) \qquad \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{16}{\kappa}} \leq \hat{\epsilon} \leq \frac{1}{2} + \frac{1}{2} \sqrt{1 - \frac{16}{\kappa}}. $$

*Proof.* The complexity bound in (7.7) is tighter than the bound in (7.8) if

$$ \frac{\mathrm{tr}(\mathbf{A})}{\lambda_{\min}(\mathbf{A})} \log \left( \frac{\lambda_{\max}(\mathbf{A})}{\epsilon} \right) \geq 4 \frac{\mathrm{tr}(\mathbf{A})}{\epsilon}. $$

Substituting $\hat{\epsilon} := \frac{\epsilon}{\lambda_{\max}(\mathbf{A})}$ and rearranging the above gives

$$ (7.12) \qquad \hat{\epsilon} \log \left( \frac{1}{\hat{\epsilon}} \right) \geq \frac{4}{\kappa}. $$

To further bound the above we use the standard logarithm bound

$$ \frac{1}{1-x} \log \left( \frac{1}{x} \right) \geq 1 \quad \text{for all } x \in (0, 1), $$

which after manipulations gives

$$ x \log \left( \frac{1}{x} \right) \geq x(1-x). $$

Assuming that $\hat{\epsilon} = \epsilon/\lambda_{\max}(\mathbf{A}) \leq 1$ and using this bound with $x = \hat{\epsilon}$ we have that if (7.10) holds, then

$$ \hat{\epsilon} \log \left( \frac{1}{\hat{\epsilon}} \right) \geq \hat{\epsilon}(1 - \hat{\epsilon}) \geq \frac{4}{\kappa}. $$

Thus (7.12) holds. □

Using Corollary 7.3 we can deduce several regimes where the sublinear momentum bound in (7.7) is tighter than the bound in (7.8). As illustrated in Figure 1, when the condition number is moderate to large or the scaled precision is moderate, the sublinear bound (7.7) is often tighter. On the other hand, when $\hat{\epsilon}$ is very small, then linear rates such as (7.8) are generally preferred.

The regime where sketch-and-projecting methods are interesting is when $\epsilon$ is moderate, and $\mathbf{A}$ has large dimensions. Furthermore, in the large dimensional setting, the condition number of $\mathbf{A}$ can also be very large, thus (7.10) is likely to hold.
**8. RidgeSketch momentum experiments.** Our first experiment explores the efficiency of the momentum version of our RidgeSketch method. We then compare the different sketches described in section 4. Finally, we prove the efficiency of our method on large scale real datasets and show it is competitive with CG and direct solvers.

*Datasets.* In what follows, we test our algorithms on the datasets with different numbers of data samples $n$ and numbers of features $d$: California Housing ($n = 20,640, d = 8$), Boston ($n = 506, d = 13$), RCV1 ($n = 804,414, d = 47,236$) fetched from sklearn<fnref n="5" /> [34] and Year Prediction MSD ($n = 515,345, d = 90$) from the UCI repository.<fnref n="6" /> We converted RCV1 into a regression task by transforming multiclass labels into integers.

**8.1. Experiment 1: Comparison of different momentum settings.** In this section, we compare the sketch-and-project method with momentum for three settings: our new iteration dependent parameters given by (6.4), the constant $\beta$ setting proposed in [27] ($\gamma_k = 1, \beta_k = 0.5$), and no momentum at all ($\gamma_k = 1, \beta_k = 0$). We report iteration plots since the sketch-and-project methods with or without momentum have almost the same iteration cost. Indeed in either case, this cost is dominated by sketching $\mathbf{A}$, computing the residual $r^k$, and solving the sketched system (3.6). We report error areas (first and third quartiles) computed over 10 runs each.

*Increasing momentum.* As suggested by Theorem 6.2, our momentum algorithm, Algorithm 6.1, requires choosing the sequence $\eta_k$, after which $\gamma_k$ and $\beta_k$ are set using (6.6) and (6.4). After running several benchmark tests, we identified the following theoretical rule for setting $\eta_k$:

$$ (8.1) \qquad \eta_k = \begin{cases} 0.995 & \text{if } \beta_k < 0.5, \\ 1 & \text{if } \beta_k \ge 0.5. \end{cases} $$

We call this parameter choice increasing momentum as it allows $\beta_k$ to increase from 0 to 0.5 while the step size $\gamma_k$ decreases from 1 to 0.5, as shown in Figure 3.

Yet, we found in our experiments that this theoretical setting (8.1) closely matches the version without momentum of Algorithm 3.1. We suppose that the gain of the increasing momentum is lost by an excessively rapid drop of the step size to 0. This

<figure><img src="image_5.png" /><figcaption>FIG. 3. Momentum parameters for the theoretical rule: $\eta_k \equiv 0.5$ while $\beta_k < 0.5$, then $\eta_k \equiv 1$.</figcaption></figure>

<footnote n="5">⁵https://scikit-learn.org/stable/modules/generated/sklearn datasets.</footnote>

<footnote n="6">⁶https://archive.ics.uci.edu/ml/datasets/yearpredictionmsd.</footnote>
is why we introduce the following heuristic setting that keeps the step size to 1 and still uses the theoretical setting for momentum when $\eta_k \equiv \eta$ given by (6.12), that is,

$$ (8.2) \qquad \gamma_k \equiv 1 \quad \text{and} \quad \beta_k = 1 - \frac{2 - \eta}{(k+1)(1-\eta) + 1}. $$

We tested it for increasing momentum on the Boston and RCV1 datasets with different sketches and sketch sizes; see Figures 4 and 5. We found that this heuristic setting (8.2) had the best of both worlds, in that in the first iterations, when $\gamma_k \approx 1$ and $\beta_k \approx 0$, it benefits from the fast initial decrease of the no-momentum version. Then, in later iterations, it exploits the fast asymptotic convergence of momentum since $\beta_k \approx 0.5$.

**Regularization test.** Using the heuristic setting (8.2), we tested the impact of using a small, a medium, and a larger regularization parameter $\lambda$ on the performance of momentum; see Figure 4. In this figure, we can see that constant momentum is less effective as $\lambda$ increases, and the no-momentum variant is more effective when $\lambda$ is small. Moreover, we observe the robustness of our heuristic increasing momentum since it performs well for all regularizers.

**Sketch size.** We tested different values of the sketch size, namely $\tau = 10\%$, 50%, and 90% of $m$, and reported the run time to reach a tolerance of $10^{-4}$ for each method. In Figure 5, we observe that constant momentum is very affected by the sketch size and is always the slowest method. For intermediate sketch sizes, like $\tau = m/2$, our increasing momentum competes with no momentum. We also see that $\tau$ should not be set too small or too large. Indeed, larger sketch sizes lead to better estimates of the

<figure><img src="image_6.png" /><figcaption>FIG. 4. Effect of regularizer $\lambda$ on increasing momentum (green), constant momentum (red), and no momentum (blue) (Boston dataset with kernel, $m = 506$, $\tau = \lfloor m/4 \rfloor = 126$, count sketch).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5. Effect of sketch size on increasing momentum (green), constant momentum (red), and no momentum (blue) (California Housing with kernel, $m = 20,640$, $\lambda = 10^{-6}$, Subsample sketch).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison of sketch methods for heuristic increasing momentum on kernel ridge regression problem applied to the Boston dataset ($\lambda = 10^{-6}$, $m = 506$, and $\tau = m/4 = 126$).</figcaption></figure>

initial system (2.3) by (3.1). But if the sketch size is too large, solving the sketched system (3.6) becomes very slow.

*Conclusions*. We highlighted that momentum sketch-and-project is more efficient for small regularizers $\lambda$ as opposed to the vanilla method. Also, we showed that the run time decreases then increases as a function of the sketch size $\tau$. Thus $\tau$ should be set to an intermediate value, e.g., $\tau = m^{2/3}$, so that the cost of solving the sketched system (3.6) is manageable. Finally, the main conclusion of this experiment is the overall robustness (across values of $\lambda$ and $\tau$) and faster convergence of our heuristic increasing momentum setting.

## 8.2. Experiment 2: Comparison of different types of sketches.
In this experiment we compare the performance of different sketching methods presented in section 4 when using our heuristic increasing momentum setting (8.2). In Figure 6, we monitor both the number of iterations and the time taken, since different sketching methods take different amounts of time per iteration. We see in this figure that there is a clear ranking between sketch methods in terms of run time:

1. Subsample is the most efficient on dense data (see also Figure 7(a)).

2. Count and SubCount are the most efficient on sparse data (see Figure 7(b)) and have very similar performance.

3. Gaussian is slow because of the cost of dense matrix-matrix multiplications.

4. Hadamard is extremely slow because of the size of the padded matrix and of the preprocessing time it requires.

*Conclusions*. For dense datasets, the Subsample sketch is the fastest because it only requires slicing operations, which are very well optimized (especially for NumPy arrays). For sparse problems, the count sketch is to be preferred since it *densifies* just enough sketched matrices to extract information out of **A**. We find that computing Gaussian and Hadamard sketch is very time demanding. Furthermore, the cost associated to the padding step in Hadamard sketch is detrimental, especially for large *m*, which often makes it the slowest method.

## 8.3. Experiment 3: Comparison against direct solver and conjugate gradients.
We now compare RidgeSketch with our heuristic increasing momentum setting (8.2) with the two best sketches, Subsample and Count sketch, against a direct solver and CG [19]. The direct solver we used was LAPACK's gesv routine [3] for solving positive definite linear systems. Here we tested our code on

* a kernel ridge regression problem (2.6) on the dataset California Housing ($m = n = 20,640$),
<figure><img src="image_9.png" /><figcaption>FIG. 7. Time comparison of best RidgeSketch methods (heuristic increasing momentum) against direct and CG solvers ($\lambda = 10^{-9}$).</figcaption></figure>

• a large and sparse dataset, RCV1 ($m = d = 47,236$), with only 0.16% of nonzeros.

Figure 7 highlights a clear benefit of iterative methods like CG and sketch-and-project for solving large scale ridge problems. Moreover, this experiment on RCV1 shows that in the large scale sparse setting, count sketch is competitive as compared to CG.

**9. Acceleration.** Recently, it was shown that the convergence rate of sketch-and-project can be improved by using acceleration [45, 15]. See Algorithm 9.1 for our pseudocode of the accelerated sketch-and-project method. In [15] it was shown that, by using specific parameter settings, the accelerated method enjoys a linear convergence with a rate that can be an order of magnitude better than rate given in Theorem 5.1.

Despite this strong theoretical advantage of the accelerated method, it is not clear if this translates into a practical advantage because (1) the additional overhead costs of the method may outweigh the benefits of the improved iteration complexity and (2) the accelerated method relies on knowing beforehand spectral properties of the matrix $A$ that are expensive to compute. Here we show that (1) can be remedied by a careful implementation and (2) is indeed a fundamental issue that prevented us from developing a practical method.

**9.1. From theory to practical implementation of acceleration.**

*Additional overhead: Pseudocode and efficient implementation.* The accelerated version in Algorithm 9.1 has three ($w^k, z^k, v^k$) sequences of iterates. The bottleneck costs of Algorithm 9.1 are the same as the standard sketch-and-project method in Algorithm 3.1, which are the sketching operations on line 12. Indeed, the only additional computations in Algorithm 9.1 as compared to Algorithm 3.1 are lines 14 and 16 which cost $O(m)$. The other additional overhead is how to monitor the residual so as to know when to stop the algorithm. We found that for the residual to be efficiently maintained and updated, we had to monitor three residual vectors $r_z^k := Az^k - b$, $r_w^k := Aw^k - b$, and $r_v := Av^k - b$. These residual vectors can be updated efficiently since

$$r_w^{k+1} = Aw^{k+1} - b = A(z^k - g^k) - b = r_z^k - Ag^k = r_z^k - AS_k\delta^k.$$
Since we have already precomputed $\mathbf{AS}_k$ and $\delta^k$, the additional cost is $O(m\tau)$. Furthermore from lines 14 and 16 we have that

$$r_z^{k+1} = \alpha \mathbf{A} v^k + (1-\alpha) \mathbf{A} w^k - b = \alpha r_v^k + (1-\alpha) r_w^k$$

<caption>**Algorithm 9.1** The sketch-and-project method with acceleration.</caption>

1: **parameters:** distribution over random $m \times \tau$ matrices in $\mathcal{D}$, tolerance $\epsilon > 0$, acceleration parameters $\mu \in [0, 1]$, $\nu \in [1, \frac{1}{\mu}]$
2: Set $w^0 = v^0 = 0 \in \mathbb{R}^m$ ▷ weights initialization
3: Set $r_z^0 = r_v^0 = r_w^0 = \mathbf{A}w^0 - b \in \mathbb{R}^m$ ▷ residual initialization
4: Set $\beta = 1 - \sqrt{\frac{\mu}{\nu}}$ ▷ weights initialization
5: Set $\gamma = \sqrt{\frac{1}{\mu\nu}}$ ▷ residual initialization
6: Set $\alpha = \frac{1}{1+\sqrt{\frac{\mu}{\nu}}}$ ▷ weights initialization
7: $k = 0$
8: **while** $\|r_v^k\|_2 / \|r_v^0\|_2 \le \epsilon$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;Sample $\mathbf{S}_k \sim \mathcal{D}$ i.i.d
10: &nbsp;&nbsp;&nbsp;&nbsp;Compute and store $\mathbf{AS}_k$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\delta^k = \text{least\_norm\_solution} (\mathbf{S}_k^\top \mathbf{AS}_k, \mathbf{S}_k^\top r_z^k)$$ ▷ solve sketched system
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$g^k = \mathbf{S}_k \delta^k$$
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\mathbf{A}g^k = (\mathbf{AS}_k)\delta^k$$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$z^k = \alpha v^k + (1-\alpha)w^k$$ ▷ update the iterates
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$w^{k+1} = z^k - g^k$$
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$v^{k+1} = \beta v^k + (1-\beta)z^k - \gamma g^k$$
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$r_z^k = \alpha r_v^k + (1-\alpha)r_w^k$$ ▷ update the residuals
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$r_w^{k+1} = r_z^k - \mathbf{A}g^k$$
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$r_v^{k+1} = \beta r_v^k + (1-\beta)r_z^k - \gamma \mathbf{A}g^k$$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$k = k + 1$$
21: **end while**
22: **return** $w^t$ ▷ return weights vector

and

$$r_v^{k+1} = \beta \mathbf{A} v^k + (1-\beta) \mathbf{A} z^k - b - \gamma \mathbf{A} g^k = \beta r_v^k + (1-\beta) r_z^k - \gamma \mathbf{A} \mathbf{S}_k \delta^k .$$

Thus the residuals $r_v^k$ and $r_z^k$ can be updated at an additional $O(m)$ cost to perform the above vector additions and scalar multiplications.

*Setting the acceleration parameters with spectral properties.* The main issue with the accelerated version is that it introduces two new hyperparameters $\mu$ and $\nu$ which have to be estimated. In theory [15], by setting these two parameters according to

$$(9.1) \quad \mu := \inf_{x \in \text{Range}(\mathbf{A}^\top)} \frac{\langle \mathbb{E}[\mathbf{Z}|x,x] \rangle}{\langle x,x \rangle} \quad \text{and} \quad \nu := \sup_{x \in \text{Range}(\mathbf{A}^\top)} \frac{\langle \mathbb{E}[\mathbf{Z}\mathbb{E}[\mathbf{Z}]^\top \mathbf{Z}]|x,x \rangle}{\langle \mathbb{E}[\mathbf{Z}]|x,x \rangle},$$

where

$$(9.2) \qquad \mathbf{Z} := \mathbf{A}^\top \mathbf{S}^\top (\mathbf{S}^\top \mathbf{A} \mathbf{S})^\dagger \mathbf{S}^\top \mathbf{A},$$

we can guarantee an accelerated rate of convergence. The issue is that the theory in [15] requires that these parameters be set exactly using (9.1) and computing (9.1)
is more costly then solving the original linear system! So this leads us to the following practical question.

Is there a rule of thumb setting for the acceleration parameters $\mu$ and $\nu$ such that the accelerated sketch-and-project method is consistently better than the standard sketch-and-project method?

In [45] the authors propose some settings for $\nu$ and $\mu$ when the sketch size is large. But there is currently no practical rule for setting these parameters in general. In theory, we know that

$$0 \le \mu \le \frac{1}{\nu} \le 1,$$

as proven in Lemma 2 in [15]. Furthermore the extreme case where $\mu = \nu = 1$ corresponds to the standard sketch-and-project method, as can be seen by induction on Algorithm 9.1 since $z^k = v^k = w^k$ for all iterations, and $w^k$'s are thus equivalent to the $w^k$'s in Algorithm 3.1. We now look at some other extreme cases to better understand these parameters.

*Single row sampling.* For this special case of Subsample sketchs with $\tau = 1$, that is, $S = e_i$, where we recall that $(e_i)_{1 \le i \le m}$ are the canonical basis vectors of $\mathbb{R}^m$, with probability $\frac{1}{m}$ we know that

$$(9.3) \qquad \mu = \frac{\lambda_{\min}(\mathbf{A})}{\mathrm{tr}(\mathbf{A})} \quad \text{and} \quad \nu = \frac{\mathrm{tr}(\mathbf{A})}{\min_{i=1,\dots,m} \mathbf{A}_{ii}}.$$

Consequently, if the eigenvalues of $\mathbf{A}$ are concentrated with $\lambda_{\min}(\mathbf{A})$ close to $\lambda_{\max}(\mathbf{A})$, then we have that $\mu \approx \frac{1}{m}$ and $\nu \approx m$. Alternatively, if the eigenvalues of $\mathbf{A}$ are far apart, then it may be that $\mu \approx 0$ and $\nu \approx \infty$.

*No sketching.* When $S = I$ then $Z = A$ since $A$ is invertible. Consequently

$$(9.4) \qquad \mu = \inf_{x \in \mathrm{Range}(\mathbf{A}^{\top})} \frac{\langle \mathbf{A}x, x \rangle}{\langle x, x \rangle} = \lambda_{\min}(\mathbf{A}), \qquad \nu = \sup_{x \in \mathrm{Range}(\mathbf{A}^{\top})} \frac{\langle \mathbf{E}[\mathbf{A}]x, x \rangle}{\langle \mathbf{E}[\mathbf{A}]x, x \rangle} = 1.$$

In either of these two extremes, we need the smallest eigenvalue of $\mathbf{A}$ to set $\mu$ and $\nu$, which is a prohibitive cost. In subsection 9.2 we show that finding a setting for $\mu$ and $\nu$ that outperforms the standard sketch-and-project method is difficult, akin to finding a needle in a haystack.

**9.2. Experiments setting the acceleration parameters.** Here we would like to verify if there exists a default setting for the acceleration parameters $\mu$ and $\nu$ that results in consistently faster execution than the nonaccelerated version. In Figure 8(a), we show the results of an extensive grid search for trying to identify a suitable $\mu$ and $\nu$ setting. This figure shows the time taken to reach a $\epsilon = 10^{-4}$ solution of the relative residual for different pairs of $\mu$ and $\nu$ such that

$$0 \le \mu \le \frac{1}{\nu} \le 1.$$

The problem we considered here is kernel ridge regression, with an RBF kernel with $\sigma = 0.5$. The data $\mathbf{X}$ is a random $n \times d$ sparse CSC matrix with density 0.25, and the regularizer is set to $\lambda = 1/n$. From Figure 8(a), there is no clear pair of parameters $(\mu, \nu)$ leading to an improvement in convergence. Even if a finer grid search might allow us to find optimal parameters, the gain in convergence is so marginal
<figure><img src="image_10.png" /><figcaption>FIG. 8. Grid time plot (left) and single run (right) to solve an RBF kernel regression ($\sigma = 0.5$, $\lambda = 1/n$) for random sparse CSC matrix $\mathbf{X} \in \mathbb{R}^{n \times d}$ ($n = m = 2000$, $d = 500$) with density 0.25. We used a Subsample sketch method and the sketch size is set to $\tau = [m^{2/3}] = 158$.</figcaption></figure>

that it makes acceleration impractical compared to the version without acceleration; see Figure 8(b).

We conclude that there is yet no known way to set these acceleration parameters in practice: the theory might be too loose to set them and looking empirically on a grid search for optimal points is too cumbersome.

## Appendix A. Auxiliary lemmas.

LEMMA A.1. The linear systems (2.1) and (2.2) have the same solution.

*Proof.* Note that

$$ (\mathbf{X}^{\top}\mathbf{X} + \lambda\mathbf{I})\mathbf{X}^{\top} = \mathbf{X}^{\top}(\mathbf{X}\mathbf{X}^{\top} + \lambda\mathbf{I}) . $$

Left multiplying the above by the inverse of $(\mathbf{X}^{\top}\mathbf{X} + \lambda\mathbf{I})$ and right multiplying by the inverse of $(\mathbf{X}\mathbf{X}^{\top} + \lambda\mathbf{I})$ gives

$$ \mathbf{X}^{\top}(\mathbf{X}\mathbf{X}^{\top} + \lambda\mathbf{I})^{-1} = (\mathbf{X}^{\top}\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^{\top} . $$

Finally right multiplying by $y$, on the left we have the solution to (2.2) and on the right the solution to (2.1). □

## Appendix B. RidgeSketch package.
Our RidgeSketch Python package is designed to be easily augmented by new contributions. Users are encouraged to add new sketches, new parameter settings (e.g., for momentum), and new datasets (see datasets/data_loaders.py). They can also easily compare methods using the command

```
$ python benchmarks.py [options] [name of config]
```

RidgeSketch comes with two tutorial Jupyter Notebooks: the first one explains how to fit a RidgeSketch model, and the second one indicates how to add new sketches and benchmarks. Next we present some snippets of code.

### B.1. Solving the ridge regression problem.
In Code 1, we provide an example creating a ridge regression model and solving the fitting problem with our Subsample sketch solver.
```python
from datasets.data_loader import BostonDataset
from ridge_sketch import RidgeSketch

dataset = BostonDataset()
X, y = dataset.load_X_y()  # data loader

model = RidgeSketch(
    alpha=1.0,  # regularizer
    algo_mode="mom",  # momentum version
    solver="subsample",  # sketch, cg, direct or sklearn solver
    sketch_size=10)
model.fit(X, y)  # solve the fitting problem
```

<caption>CODE 1. Example of how to build a RidgeSketch model and fit it to a dataset.</caption>

## B.2. Adding sketching methods.
We built our package so that new sketching methods can easily be added by contributors. To do this, one should open the `sketching.py` file and create a new sketching subclass that inherits from the Sketch class. In Code 2 we provide the example of the Subsample sketching method.

```python
import numpy as np

class SubsampleSketch(Sketch):
    def __init__(self, A, b, sketch_size):
        """Initializes the Sketch object and defines its
        attributes"""
        super().__init__(A, b, sketch_size)
        self.sample_indexes = None

    def set sketch(self):
        """Generates subsampling indices representing the
        sketching matrix"""
        self.sample_indexes = np.random.choice(
            self.m, size self sketch_size, replace=False). tolist()

    def sketch(self, r):
        """Generates sketched system"""
        self.set sketch()
        SA = self.A. get_rows(self.sample_indexes)
        SAS = SA[:, self.sample_indexes]
        rs = r[ self sample_indexes]
        return SA, SAS, rs

    def update_iterate(self, w, lmbda, step_size=1.0):
        """Updates the weights by solving the sketched system"""
        w[ self sample_indexes] -= step_size * lmbda
        return w
```

<caption>CODE 2. Implementation of the Subsample sketch subclass.</caption>
[1] N. AILON AND B. CHAZELLE, *The fast Johnson-Lindenstrauss transform and approximate nearest neighbors*, SIAM J. Comput., 39 (2009), pp. 302–322.

[2] N. AILON AND E. LIBERTY, *Fast dimension reduction using Rademacher series on dual BCH codes*, Discrete Comput. Geom., 42 (2009), pp. 615–630.

[3] E. ANDERSON, Z. BAI, C. BISCHOF, S. BLACKFORD, J. D. J. DONGARRA, J. D. CROZ, A. GREENBAUM, S. HAMMARLING, A. MCKENNEY, AND D. SORENSEN, *LAPACK Users' Guide*, 3rd ed., SIAM, Philadelphia, 1999.

[4] H. AVRON, P. MAYMOUNKOV, AND S. TOLEDO, *Blendenpik: Supercharging LAPACK's Least-Squares Solver*, SIAM J. Sci. Comput., 32 (2010), pp. 1217–1236, http://dx.doi.org/10.1137/090767911.

[5] Z.-Z. BAI AND W.-T. WU, *On greedy randomized Kaczmarz method for solving large sparse linear systems*, SIAM J. Sci. Comput., 40 (2018), pp. A592–A606.

[6] C. BOUTSIDIS AND A. GITTENS, *Improved matrix algorithms via the subsampled randomized Hadamard transform*, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 1301–1340.

[7] M. CHARIKAR, K. CHEN, AND M. FARACH-COLTON, *Finding frequent items in data streams*, in Proceedings of the 29th International Colloquium on Automata, Languages and Programming, Springer-Verlag, London, 2002, pp. 693–703.

[8] K. L. CLARKSON AND D. P. WOODRUFF, *Low-rank approximation and regression in input sparsity time*, J. ACM, 63 (2017), pp. 1–45.

[9] G. CORMODE, *Count-min sketch*, in Encyclopedia of Database Systems, Springer, New York, 2009, pp. 511–516, https://doi.org/10.1007/978-0-387-39940-9_87.

[10] G. CORMODE AND S. MUTHUKRISHNAN, *An improved data stream summary: The count-min sketch and its applications*, J. Algorithms, 55 (2005), pp. 58–75.

[11] J. A. DE LOERA, J. HADDOCK, AND D. NEEDELL, *A sampling Kaczmarz-Motzkin algorithm for linear feasibility*, SIAM J. Sci. Comput., 39 (2017), pp. S66–S87.

[12] K. DU AND H. GAO, *A new theoretical estimate for the convergence rate of the maximal weighted residual Kaczmarz algorithm*, Numer. Math. Theory Methods Appl., 12 (2019), pp. 627–639.

[13] B. J. FINO AND V. R. ALGAZI, *Unified matrix treatment of the fast Walsh-Hadamard transform*, IEEE Trans. Computers, 25 (1976), pp. 1142–1146.

[14] G. H. GOLUB, P. C. HANSEN, AND D. P. O'LEARY, *Tikhonov regularization and total least squares*, SIAM J. Matrix Anal. Appl., 21 (1999), pp. 185–194.

[15] R. GOWER, F. HANZELY, P. RICHTARIK, AND S. STICH, *Accelerated stochastic matrix inversion: General theory and speeding up BFGS rules for faster second-order optimization*, in Advances in Neural Information Processing Systems 31, 2018, pp. 2292–2300.

[16] R. M. GOWER AND P. RICHTARIK, *Randomized iterative methods for linear systems*, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 1660–1690.

[17] R. M. GOWER AND P. RICHTARIK, *Stochastic Dual Ascent for Solving Linear Systems*, arXiv:1512.06890, 2015.

[18] J. HADDOCK AND A. MA, *Greed works: An improved analysis of sampling Kaczmarz-Motzkin*, SIAM J. Math. Data Sci., 3 (2021), pp. 342–368.

[19] M. R. HESTENES AND E. STIEFEL, *Methods of conjugate gradients for solving linear systems*, J. Res. National Bureau of Standards, 49 (1952).

[20] T. S. JAYRAM AND D. P. WOODRUFF, *Optimal bounds for Johnson-Lindenstrauss transforms and streaming problems with subconstant error*, ACM Trans. Algorithms, 9 (2013), pp. 26:1–26:17.

[21] W. JOHNSON AND J. LINDENSTRAUSS, *Extensions of Lipschitz mappings into a Hilbert space*, in Conference in Modern Analysis and Probability (New Haven, CT, 1982), Contemp. Math. 26, AMS, Providence, RI, 1984, pp. 189–206.

[22] M. S. KACZMARZ, *Angenäherte auflösung von systemen linearer gleichungen*, Bulletin International de l'Académie Polonaise des Sciences et des Lettres. Classe des Sciences Mathématiques et Naturelles. Série A, Sciences Mathématiques, 35 (1937), pp. 355–357.

[23] T. N. KRISHNAMURTI, *Numerical weather prediction*, Annu. Rev. Fluid Mech., 27 (1995), pp. 195–224.

[24] Y. T. LEE AND A. SIDFORD, *Efficient accelerated coordinate descent methods and faster algorithms for solving linear systems*, in Proceedings of the Annual IEEE Symposium on Foundations of Computer Science, 2013, pp. 147–156, https://doi.org/10.1109/FOCS.2013.24.

[25] D. LEVENTHAL AND A. S. LEWIS, *Randomized methods for linear constraints: Convergence rates and conditioning*, Math. Oper. Res., 35 (2010), pp. 641–654, https://doi.org/10.1287/moor.1100.0456.
[26] J. LIU AND S. J. WRIGHT, *An accelerated randomized Kaczmarz algorithm*, Math. Comp., 85 (2016), pp. 153–178.

[27] N. LOIZOU AND P. RICHTARIK, *Momentum and stochastic momentum for stochastic gradient, Newton, proximal point and subspace descent methods*, Comput. Optim. Appl., 77 (2020), pp. 653–710.

[28] A. MA, D. NEEDELL, AND A. RAMDAS, *Convergence properties of the randomized extended Gauss-Seidel and Kaczmarz methods*, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 1590–1604.

[29] M. W. MAHONEY, *Randomized algorithms for matrices and data*, Found. Trends Mach. Learn., 3 (2011), pp. 123–224.

[30] M. S. MORSHED, S. AHMAD, AND M. NOOR-E-ALAM, *Stochastic Steepest Descent Methods for Linear Systems: Greedy Sampling & Momentum*, arXiv:2012.13087, 2020.

[31] M. S. MORSHED, M. S. ISLAM, AND M. NOOR-E-ALAM, *Accelerated sampling Kaczmarz Motzkin algorithm for the linear feasibility problem*, J. Global Optim., 77 (2020), pp. 361–382.

[32] M. S. MORSHED AND M. NOOR-E-ALAM, *Sketch & Project Methods for Linear Feasibility Problems: Greedy Sampling & Momentum*, arXiv:2012.02913, 2020.

[33] I. NECOARA, *Faster randomized block Kaczmarz algorithms*, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 1425–1452.

[34] F. PEDREGOSA, G. VAROQUAUX, A. GRAMFORT, V. MICHEL, B. THIRION, O. GRISEL, M. BLONDEL, P. PRETTENHOFER, R. WEISS, V. DUBOURG, J. VANDERPLAS, A. PASSOS, D. COURNAPEAU, M. BRUCHER, M. PERROT, AND E. DUCHESNAY, *SCIKIT-learn: Machine learning in Python*, J. Mach. Learn. Res., 12 (2011), pp. 2825–2830.

[35] S. PETRA AND C. POPA, *Single projection Kaczmarz extended algorithms*, Numer. Algorithms, 73 (2016), pp. 791–806.

[36] M. PILANCI AND M. WAINWRIGHT, *Randomized sketches of convex programs with sharp guarantees*, Proc. IEEE Trans. Inform. Theory, 61 (2015), pp. 5096–5115.

[37] M. PILANCI AND M. J. WAINWRIGHT, *Iterative Hessian sketch: Fast and accurate solution approximation for constrained least-squares*, J. Mach. Learn. Res., 17 (2016), pp. 1–33.

[38] B. T. POLYAK, *Some methods of speeding up the convergence of iteration methods*, USSR Comput. Math. Math. Phys., 4 (1964), pp. 1–17.

[39] P. RICHTÁRIK AND M. TAKÁČ, *Stochastic reformulations of linear systems: Algorithms and convergence theory*, SIAM J. Matrix Anal. Appl., 41 (2020), pp. 487–524.

[40] G. SAUNDERS, A. GAMMERMAN, AND V. VOVK, *Ridge regression learning algorithm in dual variables*, in Proceedings of the 15th International Conference on Machine Learning, Morgan Kaufmann, San Francisco, CA, 1998, pp. 515–521.

[41] O. SEBBOUH, R. M. GOWER, AND A. DEFAZIO, *On the Convergence of the Stochastic Heavy Ball Method*, arXiv:2006.07867, 2020.

[42] S. SHALEV-SHWARTZ AND S. BEN-DAVID, *Understanding Machine Learning: From Theory to Algorithms*, Cambridge University Press, Cambridge, UK, 2014.

[43] T. STROHMER AND R. VERSHYNIN, *A randomized Kaczmarz algorithm with exponential convergence*, J. Fourier Anal. Appl., 15 (2009), pp. 262–278.

[44] J. A. TROPP, *Improved analysis of the subsampled randomized Hadamard transform*, Adv. Adapt. Data Anal., 3 (2011), pp. 115–126.

[45] S. TU, S. VENKATARAMAN, A. C. WILSON, A. GITTENS, M. I. JORDAN, AND B. RECHT, *Breaking locality accelerates block Gauss-Seidel*, in Proceedings of the 34th International Conference on Machine Learning, Sydney, Australia, 2017, pp. 3482–3491.

[46] S. VAN DER WALT, S. C. COLBERT, AND G. VAROQUAUX, *The Numpy array: A structure for efficient numerical computation*, Comput. Sci. Eng., 13 (2011), pp. 22–30.

[47] P. VIRTANEN, R. GOMMERS, T. E. OLIPHANT, M. HABERLAND, T. REDDY, D. COURNAPEAU, E. BUROVSKI, P. PETERSON, W. WECKESSER, J. BRIGHT, ET AL., *SCIPY 1.0: Fundamental algorithms for scientific computing in Python*, Nature Methods, 17 (2020), pp. 261–272.

[48] V. VOVK, *Kernel ridge regression*, in Empirical Inference, Springer, New York, 2013, pp. 105–116.

[49] S. WANG, A. GITTENS, AND M. W. MAHONEY, *Sketched ridge regression: Optimization perspective, statistical perspective, and model averaging*, J. Mach. Learn. Res., 18 (2017), pp. 218:1–218:50.

[50] F. WOOLFE, E. LIBERTY, V. ROKHLIN, AND M. TYGERT, *A fast randomized algorithm for the approximation of matrices*, Appl. Comput. Harmon. Anal., 25 (2008), pp. 335–366.

[51] S. J. WRIGHT, *Coordinate descent algorithms*, Math. Program., 151 (2015), pp. 3–34.

[52] R. YUAN, A. LAZARIC, AND R. M. GOWER, *Sketched Newton-Raphson*, arXiv:2006.12120, 2020.