{'fraction_non_alphanumeric': 0.048260266228092714, 'fraction_numerical': 0.020815130801254047, 'mean_word_length': 4.010299987125016, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 42, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On pseudoinverse-free randomized methods for linear systems: Unified framework and acceleration

We present a new framework for the analysis and design of randomized algorithms for solving various types of linear systems, including consistent or inconsistent, full rank or rank-deficient. Our method is formulated with four randomized sampling parameters, which allows the method to cover many existing randomization algorithms within a unified framework, including the doubly stochastic Gauss-Seidel, randomized Kaczmarz method, randomized coordinate descent method, and Gaussian Kaczmarz method. Compared with the projection-based block algorithms where a pseudoinverse for solving a least-squares problem is utilized at each iteration, our design is pseudoinverse-free. Furthermore, the flexibility of the new approach also enables the design of a number of new methods as special cases. Polyak's heavy ball momentum technique is also introduced in our framework for improving the convergence behavior of the method. We prove the global linear convergence rates of our method as well as an accelerated linear rate for the case of the norm of expected iterates. Finally, numerical experiments are provided to confirm our results.


Introduction
Randomized iterative methods such as the randomized Kaczmarz method have been very popular recently for solving large-scale linear systems. They are preferred mainly because of their low iteration cost and low memory storage so that they scale better with the size of the problems. Randomized methods are now playing a major role in areas like numerical linear algebra, scientific computing, and optimization. Their applications include tensor recovery [12], image reconstruction [26], signal processing [8], optimal control [47], partial differential equations [46], and machine learning [11].
Let us consider the following linear system where A ∈ R m×n and b ∈ R m . The Kaczmarz method [28], also known as algebraic reconstruction technique (ART) [20,26], is one of the popular methods for solving (1). At each step, the method cyclically projects the current iterate onto the solution space of a single constraint and converges to a certain solution of the consistent linear systems.
However, the rate of convergence is hard to obtain. In a seminal paper [55], Strohmer and Vershynin studied the randomized Kaczmarz (RK) method and proved its linear convergence for solving consistent linear systems. Their elegant result has inspired a lot of subsequent work. We refer the interested reader to [3-5, 13-16, 21, 22, 25, 27, 32, 45, 52-55, 57] and the references therein.
In [22], Gower and Richtárik developed a versatile randomized iterative method for solving consistent linear systems which includes the RK algorithm as a special case. The method updates with the following iterative strategy: where S ∈ R m×q are random matrices drawn from a distribution D at each step, B ∈ R n×n is a user-given symmetric positive matrix, α ∈ (0, 2) is a stepsize parameter, and the symbols † and ⊤ denote the Moore-Penrose pseudoinverse and the transpose of either a vector or a matrix, respectively. By varying these two parameters S and B, one can recover a number of popular methods as special cases, including the RK method, the Gaussian Kaczmarz method, the randomized Newton method, random Gaussian pursuit, and variants of all these methods using blocks and importance sampling. Their method is the first algorithm uncovering the close relations between these methods and is also known as the sketch-andproject method. The method is shown to enjoy linear convergence from which existing convergence results can also be obtained.
However, the main disadvantage of (2) is that each iteration is expensive since we need to apply the pseudoinverse to a vector, or equivalently, we have to solve a least-squares problem at each iteration. Moreover, it is also difficult to parallelize. In this paper, we aim to develop randomized algorithms with accelerated convergence rates that do not employ the pseudoinverse. Indeed, we will utilize Polyak's heavy ball momentum [49] technique to accelerate the convergence rate of the method.
The heavy ball method. Consider the following optimization problem: min x∈R n f (x), where f is a differentiable convex function. Gradient descent (GD) is widely used method for solving the above problem. Starting with an arbitrary point x 0 , the iteration scheme of the original GD method reads as where γ k > 0 is a stepsize and ∇f x k denotes the gradient of f at x k . For convex functions with L-Lipschitz gradient, GD needs O(L/ε) steps to guarantee the error of the solution is within ε. When, in addition, f is µ-strongly convex, it converges linearly with convergence rate of O(log(ε −1 )(L/µ)) [43]. To improve the rate of convergence, Polyak proposed to modify GD by the introduction of a (heavy ball) momentum term, ω x k − x k−1 . This leads to the gradient descent method with momentum (mGD), popularly known as the heavy ball method: Polyak [49] proved that for twice continuously differentiable objective function f (x) with µ strong convexity constant and L-Lipschitz gradient, mGD achieves a local accelerated linear convergence rate of O log(ε −1 ) L/µ (with the appropriate choice of the stepsize parameters γ k and momentum parameter ω).
1.1. Our contribution. In this paper, we propose the following more general randomized algorithms framework for solving the systems of linear equations (3) x where α > 0 is a stepsize parameter, ω ≥ 0 is a momentum parameter , and (S 1 , S 2 , T 1 , T 2 ) is a random parameter matrix pair which is sampled independently in each iteration from a certain distribution D. The main contributions of this work are as follows.
1. We propose a novel randomized algorithmic framework for solving different types of linear systems, including consistent or inconsistent, full rank or rank-deficient.
We show the method enjoys a global linear convergence rate. We also study the convergence of the quantity E[x k − x 0 * ] 2 2 . In this case, we show that by a proper combination of the stepsize parameter α and the momentum parameter ω, the proposed method enjoys an accelerated convergence rate.
2. The convergence of the algorithm in some special cases is also studied. Particularly, global linear convergence rates are established under various formulations and assumptions about the sampling matrices and the linear systems.
3. We recover several known methods as a special case of our general framework which enables us to draw previously unknown links between those methods. The flexibility of the method also allows us to adjust the parameter matrices to obtain completely new methods. Particularly, we study the methods based on Gaussian sampling.

Related work.
This result is related to the recent work of Du and Sun in [15], where two novel pseudoinverse-free randomized methods for solving a linear system are proposed.
In particular, they studied the framework in (3) with S 1 = S 2 , T 1 = T 2 , and ω = 0. When A is a full column rank matrix, they [15,Theorem 4] showed that the method converges linearly in the mean square sense. We emphasize that our convergence result is more general and works for all types of matrices (consistent or inconsistent, full rank or rank-deficient); see Theorems 3.1 and 4.1. In addition, we use Polyak's heavy ball momentum for accelerating the convergence of the method in our framework.
There is a variant of the randomized block Kaczmarz (RBK) that avoids the issues of pseudoinverse, we emphasize the results of Necoara [42]. This leads to the following iteration: where the weights ω i ∈ [0, 1] such that i∈J ω i = 1, J is drawn from a certain distribution D at each step, α > 0 is the stepsize parameter, and a ⊤ 1 , . . . , a ⊤ m denote the rows of A. Let I :J denote a column concatenation of the columns of the m × m identity matrix I indexed by J, and the diagonal matrix D J := diag( √ ω i / a i 2 , i ∈ J). Then the iteration (4) can alternatively be written as where S = I :J D J , which can be viewed as randomized matrices selected from the distribution D at each step. We note that the iteration scheme (5) can also be recovered by our framework in the same special cases; see Section 5.4 for more details. Indeed, it is also an interesting topic that one can equip (5) with more general sampling matrices S for obtaining new classes of randomized algorithms.
Recently, Loizou and Richtárik [37] studied several kinds of stochastic optimization algorithms enriched by Polyak's heavy ball momentum for solving convex quadratic problems.
They proved the global, non-asymptotic linear convergence rates of the proposed methods.
The momentum variants of the sketch-and-project method has been investigated in [38,51].
In [39], Polyak's momentum technique was incorporated into the stochastic steepest descent methods for solving linear systems. Their momentum framework integrated several momentum algorithms, such as the Kaczmarz method, the block Kaczmarz method, and the coordinate descent, into one framework. In [13], the authors studied the Douglas-Rachford (DR) method [1,2,10,17,24,33,34] enriched with randomization and Polyak's heavy ball momentum for solving linear systems. They showed the power of randomization in simplifying the analysis of the DR method and made the divergent r-sets-DR method converge linearly.
We note that a more widely used, and better understood alternative to Polyak's momentum is the momentum introduced by Nesterov [43,44], leading to the famous accelerated gradient descent (AGD) method [6]. Recently, variants of Nesterov's momentum have also been introduced for the acceleration of stochastic optimization algorithms [30]. In [35], Liu  1.3. Organization. The remainder of the paper is organized as follows. In Section 2, we give some notations and a useful lemma. In Section 3, we present the pseudoinversefree randomized (PFR) method and analysis its convergence properties. In Section 4, we propose the PFR with momentum (mPFR) and show its accelerated linear convergence rate. We also study the convergence of the algorithm in some special cases. In Section 5, we mention how by selecting the parameters (T 1 , T 2 , S 1 , S 2 ) of our method, we recover several existing methods as well as obtain new methods. In Section 6, we perform some numerical experiments to show the effectiveness of the proposed methods. Finally, we conclude the paper in Section 7. Proofs of all main results can be found in the "Appendix". For any vector x ∈ R n , we use x i , x ⊤ and x 2 to denote the i-th entry, the transpose and the Euclidean norm of x, respectively. We use e i to denote the column vector with a 1 at the i-th position and zeros elsewhere. For any matrix A, we use a i , A j , a i,j , A ⊤ , A † , A 2 , A F and Range(A) to denote the i-th row, the j-th column, the (i, j)-th entry, the transpose, the Moore-Penrose pseudoinverse, the spectral norm, the Frobenius norm and the column space of A, respectively. Let A ∈ R m×n , we use

Notation and preliminaries
where r is the rank of A and we use σ min (A) to denote the smallest nonzero singular values of A. We see that For a square matrix B ∈ R n×n , we use Tr(B) := Throughout this paper, we use x * to denote a certain solution of the linear system (1), and for any x 0 ∈ R n , we set We know that x * 0 is the orthogonal projection of x 0 onto the set and x * LS is the least-squares or the least-norm least-squares solution of the linear system. In this paper, by "momentum" we refer to the heavy ball technique originally developed by Polyak [49] to accelerate the convergence rate of gradient-type methods.
Throughout this paper, we will always assume that the sampling matrices satisfy the following assumption.
Assumption 2.1. The random parameter matrix pair (S 1 , S 2 , T 1 , T 2 ) is sampled independently in each iteration from a distribution D and satisfies We note that in Assumption 2.1, the random parameter matrices S 1 , S 2 ∈ R m×p and T 1 , T 2 ∈ R n×s can be independent or not. We also emphasize that we do not restrict the numbers of columns of S 1 , S 2 , T 1 and T 2 ; indeed, we allow p and s to vary (and hence p and s are random variables).
We will utilize the following result to estimate the upper bound of the stepsize parameter for the Gaussian-type randomized algorithms.
Lemma 2.1. Suppose that S ∈ R m×p is a Gaussian matrix, i.e. has independent N (0, 1)

The pseudoinverse-free randomized method for linear systems
Let us first study the pseudoinverse-free randomized (PFR) method for solving linear systems. Given an arbitrary initial guess x 0 ∈ R n , the k-th iterate of the PFR method is defined as where the stepsize parameter α > 0. Note that (6) is equivalent with the no-momentum variant of (3), i.e. ω = 0, and under Assumption 2.1, we have which is the update of the Landweber iteration [31]. The PFR method is formally described in Algorithm 1.
3.1. Convergence of iterates: Linear rate. In this subsection, we analyze some convergence properties of Algorithm 1. We emphasize that the matrix pair (S 1 , S 2 , T 1 , T 2 ) selected from the distribution D satisfies the Assumption 2.1. The following theorem shows that Algorithm 1 converges linearly in expectation for different types of linear systems (consistent or inconsistent, full rank or rank-deficient).

Algorithm 1 The pseudoinverse-free randomized (PFR) algorithm
Input: Distribution D from which to sample matrices, A ∈ R m×n , b ∈ R m , stepsize/relaxation parameter α > 0, k = 0, and an initial x 0 ∈ R n . 1: Randomly select a matrix pair (S 1 , S 2 , T 1 , T 2 ) from the distribution D.
. 3: If the stopping rule is satisfied, stop and go to output. Otherwise, set k = k + 1 and return to Step 1. Output: The approximate solution.
Theorem 3.1. Let x * LS be the least-squares or the least-norm least-squares solution of the linear system Ax = b and {x k } ∞ k=0 be the sequence generated by Algorithm 1. Denote by , then the k-th Next, we consider the convergence of x k − x 0 * . Particularly, we will consider the conver- Theorem 3.2. Let x 0 ∈ R n be an arbitrary initial vector and x 0 3.2. Convergence direction. Inspired by the recent work of Steinerberger [53], this section aims to consider the convergence direction of Algorithm 1. The following result shows different convergence rates of Algorithm 1 along different singular vectors of A.
Theorem 3.3. Let x 0 ∈ R n be an arbitrary initial vector and x 0 This exhibits that Algorithm 1 decays exponentially at different rates depending on the singular values and accounts for the typical semiconvergence phenomenon. That is, the residual Ax k − b 2 2 decays faster at the beginning, but then gradually stagnates. Recently, such semiconvergence phenomenon has been studied in the literature for randomized methods. In [53], Steinerberger studied semiconvergence phenomenon for the randomized Kacazmarz method. In [56], Wu and Xiang exploited the semiconvergence phenomenon for the sketch-and-project method [22], where they generalized the study in [27] and split the total error into the low-and high-frequency solution spaces. Very recently, the semiconvergence phenomenon has also been studied for randomized DR method in [13].

Momentum acceleration
In this section, we provide the momentum induced PFR method, i.e. the mPFR method for solving a linear system. The method is formally described in Algorithm 2.

Algorithm 2
The pseudoinverse-free randomized algorithm with momentum (mPFR) Input: Distribution D from which to sample matrices, A ∈ R m×n , b ∈ R m , stepsize/relaxation parameter α > 0, the heavy ball/momentum parameter ω, k = 1, and initial vectors x 1 , x 0 ∈ R n . 1: Randomly select a matrix pair (S 1 , S 2 , T 1 , T 2 ) from the distribution D.
2: Update which can be viewed as the momentum variant of the Landweber iteration. In the rest of this section, we will study the convergence properties of Algorithm 2. We note that the matrix pair (S 1 , S 2 , T 1 , T 2 ) selected from the distribution D satisfies Assumption 2. Theorem 4.1. Let x * LS be the least-squares or the least-norm least-squares solution of the linear system Ax = b and the initial vectors , ω ≥ 0 and that the expressions Then the iteration sequence of residuals {r k } ∞ k=1 generated by Algorithm 2 satisfies if ω = 0, and τ = q − γ 1 ≥ 0.
Let us explain how to choose the parameters α and ω such that γ 1 + γ 2 < 1 is satisfied in Theorem 4.1. Indeed, let , then τ 2 > 0 and the condition γ 1 + γ 2 < 1 now is satisfied for Next, we compare the convergence rates obtained in Theorems 3.2 and 4.1. From the definition of γ 1 and γ 2 , we know that convergence rate q(ω) in Theorem 4.1 can be viewed as a function of ω. Note that since ω ≥ 0, we have Clearly, the lower bound on q is an increasing function of ω. Also, for any ω the rate is always inferior to that of Algorithm 1 in Theorem 3.1.

4.2.
Accelerated linear rate for expected iterates. In this subsection, we study the convergence of the quantity E x k − x 0 * 2 2 to zero for Algorithm 2. We show that by a proper combination of the relaxation parameter α and the momentum parameter ω, Algorithm 2 enjoys an accelerated linear convergence rate in mean.

Theorem 4.2.
Suppose that x 1 = x 0 ∈ R n are arbitrary initial vectors and x 0 Then there exists a constant c > 0 such that for all k ≥ 0 we have In contrast, based on Theorem 4.2 we have, for which is a quadratic improvement on the above result.

Convergence analysis: Special cases.
In this subsection, we study the convergence results of Algorithm 2 (Algorithm 1) for some special cases where the coefficient matrix A is full column rank or the sampling matrices T 1 , T 2 , S 1 , S 2 are scalar matrices.
The full column rank case. If the matrix A is full column rank, then we can study the convergence rate of the quantity E x k − A † b 2 2 . We show that for a range of stepsize parameters α > 0 and momentum terms ω ≥ 0, the method enjoys a global linear convergence rate.
Theorem 4.4. Suppose that A has full column rank. Let x 1 = x 0 ∈ R n be the initial vectors and denote , ω ≥ 0 and that the expressions Furthermore, if the linear systems are consistent, then we can show that the quantity E x k − A † b 2 2 converges to zero linearly. Besides, one can also choose a larger stepsize.
Theorem 4.5. Consider a consistent system Ax = b with A being full column rank. Let and that the expressions The scalar matrices case. If the sampling matrices S 1 , S 2 or T 1 , T 2 are scalar matrices, then we can study the convergence rate of the quantity E r k − r * 2 2 or E x k − A † b 2 2 for Algorithm 2 (Algorithm 1).
Theorem 4.6. Suppose that the sampling matrix S 1 = S 2 = √ ρI and ρ > 0. Let x * LS be the least-squares or the least-norm least-squares solution of the linear system Ax = b and the initial vectors , ω ≥ 0 and that the expressions if ω = 0, and τ = q − γ 1 ≥ 0.
If T 1 , T 2 are scalar matrices, then we have the following two results for Algorithm 2 (Algorithm 1).
Theorem 4.7. Suppose that the linear system Ax = b is consistent and the sampling matrix , ω ≥ 0 and that the expressions Remark 4.8. If we choose the initial vectors k=0 generated by Algorithm 2 now converges to the least-squares or the least-norm least-squares solution This implies that the quantity E , and hence harder to push to zero. As a corollary, the convergence rate of E x − x * Finally, we present the following result for the case where Theorem 4.10. Let x * LS be the least-squares or the least-norm least-squares solution of the linear system Ax = b and the initial vectors x 1 = x 0 ∈ R n . Denote by r * = Ax * LS − b and r k = Ax k − b. Suppose that for any random parameter matrix pair (S 1 , S 2 , T 1 , T 2 ) sampled from the distribution D, it satisfies that , ω ≥ 0 and that the expressions Table 1. Special cases of Algorithm 2, where S is an m × p Gaussian matrix, T is an n × s Gaussian matrix, and η is a Gaussian vector. mSGC works for the case where A is a symmetric matrix.

Special cases and new efficient methods
Our framework flexibility allows us to adjust the parameter matrices and leads to a number of popular methods. In this section, we briefly mention how by selecting the parameters (S 1 , S 2 , T 1 , T 2 ) of our method, we recover several existing methods and lead to completely new methods; see Table 1 for a quick summary.

Randomized
Kaczmarz method with momentum. The randomized Kaczmarz method is one special case of Algorithm 1. Choosing and Algorithm 1 recovers the RK iteration For this case, we have β 3 = 1 and Theorem 4.7 yields the convergence estimate of [9, Theorem 1]: where the stepsize parameter α ∈ (0, 2). Furthermore, Algorithm 2 recovers the RK with momentum (mRK) Theorem 4.7 also yields the following convergence result of mRK for solving the consistent linear systems: , and the stepsize parameter 0 < α < 2. We note that when The above result has first been established in [37,Theorem 1]. However, they only treated the consistent linear systems. In this paper, our convergence results are more general, see Theorems 4.1 and 4.2, which work for various types of linear systems. We note that the method studied in [37] is in the setting of quadratic objectives, where the equivalence between stochastic gradient descent, stochastic Newton method, and stochastic proximal point method can be established.

5.2.
Randomized Gauss-Seidel method with momentum. The randomized Gauss-Seidel (RGS) method, also known as the randomized coordinate descent (RCD) method, is one special case of Algorithm 1. Choosing S 1 = S 2 = I and and Algorithm 1 recovers the RGS/RCD iteration For this case, we have β 2 = 1 and Theorem 4.6 yields the convergence estimate of [23]: where the stepsize parameter α ∈ (0, 2). In addition, Algorithm 2 derives the following momentum variant of RGS/RCD (mRGS/mRCD): Theorem 4.6 yields the following convergence result for mRGS/mRCD: , and 0 < α < 2. We note that for solving the consistent linear system, the mRCD has also been studied in [37].

Now we have
and By Algorithm 2, we can obtain the DSGS with momentum (mDSGS) Let us consider the convergence properties of mDSGS. If A is full column rank and the linear system is consistent, then Theorem 4.5 yields the following global convergence rate for mDSGS: , and the stepsize parameter 0 < α < 2 n . When ω = 0, the above convergence result recovers the result obtained in [50, Theorem 1] for DSGS.
If A is not full column rank and the linear system is consistent, then Theorem 4.10 yields the following global convergence rate for mDSGS: , and 0 < α < . To the best of our knowledge, mDSGS has never been analyzed before in any setting.

Randomized block Kaczmarz with momentum.
Our framework also extends to block formulations of the RK method studied in [15]. However, to the best of our knowledge, no momentum variants of this method were analyzed before.
and if p = 1, then By Algorithm 1, the randomized block Kaczmarz (RBK) method updates with the following iterative strategy: Assume that a i 2 = 1, i = 1, . . . , m, the RBK algorithm (11) can be viewed as a special case of the randomized algorithm (4) proposed in [42]. Indeed, let J be the set consisting of the uniform sampling of p different numbers of [m], ω i = 1 p and choose the parameter α appropriately, (4) recovers (11). By Algorithm 2, the momentum variant of RBK (mRBK) updates with the following iterative strategy: Theorem 4.7 yields the following global convergence rate of mRBK for solving consistent linear systems: and β 3 is given by (9) and (10).

5.5.
Randomized block coordinate descent with momentum. Our framework extends to block formulations of the randomized coordinate descent method studied in [15].
To the best of our knowledge, this is the first time that the momentum variant of this method is considered.
Assume 1 ≤ s ≤ n. Let L denote the set consisting of the uniform sampling of r different numbers of [n]. Choose T = n/sI :L to be a column concatenation of the columns of the n × n identity matrix I indexed by L. Let T 1 = T 2 = T and S 1 = S 2 = I A F . Then we know that and if s = 1, then By Algorithm 1, the randomized block coordinate descent (RBCD) method updates with the following iterative strategy: By Algorithm 2, the momentum variant of RBCD (mRBCD) updates with the following iterative strategy: Theorem 4.6 shows that mRBCD converges linearly in expectation for solving different types of linear systems if ω = 0, and τ = q − γ 1 ≥ 0 and β 2 is given by (12) and (13). 5.6. Block Gaussian Kaczmarz with momentum. In this subsection, we study the method based on block Gaussian row sampling and derive a new randomized block method.
Let us first propose one special case of Algorithm 1 by using block Gaussian row sampling and refer to it as the block Gaussian Kaczmarz (BGK) method. Choose S to be an m × p random matrix whose entries are independent, mean zero, Gaussian random variables, and where the second equality follows from Lemma 2.1. It follows from Algorithm 1 that the BGK method constructs x k+1 by If we take p = 1, then we can get the Gaussian Kaczmarz proposed in [22]. By Algorithm 2, we can also obtain the accelerated BGK with momentum (mBGK) Theorem 4.7 yields the following global convergence rate of mBGK for solving consistent linear systems:

5.7.
Block Gaussian least-squares with momentum. In this subsection, we study the method based on block Gaussian column sampling. Choose T to be an n × s random matrix whose entries are independent, mean zero, Gaussian random variables, and set T 1 = T 2 = T .
. Then we know that where the second equality follows from Lemma 2.1. Algorithm 1 then has the following form , which we call the block Gaussian least-squares (BGLS) method. By Algorithm 2, we can also obtain the accelerated BGLS with momentum (mBGLS) By Theorem 4.6, we know that mBGLS yields a global convergence rate for solving various types of linear systems

The symmetric matrix cases.
Recall that all of the special cases discussed above focus exclusively on the cases T 1 = T 2 and S 1 = S 2 . In this subsection, we show that one can choose T 1 = T 2 and S 1 = S 2 .
We consider the case where A ∈ R n×n is a symmetric matrix with Tr(A) = 0. Choosing S 1 = T 2 = η with η being a random vector η ∽ N (0, I n ), and (T 1 , S 2 ) = e i a i,j , e j Tr(A) , and the index pair (i, j) is randomly selected with probability (Tr(A)) 2 . By Algorithm 2, we get the following algorithm: which we call the symmetric Gaussian coordinate with momentum (mSGC). Theorem 4.1 yields the following global convergence rate for mSGC: if ω = 0, and τ = q − γ 1 ≥ 0 with γ 1 = 1 + 3ω + 2ω 2 − (2α + αω)

Numerical experiments
In this section, we study the computational behavior of the proposed algorithms. In particular, we focus mostly on the evaluation of the performance of mRBK, mRBCD, mBGK, and mBGLS. We divide the algorithms discussed in Section 5 (see Table 1) into two groups for comparison: (1) mRK, mDGSG, mRBK, and mBGK for solving consistent linear systems; (2) mRGS, mRBCD, and mBGLS for solving inconsistent linear systems.
All the methods are implemented in Matlab R2022a for Windows 10 on a desktop PC with Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz and 16 GB memory.
6.1. Numerical setup. In our test, we use the relative solution error (RSE) or the relative residual error (RRE) as the stopping criterion. We use κ 2 (A) = A 2 A † 2 to denote the 2-norm condition number.
During our experiment, we set the stepsize parameter α = 1 for both mRK and mRGS.
For mDSGS, we set α = 1 n if the matrix A is full-rank, and set α = for mRBK, where β 3 is given by (9) and (10). For mBGK, we set where β 2 is given by (12) and (13). For mBGLS, we set α = . We note that the choices of those stepsize parameters α are not arbitrary. From the discussions in Section 5, we know that the best convergence rates of those algorithms with no-momentum are obtained precisely for those choices of α.
Therefore, the comparisons in our experiment will be with the best-in-theory no-momentum variants.
6.2. The effect of heavy ball momentum. In this subsection, we study the computational behavior of momentum variants of the proposed methods and compare them with their no momentum variants for both synthetic and real data.
Synthetic data. The synthetic data is generated based on the Matlab function randn.
Specifically, we generate the matrix A by using the Matlab function randn(m,n) and generate the exact solution x * by randn(n,1). The consistent system is constructed by setting b = Ax * . For the inconsistent system, we first generate a △b = null(A ′ ) * randn(m − r, 1), We also test our algorithms on problem instances generated with different given values of condition number κ 2 (A). To do so, we first use the Matlab function randn(m,n) to generate a Gaussian matrix. Then, we use SVD to modify the singular values to obtain a matrix with the desired values of κ 2 (A). This is done by linearly scaling the differences of all singular values with σ 2 min (A).
Real data. The real-world data are available via the SuiteSparse Matrix Collection [29].
Each dataset consists of a matrix A ∈ R m×n and a vector b ∈ R m . In our experiments, we only use the matrices A of the datasets and ignore the vector b. As before, to ensure the consistency of the linear system, we first generate the solution by x * = randn(n, 1) and then set b = Ax * . For the inconsistent system, we set b = Ax * + null(A ′ ) * randn(m − r, 1).
In our test, we use x 0 = 0 ∈ R n or x 1 = x 0 = 0 ∈ R n as an initial point. For each consistent linear system, we run mRK, mDGSG, mRBK, and mBGK (Figures 1 and 2) for several values of momentum parameters ω and plot the performance of the methods (average after 10 trials) for RSE. For each inconsistent linear system, similarly, we run mRGS, mRBCD, and mBGLS (Figures 3 and 4) for several values of momentum parameters ω and plot the performance of the methods (average after 10 trials) for RRE. We set p = 20 and s = 20 for mRBK, mBGK, mRBCD, and mBGLS, as they are always good options for sufficient fast convergence of those methods.
For ω = 0, the methods now are equivalent with their no-momentum variants. We note that in all of the presented tests the momentum parameters ω of the methods are chosen to be nonnegative constants that do not depend on parameters that are not known to the users such as σ 2 min (A) and σ 2 max (A). From Figures 1, 2, 3, and 4, it is clear that the introducing of momentum term leads to an improvement in the performance of the methods. More specifically, from the figures we observe the following: (1) The momentum technique can improve the convergence behavior of the methods.
(2) For consistent linear systems with different values of κ 2 (A) (well-or ill-conditioned), and for different types of data (synthetic or real-world), ω = 0.4 is always a good option for a sufficient fast convergence of both mRK and mDSGS; ω = 0.4 or 0.5 is always a good option for both mRBK and mBGK (Figures 1 and 2).   Figure 2. Performance of mRK, mRBK, mBGK, and mDSGS with different momentum parameters ω for consistent linear systems. The coefficient matrices are from SuiteSparse Matrix Collection [29]. All plots are averaged over 10 trials. The title of each plot indicates the test algorithm, the dimension of the matrix A, and the value of condition number κ 2 (A). We set p = 20 for both mRBK and mBGK.
6.3. Comparison among the proposed methods. In this subsection, we will compare mRK, mRBK, and mBGK for solving the consistent linear systems and compare mRGS, First we compare mRK, mRBK, and mBGK on synthetic linear systems generated with the Matlab functions randn and sprandn; see Figure 5. The poor performance of mBGK on the dense problem generated using randn can be found in Figure 5 (a). This is because that the Gaussian method almost always requires more flops to reach a solution as it requires the expensive matrix-vector product. We can observe from Figure 5 (a) that mRBK performs better than mRK and mBGK on the dense problem. In Figure 5 (b), we compare the methods on a sparse linear system generated using the Matlab sparse random matrix function sprandn(m,n,density,rc), where density is the percentage of nonzero entries and rc is the reciprocal of the condition number. On this sparse problem, mBGK is more efficient than mRK and mRBK.  Figure 4. Performance of mRGS, mRBCD, and mBGLS with different momentum parameters ω for inconsistent linear systems. The coefficient matrices are from SuiteSparse Matrix Collection [29]. All plots are averaged over 10 trials. The title of each plot indicates the test algorithm, the dimension of the matrix A, and the value of condition number κ 2 (A). We set s = 20 for both mRBCD and mBGLS.
In Figure 6, we compare mRGS, mRBCD, and mBGLS on inconsistent linear system. The high iteration cost of mBGLS resulted in poor performance on the dense problem generated using rand can be found in Figure 6 (a). On this dense problem, mRBCD is more efficient.
From Figure 6 (b), we can see that mBGLS performs better than mRGS and mRBCD for solving sparse problems. We note that since Matlab performs automatic multithreading when calculating matrix-vector products, which was the bottleneck cost in the Gaussian sampling based methods. Hence despite the higher iteration cost of the Gaussian sampling based methods, their performance, in terms of the wall-clock time, is comparable to that of other methods when the coefficient matrix is sparse.
6.4. Faster method for average consensus. Average consensus is a central problem in distributed computing and multi-agent systems [7,38]. It is raised in many applied areas, such as PageRank, coordination of autonomous agents, and rumor spreading in social networks. In the average consensus (AC) problem, we are given an undirected connected That is, the exchange of information can only occur between connected nodes (neighbours).
Recently, under an appropriate setting, the famous randomized pairwise gossip algorithm [7] for solving the AC problem has been proved to be equivalent to the RK method. One may refer to [36,38] for more details.
In our test, the linear system is the homogeneous linear system (Ax = 0) with matrix A ∈ R m×n being the incidence matrix of the undirected graph. The initial values of the nodes are chosen arbitrarily, i.e. using the Matlab function rand(n,1), and the algorithms will find the average of these values. We note that the incidence matrix A is rank deficient.
More specifically, it can be shown that rank(A) = n − 1 [36,38]. In comparing the methods, the initial vector is chosen as x 0 = c or x 1 = x 0 = c. The computations are terminated once RSE is less than 10 −12 .
We use three popular graph topologies in the literature on wireless sensor networks: the line graph, the cycle graph, and the 2-dimensional random geometric graph G(n, p).
In practice, G(n, p) is considered ideal to model wireless sensor networks because of their particular formulation. In the experiments the 2-dimensional G(n, p) is formed by placing n nodes uniformly at random in a unit square where there are only edges between nodes whose Euclidean distance is within the given radius p. To preserve the connectivity of G(n, p), a radius p = p(n) = log(n)/n is recommended [48]. Table 2 summarizes the results of the experiment, where we use mRK, mRBK, and mBGK for solving the consistent linear system Ax = 0. It is clear that the introduction of the momentum term improves the performance of the methods. It is shown that mRBK (ω = 0.5) is more efficient than other methods. Table 2. Performance of mRK, mRBK, and mBGK for the AC problem. We set p = 20 for both mRBK and mBGK. All results are averaged over 10 trials. The label "-" in the table means that the algorithms fail to find the solution within 200 seconds.

Appendix. Proof of the main results
The following lemma is crucial in our proof.
Proof. We first show that τ ≥ 0 and γ 2 = (γ 1 + τ ) τ . Indeed, since the nonnegativity of τ can be derived from γ 2 ≥ 0. It is easy to verify that τ satisfies By (15) and adding τ F k to both sides of (14), we get where the inequality follows from τ ≥ 0.
Noting that and A ⊤ r * = 0, we have Since r k − r * ∈ Range(A), we know that Using the fact that a − b 2 2 ≤ 2 a 2 2 + 2 b 2 2 , we have Note that and similarly, Combining (16), (17), (18), (19) and (20) yields Taking the expectation over the entire history, we have Now we will prove Theorem 4.1.
Proof of Theorem 4.1. Straightforward calculations yield We will now analyze the three expressions a ○, b ○, c ○ separately. The first expression can be written as We will now bound the second expression. First, we have The third expression can be bounded as By substituting all the bounds into (21), we obtain
We first prove Theorem 4.2.
Proof of Theorem 4.2. Set s k := V ⊤ E[x k − x 0 * ]. Then from Lemma 8.2 we have which can be rewritten in a coordinate-by-coordinate form as follows: where s k i indicated the i-th coordinate of s k .
If σ i (A) = 0, then (28) takes the form: If σ i (A) > 0. We use Lemma 8.3 to establish the desired bound. Since 0 < α ≤ < 0, where the last inequality can be shown to hold for Using Lemma 8.3, the following bound can be deduced where p i is a constant depending on the initial conditions (we can simply choose p i = |c 0 | + |c 1 | .) Now put the two cases together, for all k ≥ 0 we have Hence, we have This completes the proof of the theorem.
Proof of Theorem 3.3. Using similar arguments as that in the proof of Theorem 4.2 and (28), we know that and A has full column rank, we have (35) 2α Combining (30), (31), (32), (33), (34), and (35) yields Then using Lemma 8.3 and the same arguments as that in the proof of Theorem 4.1, we can get the theorem.
Proof of Theorem 4.10. Since T 1 T ⊤ 2 A ⊤ S ⊤ 1 S 2 r * = 0, then (36) becomes (r k − r * ) ≤ β r k − r * 2 2 . By using similar arguments as that in the proof of Theorem 4.6, we know the theorem holds.
Proof of Theorem 4.7. Similar to (30) and note that Ax 0 * = b, we have Since T 1 = T 2 = √ ρI, we have We also have . which together with (37) and (38) implied that Note that the iteration of Algorithm 2 now becomes and as we choose x 1 = x 0 , so x k ∈ Range(A ⊤ ) + x 0 for k = 0, 1, . . ., and therefore, Note that 2α − α 2 ρ 2 β 3 A 2 F ≥ 0, from (39), we have Then using Lemma 8.3 with ζ = 0 and the same arguments as that in the proof of Theorem 4.1, we can obtain the desired result.
8.5. Proof of Lemma 2.1. The following lemma is needed for our proof. Set ζ := U ⊤ η, then we have This implies that Hence we have where Q i denotes the i-th column of the matrix Q. Let denote the (ℓ, i)-th entry of Q as q ℓi , then q 1i q 1j q 1i q 2j . . . q 1i q mj q 2i q 1j q 2i q 2j . . . q 2i q mj . . . . . . . . . . . .