linear system; randomized PLSS; expected linear convergence; Kaczmarz method; coordinate descent

1 Introduction↩︎

Consider a large, possibly consistent or not, general linear system \[\label{1-1} \mathbf{A}\mathbf{x}=\mathbf{b},\tag{1}\] where \(\mathbf{A}\in \mathbb{R}^{m\times n}\) and \(\mathbf{b}\in \mathbb{R}^{m}\) are given, and \(\mathbf{x}\in \mathbb{R}^{n}\) is unknown. Such a linear system underpins widespread applications across image reconstruction [1], [2], signal processing [3], [4], and partial differential equations [5], etc. While the regular iterative solvers, such as Krylov subspace methods and Hermitian and skew-Hermitian methods [6]–[10], for 1 have been extensively matured, randomized solvers [11]–[21] are growing in popularity due to a revolutionary advantage. Specifically, the traditional iterative solvers, demanding full-matrix access per iteration, are costly to generate highly accurate solutions, but they are not needed or even desired since Eq.@eq:1-1 may constitute a random subset of a larger dataset, or it may be contaminated by noise. Conversely, randomized solvers address only a sequence of small projected systems, rather than 1 , while maintaining competitive convergence to yield the desired solutions. This makes them powerful for tackling large-scale, memory-constrained, or data-streaming applications. In this paper, we propose a family of randomized projection methods that solve a sequence of smaller systems and offer extraordinary advantages over state-of-the-art randomized solvers in terms of accuracy and convergence.

1.1 Notations↩︎

Let \(\mathbb{R}^{m\times n}\) be the collection of all real matrices with size \(m\times n\). Denote by \(\mathbf{A}^{(i)}\) and \(\mathbf{A}_{(j)}\), the \(i\)-th row and \(j\)-th column of a matrix \(\mathbf{A}\in \mathbb{R}^{m\times n}\), respectively. For index sets \(\tau_k\) and \(\tilde{{\tau}}_k\), \(\mathbf{A}_{\tau_k}\) and \(\mathbf{A}_{:,\tilde{{\tau}}_k}\) represent the row and column submatrices of \(\mathbf{A}\), respectively. The symbol \(\mathbf{e}_{j}\) represents the \(j\)-th column of the identity matrix \(\mathbf{I}\), with dimension depending on the context. The Frobenius norm of \(\mathbf{A}\) is defined as \(\|\mathbf{A}\|_F:= (\sum\limits_{j = 1}^n {\sum\limits_{i = 1}^m {{{\left| {{\mathbf{A}_{ij}}} \right|}^2}} } )^{\frac{1}{2}}.\) Given a symmetric positive definite \(\mathbf{B}\in \mathbb{R}^{n\times n}\), the scaled 2-norm of an \(n\)-vector \(\mathbf{p}\) is defined as \(\|\mathbf{p}\|_\mathbf{B}=(\mathbf{p}^T\mathbf{B}\mathbf{p})^{\frac{1}{2}}\). When \(\mathbf{B}=\mathbf{I}\), the norm reduces to the standard 2-norm. Denote by \(\mathbf{A}^T\), \(\mathbf{A}^{-1}\) and \(\mathbf{A}^{\dagger}\), the transpose, inverse, and pseudo-inverse of \(\mathbf{A}\), respectively, by \(\mathcal{R}(\mathbf{A})\), the column space of \(\mathbf{A}\), and by \(\mathcal{R}(A)^{\bot}\), its orthogonal complement subspace. Let \(\mathbf{b}_{\mathcal{R}(\mathbf{A})}\) and \(\mathbf{b}_{\mathcal{R}(\mathbf{A})^{\bot}}\) be the orthogonal projection of \(b\) onto \(\mathcal{R}(\mathbf{A})\) and \({\mathcal{R}(\mathbf{A})^{\bot}}\), respectively. For the matrix \(\mathbf{A}^T\mathbf{A}\), we denote its nonzero minimal eigenvalue by \(\lambda_{min}(\mathbf{A}^T\mathbf{A})\), and the maximal eigenvalue by \(\lambda_{max}(\mathbf{A}^T\mathbf{A})\). Denote by \(\mathbb{E} _{k}\), the expected value conditional of the first \(k\) iterations, i.e., \(\mathbb{E} _{k} [\cdot ]=\mathbb{E}[\cdot \mid \tau_{0}, \tau_{1},\dots,\tau_{k-1} ],\) where \(\tau_{j} (j=0,1,\dots,k-1 )\) is the \(j\)-th index set chosen at the \(j\)-th iterate, and from the law of iterated expectations, we obtain \(\mathbb{E}[\mathbb{E} _{k} [\cdot ]]=\mathbb{E}[\cdot]\).

1.2 Sketch-and-Project methods↩︎

Over the past decade, the sketch-and-project (SP) methods [22]–[24] for solving the large-scale consistent linear systems, have witnessed significant attention due to their low memory and complexity. It compresses the original system into a simplified subproblem via sketching matrices per iteration, and then projects the current iterate onto the solution space of the sketched subproblem. Let \(\mathbf{S}_k\in\mathbb{R}^{m\times r}\) be an arbitrary sketching matrix with \(r \ll m\), subject to the specified dimensions and rank. Following the Johnson-Lindenstrauss lemma [25], a solution of the reduced linear system \[\label{1-11} \mathbf{S}_k^T\mathbf{A}\mathbf{x}=\mathbf{S}_k^T\mathbf{b}\tag{2}\] is highly approximate \(\mathbf{x}\) in 1 . Viewed from the angle of projection, the next iterate \(\mathbf{x}_{k}\) of Eq.@eq:1-11 yields [26] \[\label{1-12} \mathbf{x}_{k} = \arg\min_\mathbf{x}\frac{1}{2} \|\mathbf{x} - \mathbf{x}_{k-1}\|_\mathbf{B}, \quad \text{s.t.} \quad \mathbf{S}_k^T\mathbf{A}\mathbf{x}=\mathbf{S}_k^T\mathbf{b}.\tag{3}\] As was shown in [24], [27], if \(\mathbf{S}_k\) is in the range of \(\mathbf{A}\), then \[\label{1-13} \mathbf{x}_{k}=\mathbf{x}_{k-1}+\mathbf{p}_{k},\tag{4}\] where the update \[\label{a1-5} \mathbf{p}_k = \mathbf{W} \mathbf{A}^T \mathbf{S}_k (\mathbf{S}_k^T\mathbf{A} \mathbf{W} \mathbf{A}^T\mathbf{S}_k)^{-1} \mathbf{S}_k^T \mathbf{r}_{k-1}\tag{5}\] with \(\mathbf{W}=(\mathbf{B}^T\mathbf{B})^{-1}\); otherwise, the pseudo-inverse replaces the corresponding inverse. Note that there exist infinitely many iterative schemes like 4 since the matrices \(\mathbf{S}_k\) and \(\mathbf{B}\) in the update are arbitrary. In particular, when \(\mathbf{B}=\mathbf{I}\) and \(\mathbf{S}_k=\mathbf{e}_{i_k}\) (one random column of the identity), the iterative scheme 4 becomes \[\label{1-14} \mathbf{x}_{k} =\mathbf{x}_{k-1} + \frac{\mathbf{b}_{i_k}-\mathbf{A}^{(i_{k})}\mathbf{x}_{k-1}}{\left | \left | \mathbf{A}^{(i_{k})} \right | \right |_2^{2} }(\mathbf{A}^{(i_{k})})^{T}.\tag{6}\] Obviously, they are the popular randomized Kaczmarz methods [11], [28], [29], the row-action solvers. When \(\mathbf{B}=\mathbf{A}\) and \(\mathbf{S}_k=\mathbf{A}\mathbf{e}_{j_{k}}\), i.e., \(\mathbf{A}_{(j_{k})}\) (one random column of \(\mathbf{A}\)), then 4 is of the form \[\label{1-15} \mathbf{x}_{k} =\mathbf{x}_{k-1} + \frac{\mathbf{A}_{(j_{k})}^{T}(\mathbf{b}-\mathbf{A}\mathbf{x}_{k-1})}{\left | \left | \mathbf{A}_{(j_{k})} \right | \right |_2^{2}}\mathbf{e}_{j_{k}},\tag{7}\] which are the well-known randomized coordinate descent or Gauss-Seidel methods [16], [30], [31], the column-action solvers. Note that both methods employ simple sketching matrices to minimize per-iteration computational overhead, but their convergence depends heavily on the spectral properties of the matrix \(\mathbf{A}\). When \(\lambda_{min}(\mathbf{A}^T\mathbf{A})\) is sufficiently small, the convergence factor asymptotically nears one, resulting in prohibitively large iteration complexity and severely degraded convergence performance.

To address the stagnant convergence inherent in traditional randomized methods, Brust and Saunders [27] proposed an elegant framework, called the projected linear systems solver (PLSS), for solving the consistent linear system 1 . Unlike the standard SP methods that independently regenerate the sketching matrix per iteration, the PLSS introduces an interesting process, that is, progressive sketching, which incrementally appends one column to the existing sketching matrix \(\mathbf{S}_k = [\mathbf{s}_1, \mathbf{s}_2, \dots, \mathbf{s}_k]\). More precisely, by choosing the collection of all previous residuals as the sketching matrix, the PLSS decouples the cumbersome generalized inverse updates into a simple short recurrence formula, which only needs to store and update a few vectors, effectively controlling the memory overhead. Theoretically, the method converges to an exact solution within at most \(\text{rank}(\mathbf{A})\) iterations in exact arithmetic.

It is worth pointing out that the core recurrence formulas and subspace orthogonality proofs of PLSS explicitly rely on the consistency assumption that \(\mathbf{b}\in \mathcal{R}(\mathbf{A})\). When dealing with inconsistent systems arising from observational noise or modeling deficiencies, the sketching vectors, generated from residuals, often fail to accurately capture the true error projection of the solution space. This leads to numerical ill-conditioning or complete stagnation of the denominators in the short recurrence formula, making it unavailable. Moreover, the PLSS, maintaining the completeness of the residual sketch, strictly necessitates full access to the whole matrix \(\mathbf{A}\) per iteration. This prevents the algorithm from leveraging local randomized sampling that can save operations. When encountering missing data, the PLSS does not work due to the impossibility of computing exact global residuals or executing full-matrix projections.

In this paper, we propose an efficient randomized framework, called Randomized Projected Linear Systems Solver (RPLSS), for solving Eq.@eq:1-1 whether it is under- or overdetermined, and whether it is consistent or not. The new framework lies in seamlessly integrating tailored randomized row or column selection strategies into the progressive sketching, which requires only partial matrix information per iteration. By utilizing localized sketching vectors that naturally accommodate missing or streaming data, the RPLSS maintains robust convergence and achieves substantial computational acceleration over state-of-the-art randomized solvers. Moreover, the proposed solver harmonizes the local subspace finite termination property of the PLSS with the expected linear convergence of standard randomized methods.

1.3 Our contributions↩︎

By incorporating the randomized Kaczmarz or coordinate descent methods into the RPLSS framework, we propose a family of randomized solvers for solving consistent or inconsistent linear systems. The contributions of this paper are as follows:

\(\bullet\) An RPLSS with the Gaussian Kaczmarz (RPLSS-GK) method, by taking a random linear combination consisting of residuals and standard basis vectors, is developed for solving 1 when it is consistent. The RPLSS-GK preserves the low cost while inheriting the finite-termination property of PLSS, demonstrating its robustness in row-missing scenarios.

\(\bullet\) Drawing from the double sequence idea of the randomized extended Kaczmarz method, an extended version of RPLSS-GK is derived to solve 1 when it is inconsistent. The proposed method achieves expected linear convergence without requiring any row or column paving strategies.

\(\bullet\) By randomly selecting column vectors as newly appended columns of the sketch, proportional to the square of the column norms, an RPLSS with the coordinate descent (RPLSS-CD) method is devised for solving 1 when it is inconsistent. Utilizing the orthogonal projection of PLSS to accumulate historical updates, the RPLSS-CD method effectively overcomes the slow convergence caused by single-component updates in traditional coordinate descent methods.

\(\bullet\) Combining the current residual with randomly selected standard basis vectors flexibly, an RPLSS with the Gaussian coordinate descent (RPLSS-GCD) method is proposed. This method strikes an elegant balance between local coordinate updates and global residual feedback, which substantially accelerates its convergence on ill-conditioned and inconsistent problems.

The organization of this paper is as follows. In Section 2, we first devise the RPLSS framework for solving Eq.@eq:1-1 , and then propose two fast randomized row-action solvers, including RPLSS-GK and RPLSS-GEK, for the consistent and inconsistent cases, respectively. In Section 3, we propose two randomized column-action solvers. The convergence of all methods is also established. In Section 4, we provide abundant numerical results on large-scale sparse matrices, ill-conditioned systems, and applications with severe data missing, comparing the proposed algorithms against state-of-the-art methods to validate their superiority. Finally, we conclude this paper by giving some remarks.

2 RPLSS with Gaussian Kaczmarz methods↩︎

In this section, we seamlessly integrate randomized row- or column-selection strategies into the progressive sketching process, and then propose a randomized PLSS framework for the linear system 1 . Unlike the standard PLSS framework, the proposed framework only requires partial information of the coefficient matrix per iteration, which is particularly favorable for large-scale, memory-constrained, or streaming data problems. Moreover, by incrementally storing historical information in the sketching matrix, the RPLSS also preserves the finite-termination property of PLSS. Indeed, by randomly constructing different sketch vectors \(\mathbf{s}_k\), the various effective randomized methods can be developed for solving the consistent or inconsistent Eq.@eq:1-1 .

As stated in [27], the updates \(\mathbf{p}_k=\mathbf{P}_{k-1}\mathbf{g}_{k-1}+\gamma_{k-1}\mathbf{y}_k\) are orthogonal, with \(\mathbf{P}_k=[\mathbf{p}_1,\mathbf{p}_2,\dots,\mathbf{p}_k]\) and \(\mathbf{y}_k=\mathbf{A}^T\mathbf{s}_k\). Following this orthogonality, it is trivial to see that \[\mathbf{g}_{k-1}=-\gamma_{k-1}\mathbf{\Theta}_{k-1} ^{-1}\mathbf{P}_{k-1}^T\mathbf{y}_k,\] where \(\mathbf{\Theta}_{k-1}=\mathbf{P}_{k-1}^T\mathbf{P}_{k-1}\). Moreover, since \(\mathbf{s}_k^T\mathbf{A}\mathbf{p}_k=\mathbf{s}_k^T\mathbf{r}_{k-1}\), it follows that \[\gamma_{k-1}=\frac{\mathbf{s}_k^T(\mathbf{b}-\mathbf{A}\mathbf{x}_{k-1})}{\|\mathbf{A}^T\mathbf{s}_k\|_2^2-\|\mathbf{d}_{k-1}\|_2^2},\] where \(\mathbf{d}_{k-1}=\mathbf{\Theta}_{k-1} ^{-1/2}\mathbf{P}_{k-1}^T\mathbf{A}^T\mathbf{s}_k\). Here, the \(\mathbf{\Theta}_{k-1}^{-1/2}\)denotes the inverse of the matrix square root of the diagonal matrix \(\mathbf{\Theta}_{k-1}\). Therefore, the updates \(\mathbf{p}_k\) are of the form

\[\mathbf{p}_k=\frac{\mathbf{s}_k^T(\mathbf{b}-\mathbf{A}\mathbf{x}_k)}{\|\mathbf{y}_k\|_2^2-\|\mathbf{d}_{k-1}\|_2^2}\mathbf{A}^T\mathbf{s}_k-\frac{\mathbf{s}_k^T(\mathbf{b}-\mathbf{A}\mathbf{x}_k)\mathbf{P}_{k-1}\mathbf{\Theta}_{k-1} ^{-1}\mathbf{P}_{k-1}^T}{\|\mathbf{y}_k\|_2^2-\|\mathbf{d}_{k-1}\|_2^2}\mathbf{A}^T\mathbf{s}_k.\] If \(\mathbf{W}\neq \mathbf{I}\), then we let \(\mathbf{W}=\mathbf{R^T_{\mathbf{W}}}\mathbf{R_{\mathbf{W}}}\), that is, Cholesky decomposition. Define \[\mathbf{x}_{k-1}=\mathbf{x}_{k-1}^\mathbf{W}\mathbf{R}_\mathbf{W}^{T},\quad\mathbf{p}_k=\mathbf{p}_k^\mathbf{W}\mathbf{R}_\mathbf{W}^{T},\quad\mathbf{A}_\mathbf{W} = \mathbf{A}\mathbf{R}_\mathbf{W}^T.\] Then it follows that \[\mathbf{p}_k^\mathbf{W} = \gamma_{k-1}^\mathbf{W} \Bigl( \mathbf{y}_k^\mathbf{W} - \mathbf{P}_{k-1}^\mathbf{W} \bigl((\mathbf{P}_{k-1}^{\mathbf{W}})^T\mathbf{P}_{k-1}^\mathbf{W}\bigr)^{-1} (\mathbf{P}_{k-1}^{\mathbf{W}})^T\mathbf{y}_k^\mathbf{W} \Bigr),\] where \[\mathbf{y}_k^\mathbf{W} = \mathbf{A}_\mathbf{W}^T\mathbf{s}_{k},\, \gamma_{k-1}^\mathbf{W} = \frac{\mathbf{s}_k^T(\mathbf{b}-\mathbf{A}\mathbf{x}_{k-1})}{\|\mathbf{y}_k^\mathbf{W}\|^2 - \|\mathbf{d}_{k-1}^\mathbf{W}\|^2},\, \mathbf{d}_{k-1}^\mathbf{W} = \bigl((\mathbf{P}_{k-1}^{\mathbf{W}})^T\mathbf{P}_{k-1}^\mathbf{W}\bigr)^{-1/2}(\mathbf{P}_{k-1}^{\mathbf{W}})^T\mathbf{y}_k^\mathbf{W}.\] From the above results, the update \(\mathbf{p}_k\) can be rewritten as \[\mathbf{p}_k=\gamma_{k-1}^{\mathbf{W}}(\mathbf{W}\mathbf{A}^T\mathbf{s}_k-\mathbf{P}_{k-1}(\mathbf{\Theta}_{k-1}^{\mathbf{W}}) ^{-1}\mathbf{P}_{k-1}^T\mathbf{A}^T\mathbf{s}_k),\] where \(\mathbf{d}_{k-1}^{\mathbf{W}}=\mathbf{\Theta}_{k-1} ^{-1/2}\mathbf{P}_{k-1}^T\mathbf{A}^T\mathbf{s}_k\) and \(\mathbf{\Theta}_{k-1}^{\mathbf{W}}=\mathbf{P}_{k-1}^T\mathbf{W}^{-1}\mathbf{P}_{k-1}\). Remarkably, the PLSS with residual sketches requires the whole matrix-vector product \(\mathbf{A}^T \mathbf{r}_{k}\) per iteration. When dealing with large-scale problems, this full-dimensional operation incurs high computational costs and slow convergence.

By replacing the deterministic operation with a randomized subsampling of the rows or columns of \(\mathbf{A}\), such as, the sketch vector \(\mathbf{s}_k\), which consists of a random linear combination of the previous residual and the standard basis vectors, we introduce a randomized PLSS (RPLSS) framework for Eq.@eq:1-1 . While sacrificing a marginal degree of precision, this framework significantly reduces the per-iteration cost, rendering it highly effective for large-scale applications or environments where the complete data is inaccessible. Now, we present the following RPLSS framework.

As seen, the RPLSS framework gives rise to a family of different algorithms depending on the choice of the sketch vector \(\mathbf{s}_k\). Here, we mainly consider two different the sketch vectors \(\mathbf{s}_k=\sum\limits_{i\in \tau_k}\mathbf{r}^{(i)}_{k-1}\mathbf{e}_i\) or \(\mathbf{s}_k=\sum\limits_{i\in\tau_k}\mathbf{A}\mathbf{e}_i\). The former is a random linear combination of the previous residual and the standard basis vectors, while the latter is a random submatrix that consists of some columns of \(\mathbf{A}\).

From Theorem 3.1, Corollaries 3.2 and 3.3 presented in [27], we make the following statement.

2.1 RPLSS-GK↩︎

The randomized Kaczmarz-type methods, a leading class of randomized solvers, have received significant attention due to their simplicity and efficiency. For instance, Needell and Tropp [14] proposed a randomized block Kaczmarz method that greatly accelerates the convergence of the standard RK, but suffers from the high computational expense of evaluating pseudo-inverses. Then, Necoara [32] pioneered the randomized average block Kaczmarz (RaBK) method without any pseudo-inverses. Moreover, Bai and Wu [28] proposed a fast greedy randomized Kaczmarz (GRK) method, which adaptively selects rows based on residual magnitudes.

Overall, these randomized Kaczmarz solvers demonstrate that utilizing only a subset of row information can yield efficient convergence when solving large-scale and consistent linear systems. Motivated by this, we devise a similar strategy that randomly samples a portion of the residual vector to form the sketch vector \(\mathbf{s}_k\), thereby circumventing the full-dimensional residual operations mandated by PLSS per iteration. Here, we set \(\mathbf{s}_k=\sum\limits_{i\in \tau_k}\mathbf{r}^{(i)}_{k-1}\mathbf{e}_i\), a linear combination of previous residual and standard basis vectors, where \(\tau_k\) is a randomly selected set of row indices. By integrating this strategy into the RPLSS framework, we propose an efficient RPLSS variant based on the Gaussian Kaczmarz method for solving the consistent linear system 1 .

We now state some necessary lemmas before proceeding to the convergence analysis.

From the above lemmas, we will establish the following results regarding the convergence rate of Algorithm 2. To simplify the analysis, we let \(\mathbf{W}=\mathbf{I}\).

2.2 RPLSS-GEK↩︎

When the linear system 1 is inconsistent, Needell [12] showed that RK fails to converge to the least-norm least-squares solution. To solve this problem, Zouzias and Freris [18] first proposed a randomized extended Kaczmarz (REK) method. More concisely, the REK solver iteratively updates two vectors \(\mathbf{z}_k\) and \(\mathbf{x}_k\) (\(\mathbf{A}^*\mathbf{z}=0\) for \(\mathbf{z}_k\) and \(\mathbf{A}\mathbf{x}=\mathbf{b}-\mathbf{z}_k\) for \(\mathbf{x}_k\)) by \[\mathbf{z}_{k+1} =(I_m-\frac{\mathbf{A}_{(j_k)}\mathbf{A}_{(j_k)}^T}{\|\mathbf{A}_{(j_k)}\|_2^2})\mathbf{z}_{k},\] \[\mathbf{x}_{k+1} =\mathbf{x}_k+\frac{{\mathbf{b}}^{(i_{k})}-\mathbf{z}_k^{(i_k)}-\mathbf{A}^{(i_{k})}\mathbf{x}_{k}}{\| \mathbf{A}^{(i_{k})}\|_{2} ^{2} }(\mathbf{A}^{(i_{k})})^{T}.\] To further improve its convergence, a block version of REK has been proposed in [19]. Moreover, Du et al. devised a novel randomized extended average block Kaczmarz method [15] for finding the least-squares solution of Eq.@eq:1-1 without involving the pseudo-inverse. Ke et al. [20] derived a greedy block extended Kaczmarz (GBEK) method with the average block projection technique. Besides, Bai and Wu [17] further extended the greedy randomized Kaczmarz framework to inconsistent systems.

From these randomized extended Kaczmarz solvers, we see that double‑sequence iteration is capable of circumventing the convergence failure of single‑sequence iteration on large‑scale inconsistent linear systems. Moreover, when the linear system 1 is inconsistent, PLSS also fails to converge to the least-squares solution. To address this drawback, we develop an extended version of the RPLSS Gaussian Kaczmarz method, which addresses the slow convergence issue of conventional extended Kaczmarz methods. Specifically, referencing the double sequence iteration of REK, we maintain two sequences: one for \(\mathbf{z}_k\) that drives \(\mathbf{A}^T\mathbf{z}\) toward zero, and one for \(\mathbf{x}_k\) that solves \(\mathbf{A}\mathbf{x} = \mathbf{b} - \mathbf{z}_k\).

In the RPLSS framework, the update for \(\mathbf{z}_k\) is obtained by solving the projected systems \[\label{2-3} \bar{\mathbf{S}}_k^T\mathbf{A}^T\mathbf{z}_k=0 ,\tag{10}\] where \(\bar{\mathbf{S}}_k\) is a sketch matrix. This yields the update \(\mathbf{z}_k=\mathbf{z}_{k-1}-\bar{\mathbf{p}}_k\) with \[\bar{\mathbf{p}}_k=\bar{\gamma}_{k-1}(\mathbf{A}\bar{\mathbf{s}}_k-\bar{\mathbf{P}}_{k-1}\bar{\mathbf{\Theta}}_{k-1} ^{-1}\bar{\mathbf{P}}_{k-1}^T\mathbf{A}\bar{\mathbf{s}}_k)\] where \(\bar{\gamma}_{k-1}=\frac{\bar{\mathbf{s}}_k^T\mathbf{A}^T\mathbf{z}_{k-1}}{\|\mathbf{A}\bar{\mathbf{s}}_k\|_2^2-\|\bar{\mathbf{d}}_{k-1}\|_2^2}\) with \(\bar{\mathbf{d}}_{k-1}=\bar{\mathbf{\Theta}}_{k-1} ^{-1/2}\bar{\mathbf{P}}_{k-1}^T\mathbf{A}\bar{\mathbf{s}}_k\). The detail of \(\bar{\mathbf{p}}_k\) is provided in Appendix 6. Similarly, the update for \(\mathbf{x}_k\) follows the RPLSS‑GK solver, except that the right side is replaced with \(\mathbf{b} - \mathbf{z}_k\) and the history of previous updates is no longer used, i.e., \(\mathbf{P}_{k-1}=\mathbf{0}\).

From what was said above, we now present the RPLSS Gaussian extended Kaczmarz (RPLSS-GEK) method as follows.

An important result for guaranteeing the convergence of Algorithm 3 is stated in the following theorem.

\[\label{2-4} \mathbb{E}\|\mathbf{z}_k-\mathbf{b}_{\mathcal{R}(\mathbf{A})^{\bot}}\|_2^2\leq\hat{\beta} ^k\|\mathbf{b}_{\mathcal{R}(\mathbf{A})}\|_2^2,\tag{11}\] where \(\hat{\beta}=\sup\limits_k\left(1-\frac{l\lambda_{min}^2(\mathbf{A}\mathbf{A}^T)}{\sin^2(\tilde{\theta}_{k-1})\lambda_{max}^2(\mathbf{A}\mathbf{A}^T)}\right)\).

From Theorem 2, we give the following theorem for guaranteeing the convergence of Algorithm 3.

3 RPLSS with coordinate descent method↩︎

Note that the column-action solver is an attractive alternative for linear systems 1 , among which the randomized coordinate descent (RCD) method, also called randomized Gauss-Seidel, is well-established. Unlike the RK, which sequentially projects the iterate onto the hyperplanes associated with individual rows, the RCD operates on the column space by successively minimizing the objective residual norm. This structural difference enables RCD to converge unconditionally to a least-squares solution for inconsistent systems, where the RK method otherwise suffers from non-convergent oscillations. However, the RCD updates only one component per iteration and discards historical information, which may hinder its convergence on ill-conditioned and large-scale problems. To develop a column-action algorithm within the randomized PLSS framework, we redefine the selection of the sketch column \(\mathbf{s}_k\) and the weight matrix \(\mathbf{W}\). This gives a class of alternative randomized PLSS algorithms based on the column actions.

3.1 RPLSS-CD↩︎

In the RPLSS framework, let \(\mathbf{s}_k=\mathbf{A}\mathbf{e}_i\) with \(\mathbf{e}_i\) being sampled with probability \(\mathbb{P}\)(column=\(i)=\) \(\frac{\|\mathbf{A}_{(i)}\|_2^2}{\|\mathbf{A}\|_F^2}\) and \(\mathbf{W}=(\mathbf{A}^T\mathbf{A})^{-1}\). Then we derive the following RPLSS coordinate descent (RPLSS-CD) method for solving Eq.@eq:1-1 .

For the convergence property of the RPLSS-CD method, we can establish the following theorem.

Next, we give a more generalized RPLSS-CD method, which is more stable than RPLSS-CD.

Regarding the convergence of the Algorithm [algo-3-2], we establish the following theorem.

3.2 RPLSS-GCD↩︎

As presented in the previous section, the RPLSS‑CD randomly selects a single column per iteration to construct the sketch. While this strategy is memory‑efficient, it updates only one coordinate at a time and ignores the global error information contained in the residual vector. To making fully use of the residuals, we construct a sketch column \(\mathbf{s}_k\) by randomly linearly combining the residuals with the standard basis vectors, which can vastly accelerate the convergence of ill-conditioned problems. Specifically, in the RPLSS framework, we set \(\mathbf{s}_k=\sum_{i\in \tau_k}\mathbf{A}\mathbf{r}^{(i)}_{k-1}\mathbf{e}_i\) and \(\mathbf{W}=(\mathbf{A}^T\mathbf{A})^{-1}\), then develop the following RPLSS Gaussian coordinate descent (RPLSS-GCD) method for solving Eq.@eq:1-1 .

Next, we show the exponential convergence of the Algorithm 4.

4 Numerical Experiments↩︎

In this section, we record some numerical results to illustrate the feasibility and validity of the proposed methods compared with some state-of-the-art randomized solvers for consistent and inconsistent linear systems 1 . We first test the methods with randomly generated coefficient matrices or those coming from the University of Florida Sparse Matrix Collection. Then we present additional numerical results on the X-ray computed tomography problem to further illustrate the effectiveness of the proposed method in comparison with the solvers mentioned above. Denote by “IT", the number of iterations, and by”CPU", the elapsed computing time in seconds. For each method, we report the mean computing time in seconds and the mean number of iterations based on their average values of \(5\) repeated tests. Additionally, all experiments have been carried out in MATLAB 2021b on a personal computer with Inter(R) Core(TM) i5-13500HX @2.5GHz and 32.00 GB memory. In Algorithms 2, 3 and 4, the parameter \(l\) is set to \(0.5\). The symbol \(\dagger\) indicates that the solver fails to converge after reaching 1000 seconds. All test methods are as follows:
\(\star\) GRK: Greedy randomized Kaczmarz method [28].
\(\star\) RaBK: Randomized average block Kaczmarz method [32].
\(\star\) FGBK: Fast greedy block Kaczmarz method [34].
\(\star\) RBK(50): Randomized block Kaczmarz method [14] with the block size set to be 50.
\(\star\) REABK: Randomized extended average block Kaczmarz method[15].
\(\star\) GBEK: Greedy block extended Kaczmarz method[20].
\(\star\) RCD: Randomized coordinate descent method [30].
\(\star\) GRCD: Greedy randomized coordinate descent method[31].

Example 1. Consider solving the Eq.@eq:1-1 with its coefficient matrix given by the MATLAB function \(sprandn\). Moreover, we randomly generate a solution \(\mathbf{x}_\star\) by the MATLAB function \(randn\), with \(\mathbf{x}_\star=randn(n,1)\). For consistent linear systems, the right-hand side vector \(\mathbf{b}\) is constructed as \(\mathbf{b}=\mathbf{A}\mathbf{x}_\star\). For inconsistent linear systems, we set \(\mathbf{b}=\mathbf{A}\mathbf{x}_\star+r\), where \(r\) is a nonzero vector belonging to the null space of \(\mathbf{A}^T\). The stopping criteria of consistent linear systems is the error norm at the current iterate \(\mathbf{x}_{k}\) satisfying \[\|\mathbf{x}_{k}-\mathbf{x}_{\star} \| _{2}<10^{-6},\] and the stopping criteria of inconsistent linear systems is the relative error norm at the current iterate \(\mathbf{x}_{k}\) satisfying \[\mathrm{RSE}:= \frac{\|\mathbf{x}_k-\mathbf{x}_{\star}\|_2^2}{\|\mathbf{x}_{\star}\|_2^2}< 10^{-6}.\] For the RaBK solver, we set \(\tau=m/\|\mathbf{A}\|_2^2\) and \(\alpha=1.95L_k\) with \[L_k = \frac{\sum_{i\in \tau_k} \bar{\omega}_i^k ( (\mathbf{A}^{(i)})^T \mathbf{x}^k - \mathbf{b}^{(i)})^2}{\left\| \sum_{i\in \tau_k} \bar{\omega}_i^k( (\mathbf{A}^{(i)})^T \mathbf{x}^k - \mathbf{b}^{(i)}) \mathbf{A}^{(i)} \right\|^2}.\]

We implement all tested methods to solve consistent linear systems 1 , and then record the obtained numerical results in Tables 1 and 2. A rough conclusion drawn from these tables is that, RPLSS-GK significantly outperforms GRK, RaBK, FGBK, and RBK(50) in all tests of consistent linear systems (overdetermined, underdetermined, and of different scales), making it the most competitive algorithm in the current table. The computation time of RPLSS-GK is about \(2-3\) times faster than FGBK, and \(1-2\) orders of magnitude faster than RaBK and RBK(50). Especially in ultra-large-scale problems with millions of rows, RPLSS-GK only takes \(3-20\) seconds, while other methods take tens to hundreds of seconds. This demonstrates RPLSS-GK is more favorable for big data problems. The convergence curves of all methods with \(300000\times 40000\) and \(50000\times 100000\) were displayed in Fig.6. It can be clearly seen from the figure that RPLSS-GK converges faster than other existing methods.

Table 1: Numerical results of Example [Example-1] when Eq.[1-1] is consistent

\(m \times n\)		\(20000\times5000\)	\(300000\times40000\)	\(600000\times50000\)	\(1000000\times100000\)
	IT	7896	\(\dagger\)	\(\dagger\)	\(\dagger\)
GRK	CPU	16.0886	\(\dagger\)	\(\dagger\)	\(\dagger\)
	IT	1022	1082	846	958
RaBK	CPU	0.3399	11.0386	17.0775	20.9910
	IT	31	84	60	53
FGBK	CPU	0.0929	2.7069	6.0688	9.6621
	IT	73	286	155	38
RBK(50)	CPU	0.1363	7.4937	5.7692	30.8101
	IT	17	29	22	18
RPLSS-GK	CPU	0.0322	0.8366	2.3663	3.4091

0.7mm

Table 2: Numerical results of Example [Example-1] when Eq.[1-1] is consistent

\(m \times n\)		\(5000\times10000\)	\(50000\times100000\)	\(100000\times250000\)	\(200000\times500000\)
	IT	269100	\(\dagger\)	\(\dagger\)	\(\dagger\)
GRK	CPU	813.1132	\(\dagger\)	\(\dagger\)	\(\dagger\)
	IT	1119	1178	1211	1212
RaBK	CPU	4.2261	22.0106	98.5854	375.4899
	IT	435	641	504	254
FGBK	CPU	0.7145	5.2255	19.7438	43.5896
	IT	8540	9081	5965	5973
RBK(50)	CPU	30.7323	160.3646	250.6066	446.2368
	IT	127	143	112	112
RPLSS-GK	CPU	0.2353	1.7733	5.0777	20.3680

0.7mm

For the inconsistent case, we report the numerical results regarding the iterations and computing time in Table 3. As shown in the table, RPLSS-GCD and RPLSS-GEK converge stably in all test scales, and the number of iterations and computing time are much better than other methods, especially showing great application potential in very large-scale problems. Among traditional methods, only the GBEK method converges on all scale matrices, the REABK algorithm fails on the largest scale case, and RCD and GRCD also fail at three and two cases, respectively. Moreover, we displayed the convergence curves of all methods with the case \(m=10000\) and \(n=2000\) in Fig.[Fig:2], which coincide with the above observations.

Table 3: Numerical results of Example [Example-1] when Eq.[1-1] is inconsistent

\(m \times n\)		\(10000\times2000\)	\(50000\times20000\)	\(200000\times80000\)	\(1000000\times100000\)
	IT	5855	29953	78568	\(\dagger\)
REABK	CPU	2.1651	34.0505	232.9871	\(\dagger\)
	IT	274	5920	818	175
GBEK	CPU	0.1826	20.8799	34.1692	51.2237
	IT	54282	\(\dagger\)	\(\dagger\)	\(\dagger\)
RCD	CPU	14.4246	\(\dagger\)	\(\dagger\)	\(\dagger\)
	IT	12708	280691	\(\dagger\)	\(\dagger\)
GRCD	CPU	3.8332	790.5998	\(\dagger\)	\(\dagger\)
	IT	2000	\(\dagger\)	\(\dagger\)	\(\dagger\)
RPLSS-CD(10)	CPU	3.7865	\(\dagger\)	\(\dagger\)	\(\dagger\)
	IT	30	54	41	23
RPLSS-GCD	CPU	0.0193	0.1991	1.4567	6.6872
	IT	40	113	65	25
RPLSS-GEK	CPU	0.0457	0.6718	2.6305	9.4435

0.4mm

Example 2. Consider solving the consistent linear system 1 with its coefficient matrix coming from the University of Florida Sparse Matrix Collection [35]. The chosen matrices are mainly used for statistical, mathematical, combinatorial, and linear programming problems. Moreover, we randomly generate a solution \(\mathbf{x}_\star\) by the MATLAB function \(randn\), with \(\mathbf{x}_\star=randn(n,1)\), such that the right-hand side vector \(\mathbf{b}\) is set to be \(\mathbf{A}\mathbf{x}_\star\). The stopping criteria of all tested methods are the error norm at the current iterate \(\mathbf{x}_{k}\) satisfying \[\|\mathbf{x}_{k}-\mathbf{x}_{\star} \| _{2}<10^{-6}.\]

Table 4: Numerical results of Example [Example-2]

Matrix		Franz11	GL7d14	relat8	testbig
\(m \times n\)		\(47104\times30144\)	\(171375\times47271\)	\(345688\times12347\)	\(17613\times31223\)
Density		0.02%	0.02%	0.03%	0.01%
GRK	IT	\(\dagger\)	\(\dagger\)	\(\dagger\)	1113996
	CPU	\(\dagger\)	\(\dagger\)	\(\dagger\)	971.3152
RaBK	IT	20122	98756	87272	\(\dagger\)
	CPU	11.4521	87.3713	38.6179	\(\dagger\)
FGBK	IT	484	1541	1766	166061
	CPU	0.6100	9.2269	12.3514	81.9880
RBK(50)	IT	1912	4083	1453	228235
	CPU	11.8855	70.2431	16.8130	834.0728
RPLSS-GK	IT	81	98	107	50
	CPU	0.2017	0.9340	1.3551	0.7952

0.3mm

Table 5: Numerical results of Example [Example-2]

Matrix		sls	kneser_10_4_1	degme	ins2
\(m \times n\)		\(1748122\times62729\)	\(349651\times330751\)	\(185501\times659415\)	\(309412\times309412\)
Density		0.006%	0.0009%	0.007%	0.003%
GRK	IT	\(\dagger\)	\(\dagger\)	\(\dagger\)	\(\dagger\)
	CPU	\(\dagger\)	\(\dagger\)	\(\dagger\)	\(\dagger\)
RaBK	IT	\(\dagger\)	\(\dagger\)	800	\(\dagger\)
	CPU	\(\dagger\)	\(\dagger\)	45.0287	\(\dagger\)
FGBK	IT	\(\dagger\)	14914	2622	\(\dagger\)
	CPU	\(\dagger\)	200.3864	43.8835	\(\dagger\)
RBK(50)	IT	429	\(\dagger\)	2554	\(\dagger\)
	CPU	24.2707	\(\dagger\)	160.1261	\(\dagger\)
RPLSS-GK	IT	290	804	103	381
	CPU	14.2922	39.6706	3.2724	11.6757

0.1mm

We run the test methods on some practical problems, and then list the numerical results in Tables 4 and 5. As shown in these tables, GRK fails to converge in most cases. Although RaBK, FGBK, and RBK converge on some problems, the number of iterations and computing time required are much higher than those of RPLSS-GK. Traditional Kaczmarz-type methods generally fail on very large-scale matrices (sls, ins2) and some singular matrices, while RPLSS-GK converges stably on all test matrices, showing stronger robustness. A rough conclusion is that the RPLSS-GK method significantly outperforms the other solvers in terms of the number of iterations and CPU time, especially for ill-conditioned and large-scale sparse matrices, where its advantages are even more pronounced. Particularly, the speed-up of RPLSS-GK compared with RaBK is at least \(13.71\) in the case of degme, and at most \(93.54\) in the case of GL7d14, compared with FGBK is at least \(3.02\) in the case of Franz11, and at most \(103.10\) in the case of testbig, and compared with RBK is at least \(1.69\) in the case of sls, and at most \(1048.88\) in the case of testbig. Fig.8 plots the convergence curves of all tested methods under the Franz11 and degme test cases. It is evident that the proposed RPLSS-GK solver achieves superior convergence performance compared with the other randomized solvers, consistent with the aforementioned observations.

Example 3. Consider solving the inconsistent linear system 1 with its coefficient matrix coming from the University of Florida Sparse Matrix Collection [35]. The chosen matrices are mainly used for combinatorial and linear programming problems. Moreover, we randomly generate a solution \(\mathbf{x}_\star\) by the MATLAB function \(randn\), with \(\mathbf{x}_\star=randn(n,1)\), such that the right-hand side vector \(\mathbf{b}\) is set to be \(\mathbf{A}\mathbf{x}_\star+r\), where \(r\in\) null\((\mathbf{A}^T)\). The stopping criteria of all tested methods are the relative error norm at the current iterate \(\mathbf{x}_{k}\) satisfying \[\mathrm{RSE}:= \frac{\|\mathbf{x}_k-\mathbf{x}_{\star}\|_2^2}{\|\mathbf{x}_{\star}\|_2^2}< 10^{-6}.\]

Table 6: Numerical results of Example [Example-3]

Matrix		ch7-9-b5	D6-6	n4c6-b7	cage12	abtaha2
\(m\)		\(423360\)	\(120576\)	\(163215\)	\(130228\)	\(37932\)
\(n\)		\(317520\)	\(23740\)	\(104115\)	\(130228\)	\(331\)
Density		0.12%	0.36%	0.01%	0.01%	1.09%
REABK	IT	39233	40961	6778	737472	27128
	CPU	353.2802	61.2662	13.6403	956.9433	12.8866
GBEK	IT	359	871	117	547	925
	CPU	6.4408	2.8610	0.9547	4.6600	1.0560
RPLSS-GEK	IT	33	113	19	160	285
	CPU	0.9433	0.6353	0.2159	1.5078	0.5514

1.6mm

In Tables 6-7, we report the numerical results regarding the iterations and computing time. As shown in these tables, one can reasonably conclude that the RPLSS-GEK solver outperforms all other tested solvers. RPLSS-GEK converges stably on all test matrices, with significantly fewer iterations and shorter computation time compared to existing methods like REABK and GBEK. It demonstrates excellent robustness and efficiency, especially for highly ill-conditioned matrices or ultra-large-scale sparse matrices (such as Delor64K). More precisely, compared with the REABK solver, RPLSS-GEK brings an amazing speed-up of \(23.37\) times to \(634.66\) times. Compared with the efficient GBEK solver, RPLSS-GEK also maintains a stable acceleration of \(1.91\) to \(8.63\) times. Additionally, we depicted the convergence curves of all tested methods with the cases abtaba2 and lp_nug20 in Fig.9. As seen, the RPLSS-GEK solver still work more effectively than the REABK and GBEK solvers, respectively, which coincides with the abovementioned conclusion.

Table 7: Numerical results of Example [Example-3]

Matrix		lp_nug20	bibd_49_3	GL7d24	rlfdual	Delor64K
\(m\)		\(15240\)	\(1176\)	\(21074\)	\(8052\)	\(64719\)
\(n\)		\(72600\)	\(18424\)	\(105054\)	\(74970\)	\(1785345\)
Density		0.03%	2.55%	0.03%	0.05%	0.0006%
REABK	IT	7632	1433	135167	67152	\(\dagger\)
	CPU	11.6339	0.6571	130.6364	130.6364	\(\dagger\)
GBEK	IT	211	52	8072	1319	\(\dagger\)
	CPU	0.4528	0.0586	32.0111	3.3230	\(\dagger\)
RPLSS-GEK	IT	32	17	1543	722	341
	CPU	0.1067	0.0197	3.7095	1.2273	9.7724

1.6mm

Example 4. Consider solving the problem of X-ray computed tomography, which is one of the classic application areas of the Kaczmarz method. The Matlab package AIR Tools II [36] is used to simulate the above problems, so as to verify the effectiveness and performance of the proposed method for different linear systems. In this section, two sets of comparative experiments are designed, corresponding to consistent and inconsistent linear systems, respectively. First, we generated a test image of size 256 by 256 based on the \(paralleltomo\) function in the AIR Tools II package and recorded the measurements of the angles 0:0.5:179. The number of sensors per projection line is 280. The right-hand side vector is \(\mathbf{b}=\mathbf{A}\mathbf{x}_{\star}\), where the unique solution of the consistent linear systems 1 is obtained by reshaping the \(256\times256\) Shepp-Logan medical phantom. We verified the effectiveness of the algorithms under three scenarios: complete projection data, \(10\%\) missing projection data, and \(30\%\) missing projection data. The stopping criteria of all tested methods are the relative error norm at the current iterate \(\mathbf{x}_{k}\) satisfying \[\mathrm{RSE}:= \frac{\|\mathbf{x}_k-\mathbf{x}_{\star}\|_2^2}{\|\mathbf{x}_{\star}\|_2^2}< 10^{-2}.\]

In Tables 8-10, we list the number of iteration steps, the computing times, and the relative solution errors for all tested methods. It appears from this table that the RPLSS-GK significantly outperforms GRK, RaBK, FGBK, and RBK(50) in terms of the number of iterations and computation time. Especially for large-scale sparse matrices, its speedup can reach tens or even hundreds of times. As matrix row information gradually becomes missing, the number of iterations and iteration time of other methods increase sharply, while the number of iterations of RPLSS-GK increases slowly and the iteration time remains basically unchanged. This indicates that RPLSS-GK has good adaptability to missing matrix rows and exhibits good robustness. This is because RPLSS-GK constructs a sketch matrix using historical residuals, which has a natural resistance to interference in row sampling, and can still maintain efficient convergence even when \(30\%\) of the matrix rows are missing. Therefore, RPLSS-GK is particularly suitable for solving large-scale linear systems with incomplete rows and uneven sampling, which are common in real-world data. The pictures of projection data are displayed in Figs.[Fig:5]-[Fig:7].

Table 8: Numerical results of Example [Example-4] for \(100520\times65536\) matrix \(A\)

Algorithm		IT	CPU	RSE
GRK		20389	\(\dagger\)	1.91e-2
RaBK		13439	23.1365	9.99e-3
FGBK		123	5.5835	9.99e-3
RBK(50)		123	7.0292	9.22e-3
RPLSS-GK		26	1.1416	9.58e-3

Table 9: Numerical results of Example [Example-4] for \(90\%\) \(100520\times65536\) matrix \(A\)

Algorithm		IT	CPU	RSE
GRK		22343	\(\dagger\)	1.98e-2
RaBK		17002	38.2633	9.99e-3
FGBK		184	6.9108	9.99e-3
RBK(50)		263	13.4510	9.99e-3
RPLSS-GK		29	1.0590	9.68e-3

Table 10: Numerical results of Example [Example-4] for \(70\%\) \(100520\times65536\) matrix \(A\)

Algorithm		IT	CPU	RSE
GRK		27802	\(\dagger\)	2.27e-2
RaBK		24596	50.5284	9.99e-3
FGBK		581	15.7611	9.99e-3
RBK(50)		681	28.7061	9.99e-3
RPLSS-GK		42	1.2160	9.94e-3

Then, we generated a test image of size 128 by 128 based on the \(paralleltomo\) function in the AIR Tools II package and recorded the measurements of the angles 0:0.5:179. The number of sensors per projection line is 140. The right-hand side vector is \(\mathbf{b}=\mathbf{A}\mathbf{x}_{\star}+r\), where the unique solution of the inconsistent linear systems 1 is obtained by reshaping the \(128\times128\) Shepp-Logan medical phantom and \(r\) is a nonzero vector in the null space generated by the MATLAB function \(null\).

In Table 11, we list the number of iteration steps, the computing times, and the relative solution errors for all tested methods. As was shown, RPLSS-GCD and RPLSS-GEK methods show better performance than REABK, GBEK, RCD, GRCD, and PLSSRCD(10) in terms of computing time. RCD, GRCD, and RPLSS-CD (10) all failed to converge or had insufficient accuracy. RPLSS-GCD achieves extreme computational efficiency and robustness in ultra-large ill-conditioned problems, making it the best method in this test.

Table 11: Numerical results of Example [Example-4] for \(50260\times16384\) matrix \(A\)

Algorithm		IT	CPU	RSE
REABK		18963	21.8153	9.99e-3
GBEK		453	8.3650	9.99e-3
RCD		82109	\(\dagger\)	5.96e+2
GRCD		77987	\(\dagger\)	1.04e-1
RPLSS-CD(10)		15932	\(\dagger\)	1.57e-1
RPLSS-GCD		27	0.5096	9.66e-3
RPLSS-GEK		240	5.0662	9.99e-3

5 Conclusions↩︎

In this paper, we propose a randomized projected linear system solver framework to solve Eq. 1 . Unlike the PLSS method that requires access to the whole coefficient matrix in each iteration, our method utilizes a random row or column selection strategy, using only partial information in each iteration. This significantly reduces memory requirements and computational cost per iteration, making the proposed method particularly suitable for large-scale, memory-constrained, or data-streaming applications. The main contributions are listed as follows.

• For consistent systems, we proposed the RPLSS with the Gaussian Kaczmarz (RPLSS‑GK) method, which constructs each sketch column as a random linear combination of residuals and standard basis vectors. This design preserves the low memory footprint of Kaczmarz‑type iterations while inheriting the finite-termination property of PLSS. Numerical experiments demonstrate that RPLSS‑GK is remarkably robust even under row‑missing scenarios.

• Inspired by the double‑sequence idea of the randomized extended Kaczmarz method, we derived an extended version of RPLSS‑GK to handle inconsistent systems. This extension achieves expected linear convergence to the least‑squares solution without requiring any row or column paving strategies, making it both theoretically rigorous and practically efficient.

• Additionally, we devise RPLSS‑CD for inconsistent systems by selecting columns with probabilities proportional to squared norms, accumulating historical updates to overcome slow single‑component convergence.

• Finally, to further balance local coordinate updates and global residual feedback, we developed the RPLSS‑GCD method. RPLSS-GCD achieves a clever trade-off by flexibly combining the current residual with randomly selected standard basis vectors, significantly accelerating the convergence speed of ill-conditioned and inconsistent problems.

Contributions: All authors contributed equally to Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Writing an original draft, Writing, Review, and Editing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement: The data used to support the findings of this study are included within this article.

Conflict of Interest: The authors declare no conflict of interest.

6 Details of \(\bar{\mathbf{p}}_k\)↩︎

Consider the constrained optimization problem \[\min_{\bar{\mathbf{p}}_k\in\mathbb{R}^n}\frac{1}{2}\|\mathbf{B\bar{\mathbf{p}}_k}\|_2^2-\mathbf{r}_{k-1}^T\mathbf{A}^T\bar{\mathbf{p}}_k\] \[subject\;to\quad\bar{\mathbf{S}}_k^T\mathbf{A}^T\mathbf{z}_k=0,\] where \(\mathbf{r}_{k-1}=\mathbf{A}^T\mathbf{z}_{k-1}\), or equivalently, \[\label{0000} \begin{bmatrix} \mathbf{B}^T\mathbf{B} & \mathbf{A}\bar{\mathbf{S}}_k \\ \bar{\mathbf{S}}_k^T\mathbf{A}^T & \mathbf{0}_{k\times k} \end{bmatrix} \begin{bmatrix} \bar{\mathbf{p}}_k \\ \boldsymbol{\lambda}_k \end{bmatrix} = \begin{bmatrix} \mathbf{A} \\ \bar{\mathbf{S}}_k^T \end{bmatrix} \left( \mathbf{A}^T\mathbf{z}_{k-1}\right).\tag{15}\] For the first term of Eq.@eq:0000 , the update \(\bar{\mathbf{p}}_k\) is of the form \[\bar{\mathbf{p}}_k=\mathbf{W}\mathbf{A}\bar{\mathbf{S}}_k(\bar{\mathbf{S}}_k^T\mathbf{A}^T\mathbf{W}\mathbf{A}\bar{\mathbf{s}}_k)^{-1}\bar{\mathbf{S}}_k^T\mathbf{A}^T\mathbf{z}_{k-1}.\] Similar to the update \({\mathbf{p}}_k\) of PLSS in [27], this update can be rewritten as \[\bar{\mathbf{p}}_k=\bar{\mathbf{P}}_{k-1}\bar{\mathbf{g}}_{k-1}+\bar{\mathbf{\gamma}}_{k-1}\bar{\mathbf{y}}_k,\] where \(\bar{\mathbf{y}}_k=\mathbf{A}\bar{\mathbf{s}}_k\). By the orthogonality \(\bar{\mathbf{P}}_{k-1}^T\bar{\mathbf{p}}_{k}=0\), it follows that \[\bar{\mathbf{g}}_{k-1}=-\bar{\gamma}_{k-1}\bar{\mathbf{\Theta}}_{k-1} ^{-1}\bar{\mathbf{P}}_{k-1}^T\bar{\mathbf{y}}_k,\] where \(\bar{\mathbf{\Theta}}_{k-1}=\bar{\mathbf{P}}_{k-1}^T\bar{\mathbf{P}}_{k-1}\). Following the fact \(\bar{\mathbf{s}}_k^T\mathbf{A}\bar{\mathbf{p}}_k=\bar{\mathbf{s}}_k^T\mathbf{r}_{k-1}\), we have \[\bar{\gamma}_{k-1}=\frac{\bar{\mathbf{s}}_k^T\mathbf{A}^T\mathbf{z}_{k-1}}{\|\mathbf{A}\bar{\mathbf{s}}_k\|_2^2-\|\bar{\mathbf{d}}_{k-1}\|_2^2},\] where \(\bar{\mathbf{d}}_{k-1}=\bar{\mathbf{\Theta}}_{k-1} ^{-1/2}\bar{\mathbf{P}}_{k-1}^T\mathbf{A}\bar{\mathbf{s}}_k\).

References↩︎