{'fraction_non_alphanumeric': 0.036727428404046775, 'fraction_numerical': 0.011896982943151349, 'mean_word_length': 4.474528185518688, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularized Momentum Iterative Hessian Sketch for Large Scale Linear System of Equations

In this article, Momentum Iterative Hessian Sketch (M-IHS) techniques, a group of solvers for large scale linear Least Squares (LS) problems, are proposed and analyzed in detail. The proposed techniques are obtained by incorporating the Heavy Ball Acceleration into the Iterative Hessian Sketch algorithm and they provide significant improvements over the randomized preconditioning techniques. Through the error analyses of the M-IHS variants, lower bounds on the sketch size for various randomized distributions to converge at a pre-determined rate with a constant probability are established. The bounds present the best results in the current literature for obtaining a solution approximation and they suggest that the sketch size can be chosen proportional to the statistical dimension of the regularized problem regardless of the size of the coefficient matrix. The statistical dimension is always smaller than the rank and it gets smaller as the regularization parameter increases. By using approximate solvers along with the iterations, the M-IHS variants are capable of avoiding all matrix decompositions and inversions, which is one of the main advantages over the alternative solvers such as the Blendenpik and the LSRN. Similar to the Chebyshev Semi-iterations, the M-IHS variants do not use any inner products and eliminate the corresponding synchronizations steps in hierarchical or distributed memory systems, yet the M-IHS converges faster than the Chebyshev Semi-iteration based solvers.

1. Introduction. We are presenting a group of solvers, named as Momentum Iterative Hessian Sketch, M-IHS, that is designed for solving large scale linear system of equations in the form of (1.1) where A ∈ R n×d is the given data or coefficient matrix, b is the given measurement vector contaminated by the noise or computation/discretization error ω, and x 0 is the vector desired to be recovered. Due to contaminated measurements, solutions can differ significantly according to the constraints imposed on the problem. In this article, we are particularly interested in the ℓ2-norm regularized Least Squares (LS) solution: x * = argmin , (1.2) which is known as the Tikhonov Regularization or the Ridge Regression. The problem in (1.2) frequently arises in various large scale applications of science and engineering. For example, it can appear in the discretization of Fredholm Integral Equations of the first kind [28]. In those cases, the data matrix might be ill conditioned and the linear system might be either over-determined or square. When the system is under-determined, although sparse solutions are more popular due to the Compressed Sensing [16], the least norm solutions occupy an important place in statistic applications such as the Support Vector Machines [55,13]. Solutions to the problem in both regimes, i.e, n ≥ d and n < d, are often required as intermediate steps of rather complicated algorithms such as the Interior Point and the ADMM that are widely used in machine learning and image processing applications [9,11,47].
Throughout the manuscript, we are going to assume that a proper regularization parameter λ estimate is available. In terms of complexity and error performance, estimating a proper regularization parameter is equally important as obtaining the regularized solution. Risk estimators such as the Discrepancy Principle, Unbiased Prediction Risk Estimate, Stein's Unbiased Risk Estimate and Generalized Cross Validation can be directly used to estimate the regularization parameter of the moderate size problems [52]. For large scale problems, these risk estimators can be adapted for the lower dimensional sub-problems that arise during the iteration of the first order iterative solvers [34]. A hybrid scheme that adaptively selects the regularization parameter along with the iterations is also suitable for the proposed M-IHS solvers. Indeed, we have developed such a technique, but we will present that study in a separate manuscript to keep the length of this document in a reasonable size.
The regularized solution in (1.2) can be obtained by using direct methods such as the Cholesky decomposition for square A, or the QR decomposition for rectangular A. However, O(nd min(n, d)) computational complexity of any full matrix decomposition becomes prohibitively large as the dimensions increase. For large scale problems, linear dependence on both dimensions is acceptable and can be obtained by using the first order iterative solvers based on Krylov Subspaces [8]. These methods require only a few matrix-vector and vector-vector multiplications for each iteration, but the number of iterations that is needed to reach a certain level of tolerance is highly sensitive to the condition number of the coefficient matrix. If the largest and the smallest singular values of A are known, the optimal and un-improvable convergence rate of the first order iterative solvers is O(1/k 2 ) which can be obtained by using Nesterov's Accelerated Gradient Descent or Polyak's Heavy Ball Method (HBM), unfortunately such information on the largest and the smallest singular values of A is rarely available in practice [40,48]. In the absence of this information, the Conjugate Gradient (CG) technique achieves the same rate by adaptively tuning the momentum parameters through additional calculations at each iteration [12]. For the problem in (1.2), techniques such as the LSQR and LSMR, that are based on Golub-Kahan-Lanczos Bidiagonalization, produce more stable results than the CG technique with the same convergence rate [44,23]. Other techniques that are based on the Krylov subspace approach with different convergence behaviours can be added to this list [5,25]. When these techniques are applied to the problem in (1.2), their common convergence rate is characterized by the following inequality: where x * is the optimal solution of (1.2), x 1 is the initial guess, x i is the i-th iterate of the solver and the condition number κ(·) is defined as the ratio of the largest singular value to the smallest singular value of its argument. Since for ill conditioned matrices κ(A T A + λI d ) becomes large, the rate of convergence can be extremely slow.
The computational complexity of Krylov subspace-based iterative solvers is O(nd) for each iteration, which is significantly less than O(nd min(n, d)) if the number of iterations can be significantly fewer than min(n, d). However, in applications such as big data where A is very large dimensional, the computational complexity is not the only metric for feasibility of the algorithms. For example, if the coefficient matrix is too large to fit in a single working memory and it could be merely stored in a number of distributed computational nodes, then at least two distributed computations of matrix-vector multiplications are required at each iteration of algorithms such as the CGLS or the LSQR (please see [31,6] for such applications). This suggests that the number of iterations should also be counted as an important metric to measure the overall complexity of an algorithm. One way to reduce the number of iterations in the iterative solvers is to use preconditioning which is a linear mapping of the solution domain for transforming an ill conditioned problem to a well conditioned problem. In the deterministic settings, finding a low-cost and effective preconditioning matrix is still a challenging task.
In addition to the number of iterations, the number of inner products at each iteration also plays an important role in the overall complexity. Each inner product calculation constitutes a synchronization point in parallel computing and therefore is undesirable for distributed or hierarchical memory systems [5]. The Chebyshev Semi-iterative (CS) technique can be preferred in this kind of applications, since it does not use any inner products and therefore eliminates some of the synchronization steps that are required by the techniques such as the CG or the GMRES. However, the CS requires prior information about the ellipsoid that contains all the eigenvalues of A, which is commonly not available in practice [26].
The Random Projection (RP) techniques that are based on the Johnson-Lindenstrauss Lemma have found diverse field of applications [30,17]. These algorithms are capable of both reducing the dimensions and bounding the number of iterations with statistical guarantees, while they are convenient for parallel and distributed computations as much as the direct methods. The development and the applications of the RP based algorithms can be found in [54,35] and references therein.
In the application of the RP to the regularized LS problem in (1.2), there are two main approaches. In the first approach, that is referred to as classical sketching, the coefficient matrix A and the measurement vector b are projected down onto a lower dimensional (m ≪ n) 1 subspace using a randomly constructed sketch matrix S ∈ R m×n . Then, the aim is to obtain an ζ-optimal solution with high probability in terms of the cost approximation [18,45]: The best lower bounds on the sketch size for obtaining an ζ-optimal cost approximation have been derived for both sparse and dense sketch matrices by Avron et al. [1]. They showed that the sketch size can be chosen proportional to the statistical dimension, which is defined as sd λ (A) = Tr A(A T A + λI d ) −1 A T . Although the cost approximation is sufficient for many machine learning problems, small distance to the optimal solution, which is defined as the solution approximation, is a more preferable metric for the problems arisen from, for example, discretization of Fredholm integrals in the applied linear algebra [8,52]. However, in [46], the classical sketching is shown 1 Without loss of generality we can assume the linear system is over-determined, if it is not, then we can take a dual of the problem to obtain an over-determined problem as examined later.
to be sub-optimal in terms of sketch size for obtaining a solution approximation.
In the second approach, that is referred to as randomized preconditioning, the algorithms with reasonable sketch sizes obtain an η-optimal solution approximation: where W is a positive definite weight matrix, by iteratively solving a number of low dimensional sub-problems constituted by (SA, ∇f (x i )) pairs. To the best of our knowledge, Rokhlin [49] is the first who uses random projection techniques proposed in [18] to construct a preconditioning matrix for CG-like algorithms. He used the inverse of R factor in the QR decomposition of the sketched matrix SA for this purpose. Later, implementation of similar ideas resulted in Blendenpik and LSRN which have been shown to be faster than some of the deterministic solvers of LAPACK [3,36]. To solve the preconditioned problems, as opposed to the Blendenpik which uses the LSQR, the LSRN uses the CS technique for parallelization purposes and deduce the prior information about the eigenvalues by using results from the random matrix theory. The main drawback of the LSRN and Blendenpik is that regardless of the desired accuracy η, one has to pay the whole cost, O(md 2 ), of a full m × d dimensional matrix decomposition which is the dominant term in the computational complexity of both algorithms. Iterative Hessian Sketch (IHS) proposed in [46] eliminates the dominant term, O(md 2 ), by using the inverse of the sketched Hessian as preconditioning matrix in the Gradient Descent method [53]. Therefore, instead of computing a full decomposition or an inversion, a linear system can be approximately solved for a pre-determined tolerance 2 . By adapting this idea into the CG technique, Accelerated IHS (A-IHS) has been obtained in [53]. Lastly, in [42], it has been showed that if the linear system is strongly over-determined, then the momentum parameters of the HBM can be robustly estimated by using Marchenko Pastur (MP) Law [19,20]. This analysis results in a prototype solver M-IHS that we study here in detail. The statistical lower bounds in the current literature suggest that the sketch size in randomized preconditioning algorithms can be chosen proportional to the rank of the problem which can be much larger than the statistical dimension. Although, some lower bounds on the sketch size that are proportional to the statistical dimension have been obtained in Kernel Ridge Regression [2], it has not obtained in regularized LS problem yet.

Contributions.
In this article, we describe a group of random projection based iterative solvers for large scale regularized LS problems. The proposed M-IHS variants can be used for any dimension regimes if the statistical dimension of the problem is sufficiently smaller than both size of the coefficient matrix A. In section 5, we give the detailed convergence analyses of the proposed techniques. Our guarantees, presented in Theorem 3.1 and Corollary 3.3, are based on the solution approximation metric as opposed to the results obtained for cost approximation metric in [1]. In Corollary 3.2, we derived the best lower bounds, in the current literature, on the sketch size of various randomized distribution for attaining a certain convergence rate with a constant probability. These guarantees can be easily extended to any other sketch types by using the Approximate Matrix Multiplication (AMM) property [15]. If tighter bounds are acquired for the AMM property in the future, our bounds can be automatically improved as well. Although our bounds for the dense sketch matrices such as Subgaussian or Randomized Orthonormal Systems (ROS) are the same as in [1], we gained slightly better results for the sparse sketch matrices. Additionally, we provide some empirical bounds for the sketch size and the rate of convergence in Corollary 3.4 which is very tight as demonstrated through various numerical experiments. Lastly, in Algorithm 6.1, we extend the idea of LSQR into the linear problems in the form of A T Ax = b which we need to solve during the iterations of all proposed Inexact M-IHS variants and of the Newton Sketch [47]. The advantages that the proposed sub-solver, referred to as AAb Solver, provides over the symmetric CG technique are the same as the advantages that are provided by the LSQR over the CGLS technique [44]. The MATLAB implementations of the proposed solvers together with the codes that generate the figures in the article can be found in the following link: https://github.com/ibrahimkurban/M-IHS.
3. The proposed solvers for regularized LS problems. The M-IHS is obtained by combining the IHS technique with the Heavy Ball Acceleration. The IHS technique iteratively minimizes the quadratic objective in (1.2) by performing the following updates: In all iterations, the same sketch matrix can be used by adding a proper step size and the convergence can be accelerated by an additional momentum term as: If we restrict our attention to the fixed parameters, the optimal momentum parameters α and β, that maximize the convergence rate, can be estimated by using random matrix theory as completely independent of spectral properties of the coefficent matrix A. Here, the linear system is assumed to be strongly over-determined, i.e., n ≫ d, but by using the dual problem, the theory can be easily extended to the strongly under-determined case of d ≫ n as well [10]. A dual of the problem in (1.2) is , and the relation between the solutions of the primal and dual problem is The dual problem in (3.2) can be minimized by using the same approach as the M-IHS: . Then, the M-IHS applied on (1.2) and the Dual M-IHS applied on (3.2) with fixed momentum parameters converge to the optimal solutions, x * and ν * , respectively, at the following rate 3 with probability at least 1 − δ: where D −1 is the diagonal matrix whose diagonal entries are σ 2 i + λ and σ i is the i-th singular value of A.  (i) Sparse Subspace Embedding [54] with single nonzero element in each column, sketch size [33,39], size

non-zero elements in each column
and SA is computable in O(s · nnz(A)); (iii) SRHT sketch matrix [15,46] with size m = Ω (sd λ (A) + log(1/ǫδ) log(sd λ (A)/δ)) /ǫ 2 and SA is computable in O(nd log(m)); [45,15] with size The proof of Corollary 3.2 can be found in subsection 5.2. For the solution approximation, we do not need the second condition in Lemma 11 of [1], hence we obtained slightly better results for the sparse subspace embeddings in item (i) and (ii) of Corollary 3.2. The number of iterations to reach a certain level of accuracy is stated in the following corollary.
will converge to the optimal solutions x * and ν * , respectively, with a convergence rate of √ β as m → ∞ while min(n, d)/m remains constant. Any sketch size m > sd λ (A) can be chosen to obtain an η-optimal solution approximation at most Remark 3.5. Although the MP law provides bounds for the singular values of the sketch matrix S in asymptotic regime, i.e., as m → ∞; the bounds become very good estimators of the actual bounds when m takes finite values. As shown in Figure 1, the rate of √ β creates a remarkable fit for the numerical convergence rate of the M-IHS variants when the momentum parameters given in Corollary 3.4 are used for the Tikhonov regularization. This is because the sigmoid-like filtering coefficients in the Tikhonov regularization can be thought of as the smoothed version of the binary coefficients in the TSVD solution and therefore the binary coefficients constitute a good approximation for the filtering coefficients of the Tikhonov regularization as noted in subsection 5.4.
In practice, the M-IHS and Dual M-IHS eliminate all the quadratic terms in the complexity expression by approximately solving the low dimensional linear systems in (3.1) and (3.4) instead of computing a matrix decomposition or an inversion. This inexact sub-solver approach constitutes a trade-off between the computational complexity and the accuracy, that is highly desirable in very large dimensional problems where solutions with relatively lower accuracy are acceptable. Unfortunately, forming this trade-off is not possible for the Blendenpik and the LSRN techniques. Inexact sub-solvers have been known to be a good heuristic way to create this trade-off and they are widely used in the algorithms that are based on the Newton Method to solve the large scale normal equations [41]. In these inexact (or truncated) Newton Methods, inner iterations are terminated at the moment that the relative residual error is below some iteration dependent threshold, named as the forcing terms [21]. In the literature, there are various techniques to choose these forcing terms that guarantee global convergence, but the number of iterations suggested by these techniques are far above the total number of iterations used in practice. We refer interested reader to [37] and we go forward with the heuristic constant threshold, ǫ sub , that checks the relative residual error of the linear system [7].
Any Krylov subspace techniques can be used to solve the sub-problems in (3.1) and (3.4), but LSQR-like solvers that are adapted for the normal equations would require computations of 4 matrix-vector multiplications per iteration. On the other hand, due to the explicit calculation of (SA) T (SA)z, the symmetric CG, that would require only 2 matrix-vector multiplications, might be unstable in the ill-conditioned problems [44]. Therefore, in section 6, we propose a stable sub-solver which is particularly designed for the problems in the form of A T Ax = b. The proposed sub-solver, referred to as AAb Solver, is based on the Golub Kahan Bidiagonalization and it uses a similar approach that the LSQR uses on the LS problem. In addition to the stability advantage over symmetric CG technique, AAb Solver produces a bidiagonal representation of sketched matrix as a byproduct of the iterations. This bidiagonal form can be used in, for example, Generalized Cross Validation [24,34] to estimate the problem related parameters including the regularization parameter and the statistical dimension. The inexact versions of the M-IHS and Dual M-IHS that use AAb Solver are given in Algorithm 3.1 and Algorithm 3.2 where RP fun represents the function that generates the desired sketched matrix such that E S T S = I m whose implementation details can be found in the relevant references in Corollary 3.2.Setting the Corollary 3.2 suggests that if the statistical dimension is several times smaller than the dimensions of A, then it is possible to choose a substantially smaller sketch size than the min(n, d). If this is the case, then the quadratic objective functions in (3.1) and (3.4) become strongly under-determined problems, which makes it possible to approximate the Hessian of the objective functions one more time by taking their convex dual as it has been done in the Dual M-IHS. This approach is similar to approximately solving the problems in (3.1) and (3.4) by using the AAb Solver, but this time we are going to apply a second dimension reduction. At the end of two Hessian sketching, the linear sub-problem whose dimensions are reduced from both sides can be efficiently solved by the AAb Solver for a pre-determined tolerance as before. Consider the following dual of the sub-problem in (3.1) , which is a strongly over-determined problem if m ≪ min(n, d). Hence, it can be approximately solved by the M-IHS updates: After M iterations, the solution of (3.1) can be recovered by using the relation in are all degraded in the LS problems with high condition numbers due to the instability issue of the symmetric CG technique [44]. Even if the regularization is used, still the performance of the solvers proposed in [53] are considerably deteriorated compared to the other randomized preconditioning techniques as shown in section 4. Further, applying the preconditioning idea of IHS to the stable techniques such as the LSQR that are adapted for the LS problem is not so efficient as the M-IHS variants, because two preconditioning systems are needed to be solved per iteration for such techniques.
The computational saving when we apply a second dimension reduction as in the Primal Dual M-IHS may not be significant due to the second gradient computations in Line 10 and 9 of the given algorithms, but the lower dimensional sub-problems that we obtain at the end of the second sketching can be used to estimate several parameters including the regularization parameter itself. Indeed, we are currently working on for j=1:M do 10: ∆z j = AAb Solver(WAS T , g i,j , λ, ǫ sub ) 12: end for 14: such an algorithm to estimate the regularization parameter of the large scale problems along with the iterations of the M-IHS variants without requiring additional accesses to coefficient matrix A.
The Primal Dual M-IHS techniques are extension of the inexact schemes. Therefore, their convergence rates depend on their forcing terms that are used to stop the inner iterations [37]. In [53], an upper bound for the error of the primal dual updates is proposed. However as it is detailed in section 8, there are several inaccuracies in the development of the bound. Therefore, finding a provably valid lower bound on the number of inner loop iterations, that guarantee a certain rate of convergence at the main loop, is still an open question for the primal dual algorithms.
The statistical dimension in Algorithm 3.1, 3.2, 3.3 and 3.4 can be estimated by using a Hutchinson-like randomized trace estimator. We refer interested reader to [4] for a detailed comparison of Hutchinson-like estimators. Alternatively, an algorithm is proposed in [1]  However, due to the third order root and the division by at least a sixth order polynomial, we can only have such small statistical dimension value when the singular values of A decay severely/exponentially and the signal-to-noise ratio is very low. Therefore, we preferred to use the heuristic trace estimator in Algorithm 3.5. The input matrix SA can be replaced with SA T or even with WA T S T and WAS T according to the algorithm used. Any estimator in [4] can be substituted for the Hutchinson Estimator and the number of samples, N , can be chosen accordingly. In the experiments that we made with various singular value profiles, small samples sizes such as 2 or 3 and ǫ tr = 0.01 is sufficient to obtain satisfactory results. Noting that the momentum pa- for j=1:M do 10: ∆z j = AAb Solver(WA T S T , g i,j , λ, ǫ sub ) 12: end for 14: 15: τ = τ + x T z i 8: end for 9: Output: τ = τ /N rameters suggested in Corollary 3.4 are not sensitive to small errors in sd λ (A), since β = sd λ (A)/m. When the sketch size m exceeds thousands, for example, an error made at the first digit of the statistical dimension estimation, distorts the momentum parameter only at the third decimal.

Numerical Experiments and Comparisons.
For a fair comparison, we have implemented all the proposed algorithms in this manuscript as well as those that are used for comparisons in MATLAB. All codes can be found in the link 4 . The coefficient matrix A ∈ R n×d is generated for various sizes as follows: we first sample the entries of A from distribution N (1 d , Γ) where Γ ij = 5 · 0.9 |i−j| so that the columns are highly correlated with each other. Then, by using the SVD, we replace the singular values with philips profile provided in RegTool. We scale the singular values to set the condition number κ(A) to 10 8 and we use the same input signal provided by RegTool [27]. In this way, we have obtained a challenging setup for any first order iterative solvers to compare their performances. In all experiments, the same setup has been used unless indicated. We counted the number of operations according to Hunger's report [29]. All results have been obtained by averaging over 32 MC simulations.
We compared the proposed algorithms with the state of the art randomized preconditioning techniques which can reach any level of desired accuracy within a bounded number of iterations. The compared methods can be briefly described as follows. The Blendenpik uses the R factor in QR decomposition of the sketched matrix SA as the preconditioning matrix for the LSQR algorithm just like the method proposed by Rokhlin et al [3,49] and it uses Randomized Orthonormal System (ROS) to generate the sketched matrix [45]. The LSRN uses the V factor in the SVD similar to the Blendenpik. In spite of its high running complexity, the Gaussian sketch matrices are preferred in the LSRN. In addition to the LSQR, the CS can be used in the LSRN as the core solver for parallelization purposes without calculating the singular values explicitly [36]. The IHS uses the sketched Hessian as the preconditioning matrix for the Gradient Descent. The Accelerated IHS (A-IHS) uses this idea for the CG algorithm in over-determined problems. The dual counter-part of the A-IHS algorithm, A-IDRP, is shown to be faster than the Dual Random Projection algorithm proposed in [55], so we did not include the DRP in the simulations. Additionally, we include a CS variant of the IHS (IHS-CS) to the comparisons. We combined the randomized preconditioning idea of IHS with the preconditioned CS method [5]. We found the bounds for the eigenvalues in the same way as the LSRN. We have solved the low dimensional sub-problems required by all IHS variants by taking the QR decomposition, but for inexact schemes, we have used the proposed AAb Solver with constant forcing term ǫ sub = 0.1. We did not include inexact versions of the accelerated algorithms proposed in [53], since their exact versions are outperformed by the Exact M-IHS variants in all settings. Except for the LSRN variants which use Gaussian sketch matrices, we used Discrete Cosine Transform in the ROS for all techniques.
We compare the operation counts required by the algorithms to obtain a certain level of accuracy for the solution approximation metric. In the first experiment, we did not include noise in the linear system to emphasize the convergence rate that the algorithms can provide in such severely ill posed problems. To make the problem more challenging, we sampled the input vector, x 0 , from uniform distribution Uni(−1, 1) for this experiment only. In such scenarios, convergence rates of Krylov subspace-based iterative solvers without preconditioning fall its minimum value. The results are shown in Figure 2. Due to high running time of the Gaussian sketches, O(mnd), the LSRN variants require more operations (for the size of the problems considered here approximately 10 times larger) than the others. Due to the lack of inner product calculations, the M-IHS requires slightly fewer operations than the Blendenpik, nonetheless it reaches the same accuracy with the LSRN-LSQR. The A-IHS algorithm has the lowest performance which is expected in the un-regularized problems, since it is adapted on the CG technique that can be unstable for the unregularized LS problems due to the high condition number [44]. The convergence of the CS-based techniques, both of the IHS and the LSRN variants, are substantially slower than the M-IHS, which suggests that the M-IHS algorithm can take the CSs place in those applications where parallel computation is viable. A similar comparison of the M-IHS with Accelerated Randomized Kaczmarz (ARK) and CGLS without preconditioning can be seen on Figure 2 of [42].
By using additive i.i.d Gaussian noise at various levels of variance, we tested the aforementioned algorithms on regularized LS problems. Results for over-determined and under-determined cases can be seen on Figure 3 and Figure 4, respectively. We Corollary 3.2 implies eligibility of sketch sizes that are smaller than the rank, m ≤ min(n, d) which suggest that the RP techniques can be used for all dimension regimes whether it is strong or not, as long as the statistical dimension of the problem is small compared to the dimensions of coefficient matrix A. This implication can be verified in Figure 5 on which we showed the performance of the Primal Dual M-IHS techniques. As discussed in section 3, the Primal Dual M-IHS techniques seem to increase the complexity in the contrary of saving, but they obtained lower dimensional sub-problems than the M-IHS and Dual M-IHS, which suggests that some problem related parameters can be deduced more cheaply during iterations of the Primal Dual M-IHS technique. Lastly, the Primal Dual M-IHS variants have a noticeably higher rate of convergence than the A-IPDS algorithm which is based on the CG technique.
5. Theoretical analysis of proposed approaches. In this section, we are going to provide a unified analysis of the regularized LS problems through the proofs by including many comments and explanations. Throughout the analysis we denote A = U ΣV T as the compact SVD with U ∈ R n×r , Σ = diag(σ 1 , . . . , r) ∈ R r×r and V ∈ R d×r where r = min(n, d).

Proof of Theorem 3.1.
To prove the theorem for the M-IHS and the Dual M-IHS, we mainly combine the idea of partly exact sketching, that is proposed in [1], with the Lyapunov analysis, that we use in [42]. In parallel to [1], we define diagonal matrix D := (Σ 2 + λI r ) −1/2 and the partly exact sketching matrix as:

The Proof for M-IHS:. Let
so that U 1 is the first n rows of an orthogonal basis for A √ λI d as required by the condition in (3.5) of the theorem. We have the following equality Since the following Hessian Sketch (HS) sub-problem is equivalent to (3.1), we can examine the following bipartite transformation to find out the convergence properties of the M-IHS: The eigenvalues and therefore the contraction ratio of the transformation can be found analytically by converting the matrix T into a block diagonal form through the same similarity transformation given in [42]: where ΨΛΨ T is the Spectral decomposition of ( A T S T S A) −1 and λ i is the i th eigenvalue. The characteristic polynomials of each block is If the following condition holds then both of the roots are imaginary and both have a magnitude √ β for all λ's. In this case, all linear dynamical systems driven by above characteristic polynomial will be in the under-damped regime and the contraction rate of the transformation T , through all directions, not just one of them, will be exactly √ β. If the condition in (5.1) is not satisfied for a λ j with j ∈ [r], then the linear dynamical system corresponding to λ j will be in the over-damped regime and the contraction rate in the direction through the eigenvector corresponding to this over-damped system will be smaller compared to the others. As a result, the overall algorithm will be slowed down substantially (see [43] for details). If the condition in (3.5) of Theorem 3.1 holds, Consequently, the condition in (5.1) can be satisfied for all λ i 's by the following choice of β that maximizes the convergence rate over step size α where the minimum is achieved at α * = The Proof for Dual M-IHS:. The proof for the dual counter-part is almost same as the proof for M-IHS except for the following modifications. Let so that U 1 is the first d rows of an orthogonal basis for A T √ λI n as required by the theorem. We have the following equality: for ∀ν ∈ R n | w = D −1 U T ν . Thus, following HS sub-problem is equivalent to (3.4) ∆w i+1 = argmin We can analyze the following bipartite transformation to figure out the convergence properties of the Dual M-IHS The rest of the proof can be completed straightforwardly by following the same analysis steps as in the proof for the M-IHS case.

5.2.
Analysis of the Condition in Theorem 3.1. Theorem 3.1 is valid if the sketch matrix S satisfies the following condition: and U 1 2 2 = σ 2 1 σ 2 1 +λ ≈ 1 for a properly chosen regularization parameter λ. In this subsection before giving the proof of Corollary 3.2, we summarize some of the current results and tools from the random matrix theory literature that is used to prove the condition in (5.2) for different applications such as low rank approximation, matrix ridge regression, kernel ridge regression, ℓ p regression, k-means clustering.
If the sketch matrix S is drawn from a randomized distribution D over matrices R m×n , then the condition in (5.2) can be met with a desired level of probability.
The sketch size given in the definition is optimal [33]. By assuming x is unit norm, the probability in Definition 5.1 can be bounded by the Markov Inequality [54,39]: which leads to a fundamental property: Definition 5.2. (JL-Moment Property Definition 6.1 in [33]) Distribution D over R m×n has (ǫ, δ, ℓ)-JL moment property if for all the unit norm vectors x ∈ R n Having JL-moment property is a stronger condition than being a JL-distribution, but if a distribution D is a (ǫ, e −Ω(ǫ 2 m) )-JL distribution, then it is called as Strong JL Distribution [32] and it automatically verifies the (ǫ, e −Ω(ǫ 2 m) , min(ǫ 2 m, ǫm))-JL moment property [33]. The JL-moment property is important to reach our goal due to the following property called the Approximate Matrix Multiplication (AMM): Theorem 5.3. (AMM Property, Theorem 6.2 in [33]) Given ǫ, δ ∈ (0, 1/2), let D be any distribution over matrices in R m×n with the (ǫ, δ, ℓ)-JL moment property for some ℓ ≥ 2. Then for any A, B ∈ R n×d real matrices, If the sketch matrix satisfies the JL-moment property, setting A = B = U 1 in the AMM property gives the condition in 5.2, but except for the sparse sketch matrices with single non-zero element in each column, the AMM property that is based on the JL-Distribution is too restrictive for the other sketch types mentioned in Corollary 3.2 to satisfy the condition. For them, the ℓ2 subspace embeddings, which is introduced in [50], can be used.
The JL-distributions are special case of (ǫ, δ, k)-OSE with k = n. While the JL distributions projects the whole R n , the OSE distributions embed only a subspace span(U ) ⊆ R n . Consequently, the subspace embeddings provide a more general moment property.
Definition 5.5. (OSE Moment Property [15]) A distribution D over S ∈ R m×n has (ǫ, δ, k, ℓ)-OSE Moment Property, if for all matrices U ∈ R n×k with orthonormal columns, The JL moment property is the special case of the OSE moment property with k = n and the relation between them is given in the following lemma.
By using the OSE Moment property to obtain a generalized AMM property is also possible: Theorem 5.7. ((ǫ, δ, k)-AMM Property [15]) If a distribution D over S ∈ R m×n has the (ǫ, δ, 2k, ℓ)-OSE moment property for some δ < 1/2 and ℓ ≥ 2, then for any A, B: Theorem 5.7 breaks the dependence between the error constant ǫ and the sketch size m in the AMM property given in Theorem 5.3. By Theorem 5.7, the sketch sizes can be chosen in accordance with the embedding size k to satisfy the condition in (5.2) as we will show next. Proof of Corollary 3.2, Item (i). Count Sketch with a single nonzero element in each column and size m ≥ 2/(ǫ ′ 2 δ) has (ǫ ′ , δ, 2)-JL moment property [51]. By Theorem 5.3: Proof of Corollary 3.2, Item (ii). Combining Theorem 4.2 of [14] and Remark 2 of [15] implies that any sketch matrix drawn from an OSNAP [39] with the conditions given in the item (ii) of Corollary 3.2 satisfies the (ǫ ′ , δ, k, log(k/δ))-OSE moment property, thus the (ǫ ′ , δ, k/2)-AMM Property. Setting A = B = U 1 and k = sd λ (A)/2 in Theorem 5.7 gives: with probability at least 1 − δ.

Convergence proof without regularization.
If the LS problem without any regularization is tried to be solved, then we need the following inequality to hold: which can be satisfied with a constant probability, if the sketch matrix is drawn from a distribution D that has the OSE Moment property defined in Definition 5.5. Then, the bound for the sketch size can be derived by using the other results presented above. The details of how to obtain the above inequality can be found in [42].

Proof of Corollary 3.4.
Consider the regularized LS solution with parameter λ and the Truncated SVD solution with parameter ⌈sd λ (A)⌉: where u i 's and v i 's are columns of U and V matrices in the SVD. The Tikhonov regularization with the closed form solution is preferred in practise to avoid the high computational cost of the SVD. The filtering coefficients of the Tikhonov regularization, i +λ , become very close to the binary filtering coefficients of the TSVD, as the decay rate of the singular values of A increases. In these cases, both the solution x * and x † become almost indistinguishable. Thereby, the diagonal matrix, ΣD, which is used in the proof of Corollary 3.2 can be approximated by the diagonal matrix Π where which is equivalent to replace the Tikhonov coefficients by the binary TSVD coefficients. Then, we have where S = SU Π ∈ R m×sd λ (A) has the same distribution as S, since U Π is an orthonormal transformation. By the Marchenko Pastur Law, the minimum and the maximum eigenvalues of this approximation converge to 1 ± sd λ (A) m −2 as m → ∞ and while sd λ (A)/m remains constant [19,20]. The rest of the proof follows from the analysis given in subsection 5.1.
6. The proposed technique to approximately solve the sub-problems. The linear sub-problems in the form of (A T A + λI)x = b, whose solutions are required by all four M-IHS variants, can be approximately solved by using the bidiag2 procedure described in [44]. The bidiag2 procedure produces an upper bidiagonal matrix as following: The upper bidiagonal decomposition R k is computed by using the Lanczos-like three term recurrence: where θ j 's and ρ j 's are chosen so that v j 2 = p j 2 = 1, respectively. Noting that P k and V k are not needed to be orthogonal in AAb Solver, therefore we do not need any reorthogonalization steps. As different from the LSQR, we choose θ 1 v 1 = b with θ 1 = b 2 so that the columns of matrix V k constitute an orthonormal basis for the k-th order Krylov Subspace: Regularization does not affect this property since the Krylov Subspace is invariant under a constant shift. In the k-th iteration of the proposed sub-solver, AAb Solver, let the solution estimate of the linear system be x k = V k y k for some vector y k ∈ R k , i.e., x k ∈ K k (A T A, b), then we have: where R k is obtained by applying a sequence of Givens rotation on R k √ λI k in order to eliminate the sub-diagonal elements coming from the regularization [22]. One instance of this procedure is  where Since R k is an upper bidiagonal matrix, the inverse always exists and φ j can be found analytically as: Furthermore, the solution at the k-th iteration, x k = V k R −1 k f k , can be obtain without computing any inversions by using the forward substitution. Define D k = V k R −1 k : and the relative residual error, to use as stopping criterion, can be found as: The first norm in (i) is zero since the linear system is always consistent. The second term in the second norm is also zero, since b ∈ span(V k ) by the initial choice of θ 1 v 1 = b. By definition, f k = R k y k gives the final results. The overall algorithm is given in Algorithm 6.1. The AAb Solver also is a Krylov Subspace method, therefore, it finds the solution at most in min(n, d, m) iterations in the exact arithmetic, but far fewer number of iterations is sufficient for our purpose.

Conclusions.
We studied the M-IHS techniques, a group of solvers for large scale LS problems, which are obtained by incorporating the Heavy Ball Acceleration into the iterations of the IHS algorithm. We obtained the optimal fixed momentum parameters that maximize the convergence rate by the help of the results in random matrix theory. We examined the effect of ℓ2 norm regularization on the optimal momentum parameters. We showed that the M-IHS variants can be used for any dimension regimes if the statistical dimension is sufficiently smaller than the dimensions of the coefficient matrix and we obtained lower bounds on the sketch size for several randomized distributions in order to get a pre-determined convergence rate with a 8. Appendices: Discussion on the Proposed Error Upper Bound for Iterations of Primal Dual Algorithms in [53]. In this appendix, we provide details of a critical discussion on the steps of the derivation that leads to an error upper bound for the iterations of the primal dual algorithms given in [53]. First, we provide a short list of minor issues that can easily be corrected.
1. During the initialization stage in Line 2 of both Algorithm 4 in page 4097 and Algorithm 5 in page 4098, the residual error vector r (0) must be set to −λy instead of −y, otherwise iterates of the both of the algorithms diverge from the optimal solution. 2. During the initialization stage in Line 2 of Algorithm 7 in 4912, the dual residual error vector r Dual must be set to −λy instead of −y and during the initialization stage of the inner loop iterations in Line 15, the primal residual error vector r (0) P must be set to −R T X T r (t+1) D ; otherwise iterates of the algorithm diverges from the optimal solution. The MATLAB codes provided in the link includes these corrections. In addition to the above mentioned minor issues, there are some major issues as well. Unfortunately, we could not obtain corrective actions on these major issues as we could have done on the minor issues mentioned above. Therefore, a lower bound on the number of inner loop iterations, that guarantee a certain rate of convergence at the main loop, is still an open question for the primal dual algorithms. In the remaining of this appendix, we will provide steps of the derivation presented in [53], along with our critical remarks on their validity.
Consider the following A-IHS updates w t+1 = w t + u t .
We are going to use exactly the same notation as [53] except for that HS subscript for the A-IHS iterates are omitted. In the primal dual algorithms, instead of exact sequence { w t }, a sequence { w t } is obtained due to the approximate minimizers that are used in place of u t . Consequently while the sequence { w t } is obtained after t exact iterations of the A-IHS algorithm, sequence { w t } is obtained after t primal dual iterations in each of which k inner loop updates are used to approximate u t 's. The details of the inner and outer loops can be found in Algorithm 7 of [53]. The aim of the Theorem 9 is to establish an upper bound for w t+1 − w *