1 Introduction↩︎

Randomized iterative methods [1]–[4], necoara2019faster?, strohmer2009randomized? form a powerful class of algorithms for solving the large-scale linear system \[\label{main-prob} Ax = b,\;A \in \mathbb{R}^{m \times n},\;b \in \mathbb{R}^m.\tag{1}\] These methods can be regarded as iterative sketching approaches in the randomized numerical linear algebra paradigm [5], [6], and also serve as extensions of stochastic approximation methods from optimization, including stochastic gradient descent (SGD) [1], [7], [8]. However, such methods may converge slowly for ill-conditioned coefficient matrices, and therefore a preconditioner is needed [9]–[14].

The randomized Kaczmarz (RK) method strohmer2009randomized? is one of the most widely used iterative methods for solving large-scale linear systems in recent years due to its simplicity and efficiency. Recently, Lok and Rebrova [15] incorporated subspace constraints into the RK method, developing the subspace-constrained RK (SCRK) method and demonstrating that it achieves superior numerical efficiency over the standard RK method. They further extended this technique to the randomized block Kaczmarz (RBK) method [16], however, under the assumption that the sampled submatrices possess full row rank [15]. In this paper, we extend this subspace-constrained preconditioning technique to a broader class of randomized iterative methods. Our analysis removes the full row rank restriction on sampled submatrices and generalizes the subspace-constrained preconditioner to arbitrary linear systems.

We next briefly illustrate the generalization from the standard RK method to the general randomized iterative method adopted in this work. For any \(i \in [m]:=\{1,\ldots,m\}\), let \(e_i\) denote the \(i\)-th unit vector, \(\top\) denotes the transpose, \(a^\top_{i}\) denote the \(i\)-th row of \(A\), \(b_i\) denote the \(i\)-th entry of \(b\), and let \(\|\cdot\|_2\) and \(\|\cdot\|_F\) denote the Euclidean and Frobenius norms, respectively. The RK method strohmer2009randomized? employs the following update formula \[\label{RK-iteration} x^{k+1} = x^k - \frac{\langle a_{i_k}, x^k\rangle - b_{i_k}}{\|a_{i_k}\|_2^2} a_{i_k}=x^k-\alpha_k A^\top e_{i_k} e_{i_k}^\top (Ax^k-b),\tag{2}\] where the \(i_k\)-th row is selected with probability \(\|a_{i}\|_2^2 / \|A\|_F^2\), and \(\alpha_k = 1 / \|a_{i_k}\|_2^2\) is the stepsize. Noting that \(e_{i_k}\) in 2 acts as the role that extracts partial information of \(A\), we extend \(e_{i_k}\) to a general randomized sketching matrix \(S_k\in \mathbb{R}^{m\times q}\) and apply the following iteration scheme \[\label{RIM-iteration} x^{k+1} =x^k-\alpha_k A^\top S_{k} S_{k}^\top (Ax^k-b).\tag{3}\] The matrix \(S_k\) is sampled from a user-defined probability space \({(\Omega,\mathcal{F},\mathbf{P})}\). This is precisely the randomized iterative method proposed in [2], [17] for solving linear systems. We note that by choosing the probability space \((\Omega, \mathcal{F}, \mathbf{P})\) in 3 appropriately, one can recover existing methods or develop new ones [17], such as the RBK method [16] and the randomized average block Kaczmarz (RABK) method necoara2019faster?.

1.1 Our contributions↩︎

For any index subset \(\mathcal{I}\subseteq[m]\), we use \(A_{\mathcal{I}}\) and \(b_{\mathcal{I}}\) to denote the row submatrix of \(A\) and the subvector of \(b\) indexed by \(\mathcal{I}\), respectively. Let \(\mathcal{I}_p \subset [m]\) be a fixed index set and \(\mathcal{I}_r = [m] \setminus \mathcal{I}_p\). The key idea of the subspace-constrained preconditioner is that we enforce the subsystem \(A_{\mathcal{I}_p}x=b_{\mathcal{I}_p}\) to hold exactly throughout the iteration. That is, every iterate \(x^k\) is required to satisfy \(A_{\mathcal{I}_p}x^k=b_{\mathcal{I}_p}\). This strategy of fixing the row block \(A_{\mathcal{I}_p}\) is particularly beneficial when such a block possesses favorable structural or spectral properties [18].

The main contributions of this work are as follows.

• We introduce a QR-like factorization of the coefficient matrix \(A^\top\) based on the fixed submatrix \(A_{\mathcal{I}_p}\). This factorization allows us not only to transform the original linear system into an equivalent preconditioned system, but also to avoid the full-rank assumptions adopted in [15]. Thus, it further extends the applicability of the subspace-constrained preconditioner to any type of linear system, including underdetermined, overdetermined, full-rank, and rank-deficient matrices.

• We introduce a projected SGD to solve the preconditioned system, where the projection enforces the subsystem \(A_{\mathcal{I}_p}x = b_{\mathcal{I}_p}\) to hold exactly throughout the iteration. We show that, when the initial point is chosen appropriately, the proposed method is equivalent to applying the randomized iterative method 3 with preconditioner \(P := I - A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p}\) to solve the reduced linear system \(A_{\mathcal{I}_r} x = b_{\mathcal{I}_r},\) where \(A_{\mathcal{I}_p}^\dagger\) denotes the pseudoinverse of \(A_{\mathcal{I}_p}\). This equivalence enables parallel implementation and greatly reduces computational cost, while the smaller matrix \(A_{\mathcal{I}_r} P\) enjoys a more favorable singular value distribution than the original coefficient matrix \(A\).

• We prove that the proposed subspace-constrained randomized iterative method (SCRIM) converges linearly in expectation, with a convergence factor that depends on both the choice of the probability space and the reduced matrix \(A_{\mathcal{I}_r} P\). We prove that the subspace-constrained preconditioning yields a tighter convergence factor than its unconstrained counterpart. In addition, we establish a decoupled upper bound on the convergence factor and, based on this bound, provide practical strategies for selecting the constraint subspace.

• We then replace the stochastic gradient direction used by SCRIM with orthogonalized search directions, motivated by the success of the orthogonal direction method [19]. Specifically, each descent direction is constructed by applying a truncated Gram-Schmidt orthogonalization process to the stochastic gradients. We demonstrate that the resulting framework constitutes an iterative-sketching-based Krylov subspace (IS-Krylov) method with subspace-constrained preconditioning, and we prove that the proposed SC-IS-Krylov method achieves a tighter upper bound for convergence. Numerical experiments validate the theoretical findings and demonstrate the efficiency of the proposed method.

1.2 Related work↩︎

1.2.1 The Kaczmarz method↩︎

The literature on solving linear systems via iterative methods is vast and has a long history [20]–[22]. Row-action methods, for instance, have been proposed as practical approaches for solving large-scale linear systems. These methods access and process one or a few rows of \(A\) at a time. A prominent example of the row-action method is the Kaczmarz method [23] and its randomized variant strohmer2009randomized?. In recent years, a large amount of work has been devoted to the development of Kaczmarz-type methods, including block Kaczmarz methods [3], [16], [17], [24], necoara2019faster?, accelerated RK methods [2], [4], [25]–[30], randomized Douglas-Rachford methods [31], [32], and others.

Very recently, Lok and Rebrova [15] introduced the subspace-constrained randomized Kaczmarz (SCRK) method and its block extensions, as briefly reviewed in Remarks 2 and 3. The core idea of SCRK is to enforce the subsystems \(A_{\mathcal{I}_p}x = b_{\mathcal{I}_p}\) to hold exactly at all iterations, such that every generated iterate \(x^k\) strictly satisfies \(A_{\mathcal{I}_p}x^k = b_{\mathcal{I}_p}\). In this paper, we extend this line of research by proposing a subspace-constrained preconditioner for a broad range of randomized iterative methods 3 . Furthermore, based on a QR-like factorization, our method does not require the full-row-rank assumption, which significantly broadens and enriches the results presented in [15].

1.2.2 Preconditioned SGD↩︎

SGD [7], [33]–[35] has recently gained popularity for solving large-scale optimization problems of the generic form \(\min_x \, F(x) := \mathbb{E}[f(x;\xi)],\) where \(\xi\) is sampled from \((\Omega, \mathcal{F}, \mathbf{P})\). However, the convergence of SGD can be slow when the problem is ill-conditioned. To mitigate this issue, one can incorporate a preconditioner that scales the gradient direction, leading to preconditioned SGD \[\label{iteration-pre-SGD} x^{k+1} = x^k - \alpha_k M \nabla f(x^k; \xi_k),\tag{4}\] where \(M\) is a symmetric positive definite preconditioner, \(\alpha_k > 0\) denotes the stepsize, and \(\nabla f(x^k;\xi_k)\) stands for the stochastic gradient evaluated at \(x^k\) based on an i.i.d. sample \(\xi_k\). The standard SGD is recovered when \(M = I\). A substantial body of work has explored preconditioned SGD [36]–[40], with the design of the preconditioner \(M\) broadly falling into two paradigms. The first employs diagonal preconditioners (e.g., Adam [41], AMSGrad [42]), which adaptively rescale coordinates using historical second-moment estimates. The second paradigm constructs curvature-aware preconditioners by approximating second-order structure [43]–[46], such as the Hessian or Fisher information matrix, to improve conditioning.

In many applications, one may require \(x\) to belong to a certain constraint set \(\mathcal{X}\). This motivates the projected variant of SGD \[\label{iteration-pro-SGD} x^{k+1} = \Pi_{\mathcal{X}} \left( x^k - \alpha_k \nabla f(x^k; \xi_k) \right),\tag{5}\] where \(\Pi_{\mathcal{X}}\) denotes the projection onto the constraint set \(\mathcal{X}\). We note that the randomized iterative method 3 can be viewed as SGD for solving the stochastic optimization problem \(\min_x \frac{1}{2} \mathbb{E}\!\left[ \|S^\top (A x - b)\|_2^2 \right]\) reformulated from the linear system. Indeed, our proposed subspace-constrained iterative framework admits a dual interpretation as both a preconditioned SGD scheme and a projected SGD scheme; see Section 3. This unified dual perspective provides a novel paradigm for designing and analyzing preconditioned randomized iterative methods for large-scale linear systems.

1.2.3 Minimal error methods↩︎

Minimal error methods [47]–[49] are iterative methods that compute the next iterate \(x^k\) by minimizing \(\|x - x^*\|_2^2\) over a prescribed search space \(\mathcal{L}_k\), i.e., \[\label{eq-MEM} x^{k+1} = \mathop{\rm argmin}\limits_{x \in \mathcal{L}_k} \|x - x^*\|_2^2,\tag{6}\] where \(x^*\) denotes a certain solution of the linear system \(Ax = b\). This approach directly targets proximity to the true solution, in contrast to residual-minimizing methods such as GMRES [50], which minimize \(\|b - Ax\|_2^2\) instead. When \(A \in \mathbb{R}^{n \times n}\) is nonsingular and the search space is chosen as \(\mathcal{L}_k = x^0 + \mathcal{K}_k(A^\top, A^\top r^0)\) with an initial point \(x^0\) and initial residual \(r^0 = b - A x^0\), the minimal error method 6 reduces to the generalized minimal error (GMERR) method [49], [51]. Here \(\mathcal{K}_k(A^\top, A^\top r^0)\) denotes the Krylov subspace of order \(k\), whose definition can be found in 26 . When \(A \in \mathbb{R}^{m \times n}\) is possibly rectangular, the minimal error method 6 with search space \(\mathcal{L}_k = x^0 + \mathcal{K}_k(A^\top A, A^\top r^0)\) coincides with the conjugate gradient normal equation (CGNE) method [52].

Note that the error-minimizing scheme amounts to projecting the true solution \(x^*\) onto the subspace \(\mathcal{L}_k\) at iteration \(k\). It faces two fundamental challenges: (1) \(x^*\) is unknown in practice and (2) the projection onto \(\mathcal{L}_k\) may be computationally expensive. Classical minimal error methods circumvent the dependence on \(x^*\) by constructing \(\mathcal{L}_k\) as an affine shift of a Krylov subspace, so that the projection can be reformulated entirely in terms of accessible quantities such as residuals and matrix-vector products. To facilitate efficient computation of the projection, these methods build an orthogonal basis for \(\mathcal{L}_k\). A canonical realization of this strategy is the orthogonal direction method [47], [49], which generates a sequence of mutually orthogonal search directions \(p^0, p^1, \dots, p^k\) and updates the iterate by minimizing \(\|x - x^*\|_2^2\) along each new direction; see Section 4.

Inspired by the intrinsic error-minimization principle and orthogonal-direction acceleration mechanism of classical minimal error methods, we further enhance the convergence performance of our subspace-constrained iterative framework. This algorithmic improvement ultimately gives rise to the proposed subspace-constrained iterative-sketching-based Krylov subspace (SC-IS-Krylov) method.

1.3 Organization↩︎

The remainder of this paper is organized as follows. In Section 2, we introduce the notation and preliminary results used throughout the paper. Section 3 develops the subspace-constrained randomized iterative method, and establishes its linear convergence in expectation. Section 4 proposes subspace-constrained iterative-sketching-based Krylov methods and analyzes their convergence properties. Numerical experiments demonstrating the effectiveness of the proposed methods are presented in Section 5. Section 6 concludes the paper. Proofs of the main results involving lengthy derivations are provided in the appendices.

2 Notation and preliminaries↩︎

2.1 Notations↩︎

For vector \(x\in\mathbb{R}^n\), we use \(x_i,x^\top\), and \(\|x\|_2\) to denote the \(i\)-th entry, the transpose, and the Euclidean norm of \(x\), respectively. We denote the linear span of vectors \(x^1,\ldots, x^k \in \mathbb{R}^n\) as \(\operatorname{span}\{ x^1,\ldots, x^k \}\). For any matrix \(A \in \mathbb{R}^{m \times n}\), we use \(a^\top_i\), \(A^\dagger\), \(A_{\mathcal{J}}\), \(A_{\mathcal{J}}^\dagger\), \(A^\top\), \(\|A\|_2\), \(\|A\|_F\), \(\operatorname{rank}(A)\), and \(\operatorname{Range}(A)\) to denote the \(i\)-th row, the Moore-Penrose pseudoinverse, the row submatrix indexed by \(\mathcal{J}\), the pseudoinverse of the row submatrix \(A_{\mathcal{J}}\), the transpose, the spectral norm, the Frobenius norm, the rank, and the column space, respectively. We use \(\lambda_{\min}(A)\) to denote the smallest eigenvalue of a symmetric matrix \(A\), and \(\sigma_{\min}(A)\) to denote the smallest nonzero singular value of a general matrix \(A\). The cardinality of the set \(\mathcal{J}\) is denoted by \(\vert \mathcal{J}\vert\). For any random variables \(\xi\), we use \(\mathbb{E}[\xi]\) to denote the expectation of \(\xi\).

2.2 A QR-like factorization↩︎

In this subsection, we introduce a QR-like factorization for \(A^\top\). We note that this factorization is introduced solely for the purpose of theoretical analysis and algorithm development, it is not explicitly computed in the actual implementation of our method.

Let \(\mathcal{I}_p \subset [m]\) be a fixed index set with cardinality \(m_p = |\mathcal{I}_p|\), and define \(\mathcal{I}_r = [m] \setminus \mathcal{I}_p\) with \(m_r= |\mathcal{I}_r|\). Without loss of generality, we assume \(\mathcal{I}_p =[m_p]\). Let \(P := I - A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p}\). Using the identity \(A_{\mathcal{I}_p}^\top(A_{\mathcal{I}_p}^\dagger)^\top A_{\mathcal{I}_r}^\top+PA_{\mathcal{I}_r}^\top=A_{\mathcal{I}_r}^\top,\) we obtain \[\label{eq-Atop-QR-like} A^\top=\begin{pmatrix} A_{\mathcal{I}_p}^\top & A_{\mathcal{I}_r}^\top \end{pmatrix} =\underbrace{\begin{pmatrix} A_{\mathcal{I}_p}^\top & PA_{\mathcal{I}_r}^\top \end{pmatrix}}_{\hat{A}^\top} \underbrace{\begin{pmatrix} I & (A_{\mathcal{I}_p}^\dagger)^\top A_{\mathcal{I}_r}^\top\\ 0 & I \end{pmatrix}}_{L^\top} =\hat{A}^\top L^\top.\tag{7}\] Substituting \(P=I-A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p}\) and employing the property \(A_{\mathcal{I}_p}A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p}=A_{\mathcal{I}_p}\), we derive \[A_{\mathcal{I}_p}(A_{\mathcal{I}_r}P)^\top = A_{\mathcal{I}_p} P A_{\mathcal{I}_r}^\top = \bigl(A_{\mathcal{I}_p}-A_{\mathcal{I}_p}A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p}\bigr)A_{\mathcal{I}_r}^\top =0,\] which verifies that \(\hat{A}^\top\) has orthogonal column blocks. Hence, the factorization \(A^\top = \hat{A}^\top L^\top\) in 7 decomposes \(A^\top\) into the product of a matrix \(\hat{A}^\top\) with orthogonal column blocks and a unit upper block triangular matrix \(L^\top\). This yields a QR-like factorization [53] for \(A^\top\), where only the block columns of \(\hat{A}^\top\) satisfy the orthogonality condition.

Now, the linear system 1 can be rewritten as the following preconditioned system \[\label{main-pre-prob} \hat{A} x = \hat{b} \;\;\text{with} \;\;\hat{b} := L^{-1} b = \begin{pmatrix} b_{\mathcal{I}_p} \\ b_{\mathcal{I}_r} - A_{\mathcal{I}_r} A_{\mathcal{I}_p}^\dagger b_{\mathcal{I}_p} \end{pmatrix}.\tag{8}\] To solve 8 , we note that the proposed method in Section 3 (see 12 ) is equivalent to applying a randomized iterative method 3 to the reduced linear system \[A_{\mathcal{I}_r} P x = \hat{b}_{\mathcal{I}_r}.\] We next present a result characterizing the singular value distribution of \(A_{\mathcal{I}_r} P\) in comparison with the original coefficient matrix \(A\).

In general, \(A^\dagger\neq \begin{pmatrix} A_{\mathcal{I}_p}^\dagger & A_{\mathcal{I}_r}^\dagger \end{pmatrix}\). However, if \(A_{\mathcal{I}_p} A_{\mathcal{I}_r}^\top=0\), the following lemma shows that \(A^\dagger= \begin{pmatrix} A_{\mathcal{I}_p}^\dagger & A_{\mathcal{I}_r}^\dagger \end{pmatrix}\). Lemma 2 is useful in our argument and we believe it is of independent interest.

The following result shows that the original linear system \(Ax=b\) and the preconditioned system \(\hat{A} x = \hat{b}\) share the same unique least-squares solution.

2.3 Preliminaries on the probability space↩︎

First, we state a basic assumption on the probability space \((\Omega, \mathcal{F}, \mathbf{P})\) used throughout this paper.

The following lemmas are crucial for our convergence analysis.

We summarize Lemma 2.3 in [54] and Lemma 2.5 in [2] in the following lemma.

3 Randomized iterative methods with subspace-constrained preconditioning↩︎

In this section, we develop subspace-constrained randomized iterative methods (SCRIM) for solving the linear system 1 . The SCRIM method can be viewed as a projected SGD method 5 for solving the following constrained stochastic optimization problem reformulated from the preconditioned linear system 8 \[\label{main-sto-prob} \begin{align} \min_{x\in\mathbb{R}^n}&\quad f(x)=\mathbb{E}[f_{S}(x)]\\ \text{subject to}&\quad x\in \mathcal{X}_p:=\{x\in\mathbb{R}^n \mid A_{\mathcal{I}_p}x=b_{\mathcal{I}_p}\}, \end{align}\tag{9}\] where \(f_{S}(x):=\frac{1}{2}\lVert S^\top(A_{\mathcal{I}_r}Px-\hat{b}_{\mathcal{I}_r})\rVert_2^2\) with \(S\) being a random sketching matrix drawn from \((\Omega,\mathcal{F},\mathbf{P})\), and \(\hat{b}_{\mathcal{I}_r}=b_{\mathcal{I}_r} - A_{\mathcal{I}_r} A_{\mathcal{I}_p}^\dagger b_{\mathcal{I}_p}\). The resulting projected SGD iteration scheme reads \[\label{eq-PSGD} \left\{ \begin{array}{ll} y^{k+1}=x^k-\alpha_kPA_{\mathcal{I}_r}^\top S_{k}S_{k}^\top(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r}),\\ x^{k+1}=\Pi_{\mathcal{X}_p}(y^{k+1})=y^{k+1}-A_{\mathcal{I}_p}^\dagger \big(A_{\mathcal{I}_p}y^{k+1}-b_{\mathcal{I}_p}\big), \end{array} \right.\tag{10}\] where \(\alpha_k\) denotes the stepsize and \(\Pi_{\mathcal{X}_p}(\cdot)\) represents the orthogonal projection onto the feasible set \(\mathcal{X}_p\). Substituting the intermediate update \(y^{k+1}\) into the projection step and note that \(A_{\mathcal{I}_p}P = 0\), yields the compact iteration \[\label{eq-bas-xk431} \begin{align} x^{k+1}=x^k-A_{\mathcal{I}_p}^\dagger\big(A_{\mathcal{I}_p}x^k-b_{\mathcal{I}_p}\big)-\alpha_k(A_{\mathcal{I}_r}P)^\top S_{k}S_{k}^\top\big(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r}\big). \end{align}\tag{11}\] If we choose the initial point \(x^0 = A_{\mathcal{I}_p}^\dagger b_{\mathcal{I}_p} \in \mathcal{X}_p\), then all subsequent iterates satisfy \(x^k \in \mathcal{X}_p\). Under this valid initialization, the update 11 simplifies to \[\label{eq-bas-xk431-2} x^{k+1}=x^k-\alpha_kPA_{\mathcal{I}_r}^\top S_{k}S_{k}^\top\left(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r}\right),\tag{12}\] which can be viewed as the randomized iterative method 3 for solving \(A_{\mathcal{I}_r}Px=\hat{b}_{\mathcal{I}_r}\).

We further show that the auxiliary reformulation involving \(\hat{b}_{\mathcal{I}_r}\) in 12 is unnecessary. For any iterate \(x^k\), we have \(Px^k=x^k-A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p}x^k=x^k-A_{\mathcal{I}_p}^\dagger b_{\mathcal{I}_p}\), which yields \(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r}=A_{\mathcal{I}_r}x^k-A_{\mathcal{I}_r}A_{\mathcal{I}_p}^\dagger b_{\mathcal{I}_p}-\hat{b}_{\mathcal{I}_r}=A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r}\). Based on this equality, 12 can be further simplified as \[\label{iter-sim} x^{k+1}=x^k-\alpha_kPA_{\mathcal{I}_r}^\top S_{k}S_{k}^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r}),\;k\geq0.\tag{13}\] This iteration scheme can be viewed as a preconditioned SGD 4 (or preconditioned randomized iterative method) with the preconditioner matrix \(P\) for the linear system \(A_{\mathcal{I}_r}x=b_{\mathcal{I}_r}\).

Note that when \(S_{k}^\top(A_{\mathcal{I}_r}x^k - b_{\mathcal{I}_r}) = 0\), the iteration scheme 13 reduces to \(x^{k+1}=x^k\), which corresponds to a null step. To avoid such null steps, we construct the randomized sketching matrix \(S_k\) to satisfy \(S_{k}^\top(A_{\mathcal{I}_r}x^k - b_{\mathcal{I}_r}) \neq 0\). In addition, it follows from Lemma 4 that \(PA_{\mathcal{I}_r}^\top S_{k}S_{k}^\top(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r})\neq0\) holds if and only if \(S_{k}^\top(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r})\neq0\). Combining this result with the identity \(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r}=A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r}\), we conclude that \(PA_{\mathcal{I}_r}^\top S_{k}S_{k}^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\neq0\) if and only if \(S_{k}^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\neq0\). Hence, \[\label{eq-Ladap} \begin{align} L_{\text{adap}}^k=\frac{\lVert S_{k}^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\rVert_2^2}{\lVert PA_{\mathcal{I}_r}^\top S_{k}S_{k}^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\rVert_2^2} \end{align}\tag{14}\] is well-defined throughout the iteration process. Now we are ready to state the proposed SCRIM, which is formally described in Algorithm 1.

3.1 Convergence analysis↩︎

For convenience, we define \[\label{eq-rho} \rho:=1- \zeta(2-\zeta)\sigma^2_{\min}(H^{\frac{1}{2}}A_{\mathcal{I}_r}P),\tag{15}\] where \(H\) is given by ?? . We have the following convergence result for Algorithm 1.

Next, we show that the convergence factor in Theorem 1 is at least as tight as that of the randomized iterative method in 3 (i.e., \(\mathcal{I}_p = \emptyset\)). To facilitate a fair comparison, we impose a unified formulation for the randomized sketching matrix. Specifically, let \(\tilde{S} \in \mathbb{R}^{m \times q}\) be drawn from the probability space \((\tilde{\Omega}, \tilde{\mathcal{F}}, \tilde{\mathbf{P}})\). For convenience, we assume \(\mathcal{I}_p = \{1, \ldots, m_p\}\), such that \(\tilde{S}\) admits the following block row partition \[\label{eq-S} \tilde{S} = \begin{pmatrix} \tilde{S}_{\mathcal{I}_p} \\ S \end{pmatrix}, \quad \tilde{S}_{\mathcal{I}_p} \in \mathbb{R}^{m_p \times q},\quad S \in \mathbb{R}^{m_r \times q}, \quad m = m_p + m_r.\tag{16}\] In this setting, \(S\) corresponds to the randomized sketching matrix employed in Algorithm 1. Now the randomized iterative method in 3 exhibits the following convergence property: \[\mathbb{E}\left[\lVert x^{k} - A^\dagger b \rVert_2^2\right] \leq \tilde{\rho}^{k} \lVert x^0 - A^\dagger b \rVert_2^2,\] where \[\label{eq-tH-trho} \tilde{H} := \mathbb{E} \left[ \frac{\tilde{S} \tilde{S}^\top}{\lVert \tilde{S}^\top A \rVert_2^2} \right] \quad \text{and} \quad \tilde{\rho} = 1 - \zeta(2-\zeta) \sigma_{\min}^2(\tilde{H}^{1/2} A).\tag{17}\] The following result verifies that Algorithm 1 achieves a tighter convergence factor than the randomized iterative method.

3.2 Selection strategies for \(\mathcal{I}_p\)↩︎

Ideally, the index set \(\mathcal{I}_p\) is chosen to minimize the convergence factor \(\rho\), which is equivalent to maximizing \(\sigma_{\min}(H^{\frac{1}{2}}A_{\mathcal{I}_r}P)\). However, since both \(P\) and \(H\) depend on \(\mathcal{I}_p\), maximizing \(\sigma_{\min}(H^{\frac{1}{2}}A_{\mathcal{I}_r}P)\) with respect to \(\mathcal{I}_p\) admits no closed-form solution and remains computationally intractable in practice. The following result establishes a tractable lower bound for \(\sigma_{\min}(H^{\frac{1}{2}}A_{\mathcal{I}_r}P)\), which provides a feasible surrogate criterion for selecting the index set \(\mathcal{I}_p\).

For any matrix \(B\), we define \(\kappa(B):=\frac{\lVert B\rVert_F}{\sigma_{\min}(B)}\) as the scaled condition number [55], [56] of \(B\). Since \(\bar{H}\) is independent of \(\mathcal{I}_p\), the only term in the lower bound ?? that depends on \(\mathcal{I}_p\) is \(\kappa(A_{\mathcal{I}_r}P)\). A natural strategy for selecting the index set \(\mathcal{I}_p\) is to minimize the scaled condition number \(\kappa(A_{\mathcal{I}_r}P)\), yet this optimization task remains computationally challenging [57]. We therefore adopt the minimization of \(\|A_{\mathcal{I}_r}P\|_F\) as a computationally feasible surrogate criterion for choosing \(\mathcal{I}_p\). Recalling that \(A_{\mathcal{I}_p}P=0\), we obtain \(AP= \begin{pmatrix} 0\\ A_{\mathcal{I}_r}P \end{pmatrix},\) which implies \(\lVert A_{\mathcal{I}_r}P\rVert_F^2=\lVert AP\rVert_F^2= \| A - A A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p} \|^2_F\). For any given \(m_p\), we thus consider the problem \[\min_{\substack{\mathcal{I}_p \subset [m],|\mathcal{I}_p| = m_p}} \left\| A - A A_{\mathcal{I}_p}^\dagger A_{\mathcal{I}_p} \right\|^2_F,\] which exactly corresponds to the row interpolative decomposition (ID) problem [58]–[61]. In the following, we introduce several standard ID-based row selection strategies for constructing the index set \(\mathcal{I}_p\).

Column-pivoted QR: Column-pivoted QR (CPQR) [53] is a deterministic greedy method for skeleton selection. Let \(X^{(0)}=A\), and let \(X_i^{(t)}\) denote the \(i\)-th row of \(X^{(t)}\). At the \((t+1)\)-th step, CPQR performs \[\label{eq-CPQR} s_{t+1} \in \arg\max_{i\in[m]} \left\|X^{(t)}_i\right\|_2^2,\quad X^{(t+1)} = X^{(t)} - X^{(t)} \frac{ \left(X^{(t)}_{s_{t+1}}\right)^\top X^{(t)}_{s_{t+1}} }{ \left\|X^{(t)}_{s_{t+1}}\right\|_2^2 }.\tag{18}\] After \(m_p\) steps, the selected pivots are collected as \(\mathcal{I}_p=\{s_1,\ldots,s_{m_p}\}.\)

SVD-based selection: The SVD-based selection [53] is motivated by the ideal case in which the selected rows span the leading right singular subspace. Let \(A=U\Sigma V^\top\) be the singular value decomposition of \(A\), and let \(V_{m_p}\) denote the leading \(m_p\) right singular vectors of \(A\). If \(\operatorname{Range}(A_{\mathcal{I}_p}^{\top}) = \operatorname{Range}(V_{m_p}),\) then \(A_{\mathcal{I}_p}^{\dagger}A_{\mathcal{I}_p}=V_{m_p}V_{m_p}^{\top}\), and hence \[\left\|A-AA_{\mathcal{I}_p}^{\dagger}A_{\mathcal{I}_p}\right\|_F^2 = \left\|A-AV_{m_p}V_{m_p}^{\top}\right\|_F^2 = \sum_{i=m_p+1}^{\min(m,n)}\sigma_i^2(A).\] In this case, the ID approximation error \(\left\|A-AA_{\mathcal{I}_p}^{\dagger}A_{\mathcal{I}_p}\right\|_F\) coincides with the error of the best rank-\(m_p\) approximation of \(A\) in the Frobenius norm [62]. To approximate such an index set, we form \(AV_{m_p}\), whose rows contain the coordinates of the rows of \(A\) along the leading right singular directions. We then apply the row-selection CPQR procedure to \(AV_{m_p}\) and take the first \(m_p\) row pivots as \(\mathcal{I}_p\). This SVD-based selection procedure serves as an idealized benchmark, but is computationally expensive because it requires computing the leading right singular vectors.

Squared-norm sampling: Squared-norm sampling (SqNorm) [63], [64] selects rows randomly according to their relative importance, measured here by the squared row norms \(\frac{\|a_i\|_2^2}{\|A\|_F^2}\), for \(i\in[m]\). We sample \(m_p\) rows without replacement from this distribution. This method is simple to implement and has low row-selection cost.

Sketched CPQR: Sketched CPQR (SkCPQR) [58], [65] is a randomized variant of CPQR that performs pivot selection on a compressed version of the matrix. Specifically, we draw a Gaussian matrix \(G\in\mathbb{R}^{n\times s}\) with \(s< n\), and form \(Y=AG\in\mathbb{R}^{m\times s}\). The rows of \(Y\) provide a lower-dimensional representation of the rows of \(A\). We then apply the row-selection CPQR procedure to \(Y\), and use the resulting row pivots as \(\mathcal{I}_p\). This strategy performs pivot selection in a lower-dimensional sketched space, thereby reducing the cost of CPQR. Related randomized sketching methods for ID can be found in [60], [61], [66], [67].

Robust blockwise random pivoting: Random pivoting [67]–[69] is an adaptive alternative to deterministic greedy pivoting. Instead of selecting the largest residual-norm pivot at each step in 18 , robust blockwise random pivoting (RBRP) [59] samples multiple pivots according to the squared row norms of the current residual matrix. This preserves the adaptive nature of greedy pivoting while introducing randomness into the pivot selection.

4 Subspace-constrained iterative-sketching-based Krylov subspace methods↩︎

In this section, we replace the stochastic gradient direction used in Algorithm 1 with search directions obtained by orthogonalizing the negative stochastic gradients via a truncated Gram-Schmidt process, motivated by the success of the orthogonal direction method [19]. This leads to the development of subspace-constrained iterative-sketching-based Krylov subspace (SC-IS-Krylov) method.

We first revisit the core framework of the orthogonal direction method [19]. Assume that we have a set of orthogonal search directions \(p^0, p^1, \ldots, p^k\). The method updates the iterate at the \(k\)-th iteration via \[\label{iter-xk431-cg} x^{k+1} = \mathop{\rm argmin}\limits_{x \in x^k + \operatorname{span}\{p^k\}} \frac{1}{2}\|x - A^\dagger b\|_2^2.\tag{19}\] When \(p^k \neq 0\), this leads to the update \(x^{k+1} = x^k + \delta_k p^k\) with \[\delta_k=\mathop{\rm argmin}\limits_{\delta\in\mathbb{R}}\frac{1}{2}\lVert x^k+\delta p^k-A^\dagger b\rVert_2^2=-\frac{\langle p^k, x^k-A^\dagger b\rangle}{\lVert p^k \rVert_2^2}.\] However, computing \(\delta_k\) requires knowledge of the unknown solution \(A^\dagger b\); therefore, the method is not implementable in general. To eliminate the dependence on \(A^\dagger b\), the author [19] required that \(p^k\) be \(A^\top A\)-conjugate to \(x^{k+1} - A^\dagger b\).

Nevertheless, we address this issue by constructing feasible search directions \(p^k\) based on the negative stochastic gradient \(d^{k}:=-(A_{\mathcal{I}_r}P)^\top S_{k}S_{k}^\top(A_{\mathcal{I}_r}Px^k-\hat{b}_{\mathcal{I}_r})\) to obtain the search directions \(p^k\). In particular, we obtain \(p^k\) via the following Gram-Schmidt orthogonalization process \[\label{iter-pk} \begin{align} p^{k}=d^k-\sum_{i=j_{k,\ell}}^{k-1}\frac{\langle d^k,p^i\rangle}{\lVert p^i\rVert_2^2}p^i, \end{align}\tag{20}\] where \(j_{k,\ell} := \max\{k - \ell + 1, 0\}\) tracks the most recent \(\ell\) iterates and \(\ell\) is a positive integer. We refer to \(\ell\) as the truncation parameter. For instance, \(\ell = +\infty\) corresponds to using all previous vectors \(p^i\). The following lemma shows that \(\langle p^k,x^k-A^\dagger b\rangle=-\lVert S_k^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\rVert_2^2\) under a mild condition.

Hence, if \(S_i^\top(A_{\mathcal{I}_r}x^i-b_{\mathcal{I}_r})\neq0,i=0,\ldots,k\), we have \[\label{def-delta} \delta_k =\frac{\lVert S_{k}^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\rVert_2^2}{\lVert p^k\rVert_2^2},\tag{21}\] where the last equation follows from Lemma 8. The proposed method is summarized in Algorithm 2. The justification for the name SC-IS-Krylov will be provided in Section 4.2.

4.1 Convergence analysis↩︎

We begin by introducing some auxiliary variables. Define \(\tilde{P}_0\) as the empty matrix, and for \(k \geq 1\), define \[\label{eq-tilde-Pk} \tilde{P}_k := \left( p^{j_{k,\ell}},\, p^{j_{k,\ell}+1},\, \ldots,\, p^{k-1} \right) \in \mathbb{R}^{n \times (k - j_{k,\ell})},\tag{22}\] where the vectors \(p^i\) for \(i = j_{k,\ell}, \dots, k-1\) are defined in Algorithm 2. Since the columns of \(\tilde{P}_k\) are mutually orthogonal and nonzero, it follows from Lemma 8 that \[\tilde{P}_k^\dagger = \begin{pmatrix} \frac{p^{j_{k,\ell}}}{\lVert p^{j_{k,\ell}} \rVert_2^2}, \frac{p^{j_{k,\ell}+1}}{\lVert p^{j_{k,\ell}+1} \rVert_2^2}, \ldots, \frac{p^{k-1}}{\lVert p^{k-1} \rVert_2^2} \end{pmatrix}^\top.\] Then, for \(k\geq1\), we can rewrite \(p^k\) in 20 as \[\label{iter-pk-2} p^k=d^k-\tilde{P}_k\tilde{P}_k^\dagger d^k=(I-\tilde{P}_k\tilde{P}_k^\dagger) d^k.\tag{23}\] In addition, we define \[\label{eq-Q95k} \mathcal{Q}_k := \left\{ S \in \Omega \mid S^\top (A_{\mathcal{I}_r} x^k - b_{\mathcal{I}_r}) \neq 0 \right\}.\tag{24}\] and \[\label{eq-q95k} q_k := \inf_{S_k \in \mathcal{Q}_k} \left\{ \left(1 - \frac{\lVert \tilde{P}_k \tilde{P}_k^\dagger d^k \rVert_2^2}{\lVert d^k \rVert_2^2} \right)^{-1} \right\}.\tag{25}\]

At the \(k\)-th iteration, we consider the product probability space \(\left( \prod_{i=0}^k \Omega, \otimes_{i=0}^k \mathcal{F}, \tilde{\mathbf{P}} \right),\) where \(\otimes\) denotes the product of \(\sigma\)-algebras and \(\tilde{\mathbf{P}}\) is the corresponding product measure, as defined in [70]. Let \(\mathcal{B}_k:= \sigma(S_0, S_1, \ldots, S_{k-1})\) be the \(\sigma\)-algebra generated by the random variables \((S_0, S_1, \ldots, S_{k-1})\), where \(\mathcal{B}_0\) is the trivial \(\sigma\)-algebra (i.e., \(\mathcal{B}_0 = \{\emptyset, \Omega\}\)). The conditional expectation given \(\mathcal{B}_k\) is denoted by \(\mathbb{E}_k[\cdot] := \mathbb{E}[\cdot \mid \mathcal{B}_k].\) The convergence result for Algorithm 2 is as follows.

Upon comparing Theorems 1 and 2, we observe that the SC-IS-Krylov method exhibits a convergence bound that is at least as tight as that of the SCRIM method. In particular, for certain probability spaces, the parameter \(q_k\) in Theorem 2 can be strictly greater than \(1\). Indeed, consider the case where the sample space is \(\Omega = \{S\}\) with \(S\) being a fixed matrix such that \(S^\top S\) is positive definite. We next show that when \(A_{\mathcal{I}_r} x^k - b_{\mathcal{I}_r} \neq 0\), it holds that \(q_k > 1\). From the definition of \(q_k\) in 25 and noting that \(\tilde{P}_k^\top = \tilde{P}_k^\top \tilde{P}_k \tilde{P}_k^\dagger\), to ensure \(\tilde{P}_k \tilde{P}_k^\dagger d^k \neq 0\), it suffices to prove that \(\tilde{P}_k^\top d^k \neq 0\). In particular, we show that the last component of \(\tilde{P}_k^\top d^k\) is nonzero, i.e., \[\langle p^{k-1}, d^k \rangle = \frac{\langle x^k - x^{k-1}, d^k \rangle}{\delta_{k-1}} = \frac{\langle x^k - A^\dagger b, d^k \rangle}{\delta_{k-1}} - \frac{\langle x^{k-1} - A^\dagger b, d^k \rangle}{\delta_{k-1}} \neq 0.\] To show this, we first evaluate the term \(\langle x^{k-1} - A^\dagger b,d^k\rangle\). Since \(d^k=-PA_{\mathcal{I}_r}^\top SS^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\) and \(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r}=A_{\mathcal{I}_r}(x^k-x^{k-1})+A_{\mathcal{I}_r}x^{k-1}-b_{\mathcal{I}_r}\), it follows that \[\begin{align} &\langle x^{k-1} - A^\dagger b, d^k \rangle\\ &= -\langle x^{k-1} - A^\dagger b, PA_{\mathcal{I}_r}^\top SS^\top A_{\mathcal{I}_r} (x^k - x^{k-1}) \rangle-\langle x^{k-1} - A^\dagger b, PA_{\mathcal{I}_r}^\top SS^\top (A_{\mathcal{I}_r} x^{k-1} - b_{\mathcal{I}_r}) \rangle\\ &=-\langle PA_{\mathcal{I}_r}^\top SS^\top A_{\mathcal{I}_r} (x^{k-1} - A^\dagger b), x^k - x^{k-1} \rangle-\| S^\top (A_{\mathcal{I}_r} x^{k-1} - b_{\mathcal{I}_r}) \|_2^2\\ &= \langle d_{k-1}, \delta_{k-1} p^{k-1} \rangle-\| S^\top (A_{\mathcal{I}_r} x^{k-1} - b_{\mathcal{I}_r}) \|_2^2. \end{align}\] From the definition of \(p^{k-1}\) in 20 , we have \[\langle d_{k-1}, \delta_{k-1} p^{k-1} \rangle=\delta_{k-1} \langle p^{k-1} + \sum_{i = j_{k-1}}^{k-2} \eta_{k-1}^i p^i, p^{k-1} \rangle= \delta_{k-1} \| p^{k-1} \|_2^2=\| S^\top (A_{\mathcal{I}_r} x^{k-1} - b_{\mathcal{I}_r}) \|_2^2,\] where the last equality follows from the definition of \(\delta_{k-1}\) in 21 . Combining these results, we obtain \(\langle x^{k-1} - A^\dagger b, d^k \rangle = 0\). Thus, we have \[\langle p^{k-1},d^k\rangle=\frac{\langle x^k-x^{k-1},d^k\rangle}{\delta_{k-1}}=\frac{\langle x^k - A^\dagger b,d^k\rangle}{\delta_{k-1}}=\frac{\| S^\top (A_{\mathcal{I}_r} x^k - b_{\mathcal{I}_r}) \|_2^2}{\delta_{k-1}}\neq0.\] This implies that \(\tilde{P}_k^\top d^k \neq 0\) as long as \(A_{\mathcal{I}_r} x^k \neq b_{\mathcal{I}_r}\). Consequently, we have \(q_k>1\).

4.2 Connection to the Krylov subspace method↩︎

We first introduce an equivalent form of the iterate \(x^{k+1}\) obtained by 19 and 20 .

Recall that, given a matrix \(B \in \mathbb{R}^{n \times n}\) and a vector \(r \in \mathbb{R}^n\), the Krylov subspace of order \(k\) is defined as [53] \[\label{def-krylov} \mathcal{K}_k(B, r) := \operatorname{span}\{r, Br, B^2 r, \ldots, B^{k-1} r\}.\tag{26}\] The following theorem demonstrates that \(x^0 + \mathcal{K}_k\bigl(\hat{A}^\top \hat{A},\; \hat{A}^\top (\hat{A} x^0 - \hat{b})\bigr)\) is a specific instance of the subspace in Lemma 9.

By Lemma 9 and Theorem 3, when \(\ell = \infty\) and \(\Omega = \{I\}\), the proposed orthogonalized search direction scheme 19 –20 reduces to a standard Krylov subspace method. Specifically, the update rule 19 takes the form \[\begin{align} x^{k+1} = \mathop{\arg\min}_{x \in \mathbb{R}^n} &\; \lVert x - A^\dagger b \rVert_2^2 \\ \text{subject to} \quad &x \in x^0 + \mathcal{K}_{k+1}\big(\hat{A}^\top \hat{A}, \hat{A}^\top \hat{r}^0\big). \end{align}\] This indicates that for general truncation parameters \(\ell\) and probability spaces \((\Omega,\mathcal{F},\mathbf{P})\), our orthogonalized search direction framework 19 –20 inspires the construction of the novel subspace-constrained iterative-sketching-based Krylov (SC-IS-Krylov) subspace method.

5 Numerical experiments↩︎

In this section, we implement SC-IS-Krylov (Algorithm 2). The case \(\ell=1\) corresponds to SCRIM (Algorithm 1) with \(\zeta=1\). We examine the five strategies discussed in Section 3.2 for selecting the constraint set \(\mathcal{I}_p\). We also compare the SC-IS-Krylov method with the IS-Krylov method proposed in [4]. All experiments are implemented in MATLAB R2024a under Windows 11 on a laptop computer equipped with an Intel Core Ultra 7 155H CPU and 32 GB memory. The code to reproduce our results can be found at https://github.com/xiejx-math/SCRIM.

5.1 Numerical setup↩︎

We consider two types of coefficient matrices in our experiments. The first type consists of synthetic Gaussian matrices with outlying singular values. Given parameters \(m,n\) and the target rank \(r\), we construct \(A=UDV^\top,\) where \(U\in\mathbb{R}^{m\times r}\) and \(V\in\mathbb{R}^{n\times r}\) have orthonormal columns, and \(D=\operatorname{diag}(\sigma_1,\ldots,\sigma_r)\) contains the prescribed nonzero singular values. In MATLAB notation, the matrices \(U\) and \(V\) are generated by \(\texttt{[U,}\sim\texttt{]=qr(randn(m,r),0)}\), and \(\texttt{[V,}\sim\texttt{]=qr(randn(n,r),0)}.\) The singular values are generated from three separated clusters. Let \(n_L\) and \(n_S\) denote the numbers of large and small outlying singular values, respectively, with \(n_L+n_S<r\). The remaining \(r-n_L-n_S\) singular values form the middle cluster. Specifically, the small, middle, and large singular values are sampled from three separated intervals \[R_S=[\beta_S,\gamma_S],\qquad R_M=[\beta_M,\gamma_M],\qquad R_L=[\beta_L,\gamma_L],\] respectively, where \(0<\beta_S\), \(\gamma_S<\beta_M\), and \(\gamma_M<\beta_L\). The sampled singular values are then sorted in nonincreasing order and used as the diagonal entries of \(D\). Given a parameter \(\kappa_M>1\), the middle interval is chosen such that \(\gamma_M/\beta_M\leq \kappa_M,\) so that the condition number of the middle singular-value cluster is bounded by \(\kappa_M\). We refer to such a matrix as an \((m,n,r,n_L,n_S,\kappa_M)\) matrix. When \(m=n=r\), this construction coincides with the \((n,n_L,n_S,\kappa_M)\) matrix model [71]. The second type consists of real-world matrices from the SuiteSparse Matrix Collection kolodziej2019suitesparse? and LIBSVM chang2011libsvm?. For these datasets, we use only the coefficient matrix \(A\).

We construct a consistent linear system by generating the ground-truth solution as \(x^\ast=\texttt{randn(n,1)}\) and setting \(b=Ax^\ast\). We initialize SC-IS-Krylov with \(x^0=A_{\mathcal{I}_p}^{\dagger}b_{\mathcal{I}_p}\), and initialize the IS-Krylov method with \(x^0=0\). The iterative process is terminated when the squared relative solution error \(\operatorname{RSE} = \frac{\|x^k-A^\dagger b\|_2^2}{\|A^\dagger b\|_2^2}\) falls below a prescribed tolerance. In the following experiments, we employ the partition sampling [2], [4], [17]. Given the block size \(q\), we consider a random partition of \([m_r]\). Let \(\varpi\) be a uniformly random permutation of \([m_r]\), and define \[\begin{align} &\mathcal{I}_i = \{\varpi(k): k=(i-1)q+1,(i-1)q+2,\ldots,iq\}, \qquad i=1,2,\ldots,t-1,\\ &\mathcal{I}_t = \{\varpi(k): k=(t-1)q+1,(t-1)q+2,\ldots,m_r\}, \qquad |\mathcal{I}_t|\leq q. \end{align}\] We select the sampling matrix \(S^\top = I_{\mathcal{I}_i}\) with probability \(\|(A_{\mathcal{I}_r})_{\mathcal{I}_i}\|_F^2/\|A_{\mathcal{I}_r}\|_F^2.\) We use SC-IS-Krylov-PS and IS-Krylov-PS to denote the SC-IS-Krylov and IS-Krylov methods with partition sampling, respectively. For a more comprehensive discussion of additional sketching strategies, we refer the reader to, e.g., [6]. In practice, the quantity \(\|S_k^\top (A_{\mathcal{I}_r} x^k - b_{\mathcal{I}_r})\|_2\) is considered zero if it is smaller than eps. Each experiment is repeated over \(20\) independent trials to ensure statistical reliability. The lightly shaded region indicates the range from the minimum to the maximum values, while the darker shaded region corresponds to the interval between the 25-th and 75-th quantiles.

5.2 Comparison of selection strategies for \(\mathcal{I}_p\)↩︎

Recent works [71]–[73] indicate that the convergence behavior of iterative methods for linear systems is closely related to the distribution of singular values, rather than to the classical condition number alone. Motivated by this observation, we first examine the singular value distribution of \(A_{\mathcal{I}_r}P\) under five selection strategies for constructing \(\mathcal{I}_p\), namely CPQR, SVD, SqNorm, SkCPQR and RBRP.

For the synthetic matrices \(A=UDV^\top\), \(U\) is obtained from the thin QR factorization \(G=UR\) of a Gaussian random matrix \(G\in\mathbb{R}^{m\times r}\). For any index set \(\mathcal{I}_p\) with \(|\mathcal{I}_p|=m_p\le r\), the submatrix \(G_{\mathcal{I}_p}\in\mathbb{R}^{m_p\times r}\) has full row rank with probability \(1\) [74]. Since \(R\) is nonsingular with probability \(1\), it follows from \(G_{\mathcal{I}_p}=U_{\mathcal{I}_p}R\) that \(\operatorname{rank}(U_{\mathcal{I}_p}) = \operatorname{rank}(G_{\mathcal{I}_p}) = m_p .\) Moreover, since \(D\) is nonsingular and \(V^\top\) has full row rank, \(DV^\top\) has full row rank. Therefore, from \(A_{\mathcal{I}_p}=U_{\mathcal{I}_p}DV^\top\), we obtain \(\operatorname{rank}(A_{\mathcal{I}_p})=\operatorname{rank}(U_{\mathcal{I}_p}) =m_p\) with probability \(1\). In our experiments, all selected matrices \(A_{\mathcal{I}_p}\) are observed to have full row rank, which implies that \(r_p:=\operatorname{rank}(A_{\mathcal{I}_p})=m_p\). By Lemma 1, the nonzero singular values of \(A_{\mathcal{I}_r}P\) satisfy \[\sigma_{i+m_p}(A) = \sigma_{i+r_p}(A) \leq \sigma_i(A_{\mathcal{I}_r}P) \leq \sigma_i(A),\;i=1,\ldots,r-m_p.\] Figure 3 illustrates these interlacing bounds for different values of \(m_p\). The nonzero singular values of \(A_{\mathcal{I}_r}P\) lie within the shaded region bounded above by \(\sigma_i(A)\) and below by \(\sigma_{i+m_p}(A)\). Compared with the original matrix \(A\), the reduced matrix \(A_{\mathcal{I}_r}P\) has a more concentrated singular-value distribution.

We next compare the five SC-IS-Krylov-PS variants obtained from the selection strategies introduced in Section 3.2, namely SVD-IS-Krylov-PS, SqNorm-IS-Krylov-PS, CPQR-IS-Krylov-PS, SkCPQR-IS-Krylov-PS, and RBRP-IS-Krylov-PS. Figure 4 shows the RSE with respect to the number of iterations and CPU time under \(m_p=150,300\) and \(450\). It can be seen from Figure 4 that the five strategies require comparable numbers of iterations. However, SVD-IS-Krylov-PS, CPQR-IS-Krylov-PS and RBRP-IS-Krylov-PS incur higher CPU time compared to the other methods, mainly due to the auxiliary matrix constructions, factorizations, and adaptive updates required in their row-selection procedures. SkCPQR-IS-Krylov-PS reduces the CPU time compared with CPQR-IS-Krylov-PS by performing pivot selection in a lower-dimensional sketched space. Indeed, SqNorm-IS-Krylov-PS is still faster across our experiments, since it only requires squared-row-norm sampling and avoids QR-type factorizations. Therefore, we focus our subsequent tests on the SqNorm-IS-Krylov-PS method.

5.3 Choice of parameters \(q\), \(\ell\), and \(m_p\)↩︎

In this subsection, we investigate the influence of the block size \(q\), the truncation parameter \(\ell\), and the number of selected constraint rows \(m_p=|\mathcal{I}_p|\) on the convergence behavior of SqNorm-IS-Krylov-PS. We note that the case \(\ell=1\) corresponds to Algorithm 1 with the SqNorm selection strategy and partition sampling, denoted by SqNorm-RIM-PS. The performance is measured by both the CPU time and the number of full iterations \(k\cdot \frac{q}{m_r}\), where \(m_r=m-m_p\).

It can be seen from the first row of Figure 5 that, for fixed values of \(m_p\) and \(q\), increasing the truncation parameter \(\ell\) reduces the number of full iterations. In particular, the case \(\ell=1\), which corresponds to SqNorm-RIM-PS, requires more full iterations than the cases \(\ell>1\). Moreover, for fixed \(q\) and \(\ell\), increasing \(m_p\) generally leads to fewer full iterations, indicating that using more selected constraint rows can improve the convergence behavior. In terms of CPU time, for a fixed \(m_p\), the SqNorm-RIM-PS method first becomes faster and then slower as the block size \(q\) increases. For the SqNorm-IS-Krylov-PS method with \(\ell=2,5,10,\) and \(30\), increasing \(\ell\) generally reduces the CPU time, which shows the benefit of using more orthogonalized historical search directions. Moreover, for a fixed \(\ell\), increasing \(m_p\) also decreases the CPU time, which is consistent with the reduction in the number of full iterations. These observations indicate that SqNorm-IS-Krylov-PS can substantially improve the efficiency of the basic SqNorm-RIM-PS method. When both CPU time and the number of full iterations are taken into account, configurations with \(q=32\) or \(64\), \(\ell=10\) or \(30\), and \(m_p=16\) or \(64\) provide favorable performance, achieving a desirable balance between convergence speed, computational cost, and storage cost.

5.4 Comparison with IS-Krylov-PS↩︎

We next compare SqNorm-IS-Krylov-PS with its non-preconditioned counterpart, namely the IS-Krylov-PS method. Table [tab:1] and Figure 6 report the number of iterations and CPU time of the two methods on real-world matrices from the SuiteSparse Matrix Collection kolodziej2019suitesparse? and LIBSVM chang2011libsvm?. The SuiteSparse test set consists of bibd_16_8, crew1, WorldCities, lp_ship04s, cr42, D_6, D_7, rel6, mk12-b2, and abtaha2. The LIBSVM test set consists of a9a, aloi, and protein. These matrices cover full-rank and rank-deficient cases, as well as overdetermined and underdetermined systems.

Table [tab:1] shows that SqNorm-IS-Krylov-PS requires fewer iterations than IS-Krylov-PS for all SuiteSparse test matrices. This indicates that the subspace-constrained preconditioning process effectively accelerates the convergence of IS-Krylov-PS. In most cases, the reduced iteration count also leads to a smaller CPU time. The only exception is rel6, where the additional preconditioning and projection costs offset the savings from fewer iterations. Figure 6 reports the corresponding results on the LIBSVM datasets. For all three datasets, SqNorm-IS-Krylov-PS requires fewer iterations and less CPU time than the IS-Krylov-PS method.

Average number of iterations and CPU time of IS-Krylov-PS and SqNorm-IS-Krylov-PS for matrices from the SuiteSparse Matrix Collection kolodziej2019suitesparse?. All computations are terminated once RSE\(<10^{-12}\).

5.5 Comparison with pinv and lsqminnorm↩︎

We compare IS-Krylov-PS and SqNorm-IS-Krylov-PS with the MATLAB built-in functions pinv and lsqminnorm. We generate full-column-rank matrices \(A\in\mathbb{R}^{m\times n}\) and vary the number of rows \(m\), while keeping the number of columns fixed. The exact solution is generated as \(x^*=\texttt{randn(n,1)}\), and the right-hand side is set to \(b=Ax^*\). The iterative methods are terminated when their solution accuracy reaches the accuracy level of the MATLAB built-in solvers.

Figure 7 compares the CPU time of pinv, lsqminnorm, IS-Krylov-PS, and SqNorm-IS-Krylov-PS as the number of rows \(m\) increases. The results show that the CPU time of pinv increases rapidly with \(m\), especially in the larger-scale settings. The function lsqminnorm is competitive for small values of \(m\), but its CPU time increases more noticeably as \(m\) grows. For larger \(m\), SqNorm-IS-Krylov-PS generally requires less CPU time than lsqminnorm. Moreover, SqNorm-IS-Krylov-PS consistently improves upon IS-Krylov-PS in all three settings, which further demonstrates the benefit of the subspace-constrained preconditioning.

6 Concluding remarks↩︎

This paper further investigated the subspace-constrained preconditioning technique based on a QR-like factorization and verified its applicability to general linear systems, including underdetermined, overdetermined, full-rank, and rank-deficient matrices. A projected randomized iterative method was subsequently adopted to solve the preconditioned linear systems, and it was confirmed that the computational cost could be effectively reduced when the initial point is chosen as a solution of the subsystem \(A_{\mathcal{I}_p} x = b_{\mathcal{I}_p}\). Furthermore, drawing on the principles of minimal error methods and the orthogonal direction method, the proposed subspace-constrained iterative framework was further accelerated, giving rise to a novel subspace-constrained iterative-sketching-based Krylov subspace (SC-IS-Krylov) method. Theoretical analysis showed that the proposed method achieves linear convergence in expectation. Numerical experiments validated our theoretical results and demonstrated the superior computational efficiency of the proposed method.

There are still many possible future avenues of research. Since linear systems arising in practical problems are often inconsistent due to noise [75]–[77], extending the subspace-constrained preconditioning technique to such cases is a valuable direction. The randomized sparse Kaczmarz method proposed in [78], [79] has been recognized as an effective strategy for obtaining sparse solutions to linear systems. A potential avenue for future research could be the extension of the subspace-constrained preconditioning technique to address sparse recovery problems.

7 Appendix. Proof of the main results↩︎

7.1 Proof of Lemma 1↩︎

We first introduce the Poincaré separation theorem, which provides an interlacing inequality for eigenvalues.

7.2 Proof of Theorem 1↩︎

To facilitate the convergence analysis, we introduce several auxiliary variables and notations associated with the \(k\)-th iteration. At the \(k\)-th iteration, we consider the product probability space \(\left( \prod_{i=0}^k \Omega, \otimes_{i=0}^k \mathcal{F}, \tilde{\mathbf{P}} \right),\) where \(\otimes\) denotes the product of \(\sigma\)-algebras and \(\tilde{\mathbf{P}}\) is the corresponding product measure, as defined in [70]. Let \(\mathcal{B}_k:= \sigma(S_0, S_1, \ldots, S_{k-1})\) be the \(\sigma\)-algebra generated by the random variables \((S_0, S_1, \ldots, S_{k-1})\), where \(\mathcal{B}_0\) is the trivial \(\sigma\)-algebra (i.e., \(\mathcal{B}_0 = \{\emptyset, \Omega\}\)). We denote the conditional expectation with respect to \(\mathcal{B}_k\) as \[\mathbb{E}_{k}[\cdot] := \mathbb{E}[\cdot \mid \mathcal{B}_k].\] Note that if \(X\) and \(XY\) are integrable random variables, and \(X\) is measurable with respect to \(\mathcal{B}_k\), then the following identities hold [70]: \[\label{prop-Econd} \begin{align} \mathbb{E}[X | \mathcal{B}_k] = X \quad \text{and} \quad \mathbb{E}[XY| \mathcal{B}_k] = X \mathbb{E}[Y| \mathcal{B}_k]. \end{align}\tag{28}\] Let \(\bar S_k\) be a sampling matrix drawn from the probability space \((\Omega,\mathcal{F},\mathbf{P})\). We define the auxiliary iterate \[\bar x^{k+1}:=x^k-\bar\alpha_k P A_{\mathcal{I}_r}^{\top} \bar S_k\bar S_k^{\top} (A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r}).\] The step-size \(\bar\alpha_k\) is chosen as \[\label{eq-bar-alpha} \bar\alpha_k = \begin{cases} (2-\zeta)\frac{ \|\bar S_k^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\|_2^2 }{ \|P A_{\mathcal{I}_r}^{\top} \bar S_k\bar S_k^\top (A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\|_2^2 }, & \text{if }\bar S_k\in\mathcal{Q}_k,\\[1mm] 0, & \text{if }\bar S_k\in\mathcal{Q}_k^c, \end{cases}\tag{29}\] where \(\mathcal{Q}_k\) is defined in 24 . For any measurable set \(\mathcal{Q}\subseteq\Omega\), we use \(\mathbf{P}(\bar S_k\in\mathcal{Q}\mid\mathcal{B}_k)\) to denote the conditional probability of the event \(\{\bar S_k\in\mathcal{Q}\}\) given \(\mathcal{B}_k\). The complement of \(\mathcal{Q}\) is defined as \(\mathcal{Q}^c = \Omega \setminus \mathcal{Q}\). The indicator function \(\mathbb{I}_{\mathcal{Q}}(\cdot)\) of the set \(\mathcal{Q}\) is defined as \[\mathbb{I}_{\mathcal{Q}}(\bar{S}_k) = \begin{cases} 1 & \text{if }\bar S_k\in\mathcal{Q};\\ 0 & \text{otherwise}. \end{cases}\]

7.3 Proof of Proposition 5 and Lemma 7↩︎

To prove the Proposition 5, we first introduce the change-of variables theorem.

7.4 Proof of Lemma 8↩︎

We first introduce a useful lemma.

7.5 Proof of Theorem 2↩︎

Let \(\hat{S}_k\) be a sampling matrix drawn from the probability space \((\Omega,\mathcal{F},\mathbf{P})\). We define the auxiliary quantities \[\label{eq-hat-xk431} \left\{ \begin{align} \hat{d}^k &=-P A_{\mathcal{I}_r}^{\top} \hat{S}_k \hat{S}_k^\top (A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r}),\\ \hat{p}^k &=\hat{d}^k-\sum_{i=j_{k,\ell}}^{k-1} \frac{\langle \hat{d}^k,p^i\rangle}{\lVert p^i\rVert_2^2}p^i,\\ \hat{x}^{k+1} &=x^k+\hat{\delta}_k\hat{p}^k, \end{align} \right.\tag{41}\] where \[\label{eq-hat-delta} \hat{\delta}_k = \begin{cases} \frac{ \|\hat{S}_k^\top(A_{\mathcal{I}_r}x^k-b_{\mathcal{I}_r})\|_2^2 }{ \|\hat{p}^k\|_2^2 }, & \text{if }\hat{S}_k\in\mathcal{Q}_k,\\[1mm] 0, & \text{if }\hat{S}_k\in\mathcal{Q}_k^c, \end{cases}\tag{42}\] and \(\mathcal{Q}_k\) is defined in 24 .

7.6 Proof of Lemma 9 and Theorem 3↩︎

References↩︎