{'fraction_non_alphanumeric': 0.05917513448894202, 'fraction_numerical': 0.01687224908982231, 'mean_word_length': 3.5819743557823975, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 38, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective Dimension Adaptive Sketching Methods for Faster Regularized Least-Squares Optimization

We propose a new randomized algorithm for solving L2-regularized least-squares problems based on sketching. We consider two of the most popular random embeddings, namely, Gaussian embeddings and the Subsampled Randomized Hadamard Transform (SRHT). While current randomized solvers for least-squares optimization prescribe an embedding dimension at least greater than the data dimension, we show that the embedding dimension can be reduced to the effective dimension of the optimization problem, and still preserve high-probability convergence guarantees. In this regard, we derive sharp matrix deviation inequalities over ellipsoids for both Gaussian and SRHT embeddings. Specifically, we improve on the constant of a classical Gaussian concentration bound whereas, for SRHT embeddings, our deviation inequality involves a novel technical approach. Leveraging these bounds, we are able to design a practical and adaptive algorithm which does not require to know the effective dimension beforehand. Our method starts with an initial embedding dimension equal to 1 and, over iterations, increases the embedding dimension up to the effective one. Finally, we prove that our algorithm improves the state-of-the-art computational complexity for solving regularized least-squares problems. Further, we show numerically that it outperforms standard least-squares solvers such as the conjugate gradient method and its pre-conditioned version on several standard machine learning datasets.


Introduction
We study the performance of a randomized method, namely, the Hessian sketch [22], in the context of regularized least-squares problems, where A ∈ R n×d is a (dense) data matrix and b ∈ R n is a vector of observations. For clarity purposes and without loss of generality (by considering instead the dual problem of (1)), we make the assumption that the problem is over-determined, i.e., n d and that rank(A) = d.
have time complexity which usually scales proportionally to the condition number κ of A (or √ κ with acceleration) in order to find a solution x with acceptable accuracy. This also becomes prohibitively large when κ 1. Besides the computational complexity, the number of iterations of an iterative solver is also a relevant performance metric in the large-scale setting, as distributed computation may be necessary at each iteration. In this regard, randomized preconditioning methods [24,2,19] involve using a random matrix S ∈ R m×n with m n to project the data A, and then improve the condition number of A based on a spectral decomposition of SA. On the other hand, the iterative Hessian sketch (IHS) introduced by [22] and considered in [20,14,15,21] addresses the conditioning issue differently. Given x 0 , x 1 ∈ R d , it uses a pre-conditioned Heavy-ball update with step size µ and momentum parameter β, given by where the Hessian A A of f (x) is approximated by H S : = A S S A and S is a sketching matrix. We refer to the update (2) as the Polyak-IHS method, and, in the absence of acceleration (β = 0), we call it the Gradient-IHS method. In contrast to preconditioning methods [24,2,19], the IHS does not need to pay the full cost O(md min{m, d}) for decomposing the matrix SA. Although solving exactly the linear system H S · z = ∇f (x t ) also takes time O(md min{m, d}), approximate solving (using for instance CG) is also efficient and faster in practice [21,20].
The choice of the sketching matrix S is critical for statistical and computational performances. A classical sketch is a matrix S with independent and identically distributed (i.i.d.) Gaussian entries N (0, m −1 ) for which forming SA requires in general O(mnd) basic operations (using classical matrix multiplication). On the other hand, it has been observed [16,17,8] and also formally proved [7,15] in several contexts that random projections with i.i.d. entries degrade the performance of the approximate solution compared to orthogonal projections. In this regard, the SRHT [1] is an orthogonal embedding for which the sketch SA can be formed in O(nd log m) time, and this is much faster than Gaussian projections. Consequently, along with the statistical benefits of orthogonal projections, this suggests to use the SRHT as a reference point for comparing sketching algorithms.
In the context of unregularized least-squares problems (ν = 0), [14] showed that the error A(x t − x * ) 2 2 of the Polyak-IHS method is smaller than (d/m) t for both Gaussian and SRHT matrices provided that m ≈ d. More recently, it has been shown in [15] that the scaling (d/m) t is exact for Gaussian embeddings in the asymptotic regime where we let the relevant dimensions go to infinity, whereas the exact scaling for the SRHT is slightly smaller than (d/m) t .
In the regularizated case (ν > 0), more relevant than the matrix rank is the effective dimension d e ≈ d i=1 σ 2 i σ 2 i +ν 2 (where the σ i 's are the singular values of A) which always satisfies d e d, and it is significantly smaller than d when the σ i 's have a fast decay. It has been shown in [21] that one can pick m ≈ d e and achieve the error rate (d e /m) t by using the well-structured approximate Hessian t as δ t : = 1 2 A(x t − x * ) 2 , and its convergence rate as lim sup t→∞ (δ t /δ 0 ) . Finally, we denote by U ∈ R n×d a matrix of left singular vectors of A, A : = A ν · I d and U ∈ R (n+d)×d a matrix of left singular vectors of A. Critical to our convergence analysis is the matrix C S : = D(U S SU − I d )D + I d .

Overview of our contributions
Our main contribution is to propose an algorithm based on the Polyak-IHS update (2) that does not require the knowledge of d e , and is still able to achieve the error rate O ((d e /m) t ). Our algorithm starts with an arbitrary m (e.g, m = 1) and, at each iteration, it uses an improvement criterion to decide if it should increase m or not. We prove that the adaptive sketch size satisfies at each iteration m d e and that our algorithm improves on the state-of-the-art computational complexity for solving regularized least-squares problems.
Our algorithmic parameters and improvement criterion depend on the extreme eigenvalues of C S , and it is then critical for optimal performance to have a sharp estimation of these. For Gaussian embeddings, we improve on the constant of the best known concentration bound [13]. Our constant is sharp in a worst-case sense, and our analysis is based on a recent extension [25] of Gordon's min-max theorem [9]. In the SRHT case, although a similar concentration bound was obtained by [6] (see Theorem 1), we provide a novel technical approach that generalizes the analysis techniques from the work of [26].
We evaluate numerically our adaptive algorithm on several standard datasets. We consider two settings: (i) the regularization parameter ν is fixed; (ii) one aims to compute the several solutions along a regularization path. The latter setting is more relevant to many practical applications [29,12] where estimating a proper regularization parameter is essential. In both cases, we show that our method is faster than the standard conjugate gradient method and one of the state-of-the-art randomized preconditioning methods [24].
Finally, we consider the underdetermined case d n. By taking the dual problem of (1) which is itself an overdetermined least-squares problems, we show that our adaptive algorithm and theoretical guarantees easily extend to this setting. We defer the presentation of these results to Appendix B.

Other related work
Other versions of the IHS have been proposed in the literature, especially in the context of unregularized least-squares. A fundamentally different version uses the same update (2) but with refreshed sketching matrices, i.e., a new matrix S is sampled at each iteration and independently of the previous ones, and the approximate Hessian H S is re-computed. Surprisingly, refreshing embeddings does not improve on using a fixed embedding: it results in the same convergence rate in the Gaussian case [14,15] and in a slower convergence rate in the SRHT case [15].
In contrast to randomized preconditioning methods, classical sketching algorithms project both A and b, and then computes x : = argmin x In the closely related work [3], the authors showed that for m ≈ d e /ε, the estimate x satisfies f ( x) (1 + ε)f (x * ). This can result in large m for even medium accuracy, whereas our method yields an ε-approximate solution with m ≈ d e under the mild requirement that the number of iterations T satisfies T ≈ log(1/ε).

Preliminaries
We start with providing deterministic convergence guarantees for the Polyak-and Gradient-IHS methods, and we relate the convergence rates to the extreme eigenvalues of the matrix C S . Similar guarantees were established in [21].
Let S ∈ R m×d be any sketching matrix with arbitrary sketch size m, and denote by γ 1 (resp. γ d ) the largest (resp. smallest) eigenvalue of C S . Since the matrix D U S SU D is positive semi-definite and D 2 < 1, it holds that C S is positive definite. Given two real numbers Λ > λ > 0, we define the S-measurable event E S : = {λ γ d γ 1 Λ}. The proofs of the two next results are essentially based on standard analyses of gradient methods [23], and they are deferred to Appendix C.1. Theorem 1. Consider the step size µ = 2/( 1 λ + 1 Λ ), where κ = Λ/λ. Then, conditional on E S , the Gradient-IHS method satisfies at each iteration Theorem 2. Consider the step size µ = 4/( 1 where κ = Λ/λ. Then, conditional on E S , the Polyak-IHS satisfies lim sup The above rates c gd (λ, Λ) and c p (λ, Λ) will play a critical role in the design of our adaptive method. For the Gradient-IHS method, it should be noted that we are able to monitor the improvement ratio between two consecutive iterates. However, for the Polyak-IHS method, we only obtain an asymptotic guarantee as t → +∞. This standard result regarding the Heavy-ball method [23] essentially follows from the fact that the iterates obey a non-symmetric linear dynamical system so that the spectral and operator norms of the linear dynamics coincide asymptotically, according to Gelfand's formula (we refer to the proof in Appendix C.1 for a formal explanation).

Sharp convergence rates for Gaussian and SRHT embeddings
According to Theorems 1 and 2, we need sharp estimates of the extreme eigenvalues of C S in order to pick optimal parameters for the Polyak-and Gradient-IHS methods, and this will be necessary to design our adaptive method.

The Gaussian case
We provide a concentration bound on the edge eigenvalues γ 1 and γ d of the matrix C S in terms of the aspect ratio de m . Our analysis is based on a generalized Gordon's Gaussian comparison theorem [9,25] and it provides sharper constants than existing results. We defer the proof to Appendix E.1.
Theorem 3. Let S ∈ R m×n be a Gaussian embedding and ρ ∈ (0, 1). If m de ρ then for any If we let d e , m → +∞ while keeping ρ : = de m fixed and take η = o(1/ √ m), then our lower and upper bounds converge to the respective limits 1− D 2 2 + D 2 2 (1− √ ρ) 2 and 1− D 2 2 + D 2 2 (1+ √ ρ) 2 . When D = D 2 · I d , these limits are exact as they correspond to the edges of the support of the Marchenko-Pastur distribution [18], so that our bounds are tight in a worst-case sense.

The SRHT case
For the SRHT, our concentration bound is given in terms of the aspect ratio C(n, d e ) · de log(de) . We defer it to Appendix E.2. We note that similar concentration bounds were obtained by the authors of [6] in the context of approximate matrix multiplication, using different analysis techniques.
Theorem 4. Let ρ ∈ (0, 1) and m C(n, d e ) · de log(de) ρ . Then it holds with probability at least As already mentioned in the previous section, the operator norm D 2 might be unknown in practice, and so would be the lower and upper bounds λ ρ and Λ ρ . However, if ρ is known and m C(n, d e ) · de log(de) ρ , one can use instead the following bounds.

An adaptive method free of the knowledge of the effective dimension
One usually does not have access to the effective dimension d e , and accurate estimates of this quantity can be efficiently obtained only in some restrictive settings (e.g., d e very small [3]).
We now propose a novel adaptive method with time-varying sketch size. Our algorithm does not require the knowledge of d e , but still achieves a fast rate of convergence while keeping m d e .
Suppose first that we are able to compute the per-iteration improvement ratio C t : = δt+1 δt . Starting from an arbitrary initial sketch size m (say m = 1) and fixing an improvement ratio threshold C, our adaptive algorithm based, for instance, on the Gradient-IHS update and Gaussian embeddings would proceed at each iteration as follows. If C t C then it accepts the update x t+1 . Otherwise it rejects the update, increases the sketch size by a constant factor (say m ← 2m) and re-computes the sketched matrix S · A. Since only updates with sufficient improvement are accepted, this method achieves a convergence rate at least as good as the chosen threshold C. Furthermore there is only a finite number of steps K where the updates get rejected. Indeed, according to Theorems 1 and 3, as soon as the sketch size becomes larger than Ω(d e /C) then all the updates are accepted with high probability. Thus, the number K of rejected steps satisfies K log(d e /C)/ log(2).
However, the main issue with the previous method is that we are not able to compute the ratio C t in practice since this requires the knowledge of Ax * . We alleviate this difficulty by considering a quantity closely related to C t , as follows. First, we introduce the error vector e t : = U A(x t − x * ), its re-scaled version E t : = C − 1 2 S e t and the re-scaled error r t : = 1 2 E t 2 . We need the following technical result whose proof is deferred to Appendix D.1.

Lemma 1.
For any sequence of iterates {x t }, it holds that δ t = 1 2 e t 2 and r t = 1 Hence, the re-scaled error r t is equal to 1 2 g t H −1 S g t which can be easily computed 1 . Consequently, we can efficiently monitor the ratio c t : = r t+1 /r t as opposed to C t in order to adapt the sketch size. Provided that c t and C t are close enough, this would yield the desired performance. We describe our proposed method in Algorithm 1.
convergence rates 0 < c p < c gd < 1, GD step size µ gd , Polyak step size µ p 0, momentum parameter β p 0 and Gaussian or SRHT embeddings. 1 Sample S ∈ R m×n and compute S A = SA. 2 Compute and cache Compute the gradients Compute the improvement ratios c gd = g + gd g + gd g t gt and c p = if c gd c gd then 11 Set x t+1 = x + gd , g t+1 = g + gd and g t+1 = g + gd .
14 Sample S ∈ R m×n and compute S A = S · A. 15 Set g t : = H −1 S g t and return to Step 4. At each iteration, Algorithm 1 uses the improvement ratio between Newton decrements instead of the ratio C t . According to Theorems 1 and 2, we can monitor exactly the improvement between two successive iterates of the Gradient-IHS method as opposed to the Polyak-IHS method. Algorithm 1 computes both updates in order to benefit either from the acceleration of the latter or from the stronger convergence guarantees of the former. If both updates do not make enough progress then the sketch size is increased.
Algorithm 1 requires to specify the target convergence rates c gd and c p . It holds then that our method converges at least as fast as c gd . Importantly, the time-varying sketch size m satisfies m = O(d e /c gd ) with Gaussian embeddings and m = O(d e log(d e )/c gd ) with the SRHT. We formalize these guarantees in Theorems 5 and 6 whose proofs are deferred to Appendices C.2 and C.3. Theorem 5. Let ρ, η ∈ (0, 1) be user's choice parameters and consider Gaussian (resp. SRHT) embeddings. Let λ ρ , Λ ρ , c gd (ρ) and c p (ρ) be as given in Definition 3.1 (resp. Definition 3.2), and choose the target convergence rates c gd = c gd (ρ) and c p = c p (ρ). Let µ gd , µ p and β p be the corresponding parameters as given in Theorems 1 and 2. Fix an arbitrary initial sketch size m initial 1. Denote by K ≡ K t the number of rejected updates at iteration t. Then, for Gaussian embeddings, with a(ρ, η) : = it holds with probability at least 1 − 16e −η 2 de/2 that at any iteration, m 2 a(ρ, η) de ρ and K log a(ρ, η) de minitialρ / log(2) + 1. Further, for the SRHT, with a(ρ) : = √ ρ , it holds with probability at least 1 − 9 de that at any iteration, m 2 a(ρ) C(n, d e ) · de log(de) ρ and K log a(ρ) C(n, d e ) de log(de) minitialρ / log(2) + 1.
Theorem 6. Consider the same conditions of Theorem 5. For Gaussian embeddings, we have for any t 1 that δt gd with probability at least 1 − 9/d e .

Time and space complexity
We consider here the SRHT for which computing SA is faster than Gaussian projections. We have the following complexity result, whose proof is deferred to Appendix C.4. Theorem 7. Consider the same conditions of Theorem 5 with the SRHT. Let ε ∈ (0, 1) be a given precision such that ε min{ ν 2 Then it holds with probability at least 1 − 9 de that the number of iterations to reach a solution . Thus the total time complexity C ε of Algorithm 1 verifies The time complexity C ε is decomposed into three terms. Sketching the data matrix takes O(nd log(d e /ρ)) time. The cost O( + nd log(1/ε) log(1/ρ) ). Thus, even with the extra factor log(d e /ρ) due to the rejected steps, our adaptive method improves on the sketching plus factor costs especially when the effective dimension d e is much smaller than the data dimension d and thus on the total complexity.
Regarding space complexity, our method requires O(d · d e log d e /ρ) space to store the sketched matrix SA whereas the other preconditioning methods needs O(d 2 log d/ρ), and this is a significant improvement when d e is much smaller than d.

Numerical experiments
We carry out numerical simulations of Algorithm 1 and we compare it to standard least-squares solvers, that is, the CG method and the randomized preconditioned CG (pCG) [24]. Numerical simulations were carried out a 512Gb desktop station and implemented in Python using its standard numerical linear algebra modules (NumPy and SciPy).
We consider two evaluation criteria: (i) the cumulative time to compute the solutions up to a given precision ε > 0 along an entire regularization path (several values of ν in decreasing order) and the memory space required by each sketching-based algorithm as measured by the sketch size m, and, (ii) the same criteria but for a fixed value of ν > 0.
We present in Figures 1 and 2 results for two standard datasets (see Appendix A for additional experiments): (i) one-vs-all classification of MNIST digits and (ii) one-vs-all classification of CIFAR10 images.
Except for very large values of the regularization parameter ν > 0 for which the regularized leastsquares problem (1) is well-conditioned so that the conjugate gradient is very efficient, we observe that our method is the fastest and requires less memory space than pCG for computing both the solutions of the entire regularization path and for a fixed value of ν. In particular, pCG uses m = d ρ for Gaussian embeddings and m = d log d ρ for the SRHT. Thus pCG is especially slow at the beginning because the factorization cost scales as O(d 3 ) and it requires memory space O(d 2 ). In contrast, our method starts with m = 1 and m does not exceed O(d e /ρ) for Gaussian embeddings and O(d e log d e /ρ) for the SRHT, as predicted by Theorem 5. Our adaptive sketch size actually remains significantly smaller than these theoretical upper bounds, and we still have a fast rate of convergence.
We observe in practice that, in Algorithm 1, the Polyak-IHS update is often rejected compared to the Gradient-IHS update, especially with the SRHT. Therefore, in addition to Algorithm 1, we consider a variant of it which does not consider the Polyak-IHS update but only the Gradient-IHS update. This variant enjoys the same convergence guarantees as presented in Theorems 5 and 6. Since it computes only a single candidate update, this variant is faster than Algorithm 1 in the case where the Polyak-IHS update is often rejected.   We fix the value of the regularization parameter ν = 10. We observe that pCG is slow at the beginning due to forming and factoring the m × d sketched matrix S · A with m ≈ d. In contrast, our methods start with m = 1 and the varying sketch size remains much smaller than that of pCG. On the other hand, we observe that the (per-iteration) convergence rate of pCG is faster and we believe this is due to using a larger sketch size (which is the best statistical lower bound known for this method). Each run is averaged over 30 independent trials.

A Additional numerical experiments
Here, we consider a synthetic dataset with A having exponential spectral decay σ j = 0.95 j for j = 1, . . . , d. The observation vector is generated as follows, b = Ax pl + η, where x pl is a planted vector with 1 √ d N (0, 1) independent entries and η is a vector of Gaussian noise 1 √ n N (0, I n ). We also consider the similar synthetic dataset but with polynomially decaying singular values σ j = 1/j for j = 1, . . . , d. Results are reported in Figure 3. For each algorithm, we start with the largest value ν = 1. For j 3, we initialize each algorithm at the previous solution x found for j + 1. For each value of ν, we stop the algorithm once ε = 10 −10 -precision is met. We observe that pCG is slow at the beginning due to forming and factoring the m × d sketched matrix S · A with m ≈ d. In contrast, our methods start with m = 1 and the varying sketch size remains much smaller than that of pCG. This leads to better time and memory space performance, except for the case of Gaussian embeddings and polynomial decays. In the latter case, our method is slowed down by Gaussian projections which are expensive. But with the SRHT, our method has the best performance. Each run is averaged over 30 independent trials. Mean standard deviations are reported in the form of error bars.

B The underdetermined case n d
A dual of the problem (1) is and one can map the optimal dual solution z * to the primal one using the relationship The dual problem fits into the primal overdetermined framework we considered in the main body of this manuscript. Note that there exists a priori the following technical difficulty. In order to formulate (7) as a regularized least-squares problem, one needs to consider the observation vector b = A † b (where A † is the pseudo-inverse of A) and then, we have One might wonder if b needs to be computed in order to apply the previous framework to the dual overdetermined case: this is not the case. Indeed, in Algorithm 1, the observation vector b only appears in the gradient formula, as ∇f (x t ) = A (Ax t − b). Equivalently, here we obtain That is, the gradient is easily computed and Algorithm 1 can be applied to the dual problem (7). Further, this yields guarantees on the error and using the dual-primal map (8), we would get guarantees on the error 1 2 x t − x * 2 . Regarding the semi-norm prediction error, we have where σ 1 is the largest singular value of A. Thus, even with the extra factor σ 2 1 which is negligible in our complexity analysis, we obtain the same guarantees as in Theorems 5, 6 and 7.
Note thatĀ =ŪΣV . Indeed, Further, the columns ofŪ are orthonormal, and the matrixΣ is diagonal with non-negative entries, so thatŪΣV is a singular value decomposition ofĀ.
Given an embedding S ∈ R m×n , denote byS the (m + d) × (n + d) block-diagonal matrix S 0 0 I d .
. We have thatĀ S SĀ = A S SA + ν 2 I d = H S . Consequently, given a step size µ ∈ R and a momentum parameter β ∈ R, the update formula (2) of the Polyak-IHS method can be equivalently written as Multiplying the update formula (9) byŪ Ā , subtractingŪ Ā x * , using the normal equationĀ b = A Ā x * and using the notation e t : =Ū Ā (x t − x * ), we obtain that e t e t−1 . (10)

C.1.1 Gradient-IHS method
For the Gradient-IHS method, we have that β = 0 so that the dynamics (10) simplifies to e t+1 = (I d − µC −1 S )e t . Using the fact that δ t = 1 2 e t 2 , we obtain that for any t 0, The eigenvalues of the matrix I d − µC −1 S are given by 1 − µ γi where the γ i 's are the eigenvalues of C S indexed in non-increasing order. Then, If λ, Λ > 0 are two real numbers such that λ γ d γ 1 Λ, then it holds that for any µ 0, where κ = Λ/λ, and this implies the claimed result.

C.1.2 Polyak-IHS method
Using (10) and the fact that δ t = 1 2 e t 2 , we immediately find by recursion that From Gelfand formula, we obtain that lim sup where ρ(M (µ, β)) is the spectral radius of the matrix M (µ, β) (the largest module of its eigenvalues), which remains to be computed. Let C S = T ΛT be an eigenvalue decomposition of the positive definite matrix C S -where Λ = diag(γ 1 , . . . , γ d ) and γ 1 . . . γ d > 0 -, and define the (2d) × (2d) permutation matrix Π as Then, it holds that

C.2 Proof of Theorem 5
Consider first Gaussian embeddings. Either the sketch size always remains smaller than a(ρ, η) · de ρ , which is equivalent to K log(a(ρ, η) de minitialρ )/ log (2). In this case, the statement of the theorem holds. Otherwise, there exists some iteration t 1 where m > a(ρ, η) · de ρ . Let t 1 be the first such iteration so that we also have m 2 · a(ρ, η) de ρ and K log(a(ρ, η) de minitialρ )/ log(2) + 1. Denote x t the current iterate, S the current sketching matrix, δ t = 1 2 A(x t − x * ) 2 and r t = According to Theorem 3, it holds with probability at least 1 − 16e −η 2 de/2 that In particular, it immediately follows from the above bounds that σmax(C S ) a(ρ, η).
The proof for the SRHT follows similar steps. Either the sketch size always remains smaller than a(ρ) · C(n, d e ) de log(de) ρ , which is equivalent to K log(a(ρ) · C(n, d e ) de log(de) minitialρ )/ log (2). In this case, the statement of the theorem holds. Otherwise, suppose that at some iteration t 1, we have m > a(ρ) · C(n, d e ) de log(de) ρ . Let t 1 be the first such iteration so that we also have m 2 · a(ρ) · C(n, d e ) de log(de) ρ and K log(a(ρ) · C(n, d e ) de log(de) minitialρ )/ log(2) + 1. Denote x t the current iterate, S the current sketching matrix, δ t = 1 2 A(x t − x * ) 2 and r t = According to Theorem 4, it holds with probability at least 1 − 9 de that In particular, it immediately follows from the above bounds that σmax(C S ) According to Theorem 1 and using that c gd (ρ/a(ρ)) = cgd(ρ) a(ρ) , the Gradient-IHS candidate x + gd of Algorithm 1 then satisfies δ + δt cgd(ρ) a(ρ) with probability at least 1 − 9 de , where we introduced the notation δ + : gd − x * ) 2 and using r + δ + σmin(C S ) and r t δt σmax(C S ) , it follows that r + r t Consequently, as soon as m > a(ρ)C(n, d e ) de log(de) ρ , it holds with probability at least 1 − 9 de that the Gradient-IHS update always verifies the improvement criterion r + rt c gd . Thus, future updates are not rejected and the sketch size remains constant.

C.3 Proof of Theorem 6
At every iteration, we observe at least the improvement r t r 0 max{c t gd ,c t p } =c t gd .
It holds that δ t σ max (C S ) · r t and δ 0 σ min (C Sinitial ) · r 0 so that Hence it suffices to control the ratio σmax(C S ) σmin(C S initial ) . For a Gaussian embedding, according to Lemma 2, we have that for any sketch size m m initial , with probability at least 1 − 8 · e −η 2 de/2 . We have (1 + de m ) 2 4 max{1, de minitial }. Using (1 + √ η) 2 4, this further yields For a SRHT embedding, according to Lemma 6, we have that for any sketch size m initial m d e , with probability at least 1 − 9/d e . On the other hand, it holds that for any embedding S, σ min (C S ) 1 − D 2 2 = ν 2 σ 2 1 + ν 2 . Thus, for Gaussian embeddings, we obtain and for SRHT embeddings,

C.4 Proof of Theorem 7
According to Theorem 6, we have for any t 1 It is immediate to verify that c gd = ρ for the SRHT. Therefore, the total number of iterations T to reach an ε-accurate solution is given by For ε min{ ν 2 σ 2 1 +ν 2 , 1 de }, this reduces to According to Theorem 5, the number of rejected steps satisfies At the beginning of the algorithm and for each rejected step, one needs to form the sketched matrix SA which costs O (nd log d e ), and, compute a factorization of H S as follows. According to the Woodbury matrix identity, the inverse of H S verifies To reduce the complexity of solving at each iteration the linear system H S · z = ∇f (x t ), one can simply cache the matrix (ν 2 I m + SA(SA) ) −1 which takes time O( d 2 e log 2 de ρ 2 d) to compute. Since K = O(log(d e /ρ)), it follows that the total sketching and factor costs scale as O(log(d e /ρ) · ( d 2 e log 2 de ρ 2 d + nd log(d e /ρ))).
The per-iteration cost is that of computing the matrix-vector products Ax t and A (Ax t − b), which is given by O(nd). Note that the other main numerical operation consists in solving the linear system H S · z = ∇f (x t ). Using the cached matrix (ν 2 I m + SA(SA) ) −1 and the identity Thus, we obtain the total time complexity which is the claimed result. By definition, we have On the other hand, we have that Observing that U S S U = C S (see the proof of Theorem 1 for the proof of this statement) and recalling the definitions E t = C − 1 2 S U A(x t − x * ) and r t = 1 2 E t 2 , it follows that 1 2 g t H −1 S g t = 1 2 E t 2 = r t , which concludes the proof.

E Proofs of concentration inequalities E.1 Gaussian concentration over ellipsoids -Proof of Theorem 3
Fix a sketch size m 1, define ρ = de m , and let S ∈ R m×n be a Gaussian embedding, i.e., a random matrix with i.i.d. entries N (0, 1/m). Let γ 1 and γ d be respectively the largest and smallest eigenvalue of the matrix C S = D(U S SU − I d )D + I, i.e., We show separately the following two results which, combined together, yields Theorem 3.
Lemma 2. For any η ∈ (0, 1), it holds with probability at least 1 − 8 · e −mρη 2 /2 that Lemma 3. Suppose that m > d e . Then, for any η ∈ (0, (1 − √ ρ) 2 /4), it holds with probability at where c(η) : = We prove Lemma 2 in Appendix E. Suppose that m 1, and denote, as previously, ρ : = de m . Define the matrixD = D D 2 . Note that D 2 F = d e and D 2 = 1. We aim to study the upper tail of the largest eigenvalue γ 1 of the matrix C S , which satisfies the variational equation where we introduced the random matrix G ∈ R m×d with i.i.d. Gaussian entries N (0, 1) and the first equality holds since SU d = 1 √ m G. We also use the notations C = D x | x = 1 and ψ(u, z) : = − 1 2 ( u 2 + m z 2 ). We introduce the auxiliary random variable where g ∈ R m and h ∈ R d are random vectors with i.i.d. entries N (0, 1). Using Theorem 9 (see Appendix E.1.3), it holds that for any c ∈ R, Consequently, it suffices to control the upper tail of Y in order to control that of γ 1 . First, we recall a few basic facts on the concentration of Gaussian random vectors (see, for instance, Theorems 3.1.1 and 6.3.2 in [28]). That is, for any η ∈ (0, 1), the following event holds with probability at least 1 − 4e −mρη 2 /2 , On the event E η , we have In equality (i), we used the fact that for a vector u with fixed norm u = t, the maximum of g u is equal to g t. In inequality (ii), we bounded h z by |h z| and then relaxed the constraint t 0 to t ∈ R. In inequality (iii), we plugged-in the value of the maximizer t * = z g + |h z|. In inequality (iv), we used the fact that for z ∈ C, z 1. In (v), we used the change of variable z =Dx with x = 1. In (vi), we used the fact that for x = 1, D h, x D h and the equality is attained for some x, and then, we used the changed of variable t = | D h, x |. In (vii), we used the fact that on E η , | g 2 −m| m η, D h √ d e (1 + √ η) and g √ m(1 + √ η). In (viii), we plugged into the objective the value of the maximizer t * = √ d e (1 + √ η). Equality (ix) follows from the definition of ρ = d e /m. Inequality (x) follows from the fact that for any η > 0, we have Thus, we have that 8 · e −mρη 2 /2 , which is the claimed inequality (11).

E.1.2 Controlling the smallest eigenvalue γ d (proof of Lemma 3)
Suppose that m > d e , and denote, as previously, ρ : = de m andD : = D D 2 . Note that D 2 F = d e and D 2 = 1. We aim to study the lower tail of the smallest eigenvalue γ d of the matrix C S , which satisfies the variational equation where we introduced the random matrix G ∈ R m×d with i.i.d. Gaussian entries N (0, 1) and the first equality holds since SU d = 1 √ m G. We also use the notations C = D x | x = 1 and ψ(u, z) : = − 1 2 ( u 2 + m z 2 ). We introduce the auxiliary random variable where g ∈ R m and h ∈ R d are random vectors with i.i.d. entries N (0, 1). Using Theorem II.1 from [25], it holds that for any c ∈ R, where, in the first inequality, we relaxed the constraint set by removing the constraint z g +h z 0 and we used the fact that z 1. In the second inequality, we used the change of variable z =Dx with x = 1. In the third inequality, we used the fact that | D h, x | D h and used the change of variable | D h, x | = t with t ∈ [0, D h ]. On the event E η , it follows that , further algebra yields that, on the event E η , Hence, combining all of the above, we find that, for any η ∈ (0, 1/2), with probability at least 1 − 8e −mρη 2 /2 ,

E.1.3 A new Gaussian comparison inequality
We start with the following well-known comparison inequality, which was first derived in [9].
Theorem 8 (Gordon's Gaussian comparison theorem). Let I, J ∈ N * , and {X ij }, {Y ij } be two centered Gaussian processed indexed on I × J, such that for any i, l ∈ I with i = l and j, k ∈ J, Then, for any {λ ij } ∈ R I×J , we have Lemma 4. Let G ∈ R m×n , Z ∈ R, g ∈ R m and h ∈ R n have independent standard Gaussian entries. Let I 1 ⊂ R n and I 2 ⊂ R m be finite sets, and ψ be a function defined over I 1 × I 2 . Then, for any c ∈ R, we have P max (x,y)∈I1×I2 y Gx + Z x y + ψ(x, y) c P max (x,y)∈I1×I2 x g y + y h x + ψ(x, y) c .
Proof. We introduce two Gaussian processes X and Y indexed over I 1 × I 2 , defined as X xy = x g y + y h x , Y xy = y Gx + Z x y , for all (x, y) ∈ I 1 × I 2 . It holds that EX xy = EY xy = 0, EX 2 xy = 2 x 2 y 2 = EY 2 xy , and E[X xy X x y ] = x x y y + y y x x , Therefore, applying Gordon's comparison theorem with I = I 1 × I 2 , J being any finite set, and λ xy = ψ(x, y) − c, we obtain that P min (x,y)∈I1×I2 y Gx + Z x y − ψ(x, y) −c P min (x,y)∈I1×I2 x g y + y h x − ψ(x, y) −c .
Using the symmetry of the Gaussian distribution, it follows that P max (x,y)∈I1×I2 y Gx + Z x y + ψ(x, y) c P max (x,y)∈I1×I2 x g y + y h x + ψ(x, y) c , and consequently, P max (x,y)∈I1×I2 y Gx + Z x y + ψ(x, y) c P max (x,y)∈I1×I2 x g y + y h x + ψ(x, y) c , Corollary 1. Let S 1 ⊂ R n and S 2 ⊂ R m be non-empty sets, and ψ : S 1 × S 2 → R be a continuous function. Then, for any c ∈ R, P sup (x,y)∈S1×S2 y Gx + Z x y + ψ(x, y) c P sup (x,y)∈S1×S2 x g y + y h x + ψ(x, y) c , Proof. According to Lemma 4, the result is true if S 1 and S 2 are finite. By monotone convergence, it is immediate to extend it to countable sets. By density arguments and monotone convergence, it also follows for any sets S 1 and S 2 .
Consequently, using the independence of f 1 and Z, we get x g y + y h x + ψ(x, y) c , which yields the claim.

E.2 SRHT matrices -matrix deviation inequalities over ellipsoids
Let S ∈ R m×n be a SRHT matrix, that is, S = RHdiag(ε) where R is a row-subsampling matrix of size m × n, H is the normalized Walsh-Hadamard transform of size n × n and ε is a vector of n independent Rademacher variables. We introduce the scaled diagonal matrixD = D D 2 . Note that D 2 F = d e and D 2 = 1. where E : = diag(e j H). Each entry of E has magnitude n − 1 2 . The function f is convex, and its Lipschitz constant is upper bounded as follows, For a Rademacher vector ε, we have Applying Lipschitz concentration results for Rademacher variables, we obtain P e j Hdiag(ε)UD d e n + 8 log(βn) n 1 nβ .
Theorem 10 (Matrix Bernstein). Let X = {X 1 , . . . , X n } be a finite set of squared matrices with dimension d. Fix a dimension m, and suppose that there exists a positive semi-definite matrix V and a real number K > 0 such that E[X I ] = 0, E[X 2 I ] V , and X I 2 K almost surely, where I is a uniformly random index over {1, . . . , n}. Let T be a subset of {1, . . . , n} with m indices drawn uniformly at random without replacement. Then, for any t m V 2 + K/3, we have P i∈T X i 2 t 8 · d e · exp − t 2 /2 m V 2 + Kt/3 , log E e θXi j   = tr exp m log E e θX I .
For the last inequality, we used the fact that e a e a −a−1 = 1 + 1+a e a −a−1 1 + 3 a 2 for all a 0. Picking θ = t/(m V 2 + Kt/3), we obtain Under the assumption t m V 2 + K/3, the parenthesis in the above right-hand side is bounded by four, which results in