## 1 Introduction  
One of the most important numerical methods to solve linear systems is the Krylov subspace method [1, 2, 3, 4, 5, 6, 7]; since the appearance of pioneering works in [8, 9], it has been studied extensively. A lot of Krylov subspace methods were reported in [10, 11, 12, 13, 14, 15, 16]. Many algorithms have been used to solve the least-squares problem with a minimum norm [17, 18]:  
\[\min\|\mathbf{x}\|\ \ \text{s.t.}\ \ \mathbf{x}\in\arg\min_{\mathbf{x}}\| \mathbf{b}-\mathbf{Ax}\|, \tag{1}\]  
where \(\|\cdot\|\) denotes the Euclidean norm, \(\mathbf{x}\in\mathbb{R}^{n}\) is unknown, and \(\mathbf{A}\in\mathbb{R}^{q\times n}\) and \(\mathbf{b}\in\mathbb{R}^{q}\) are given. We are concerned with a varying affine Krylov subspace (VAKS) solution and its corresponding iterative algorithm for solving  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]  
which is an over-determined system when \(q>n\) and an under-determined system when \(q<n\). Most results for Equation (2) are directly extendable to Equation (1), whose solution is called the minimum-norm or pseudo-inverse solution and is given by \(\mathbf{x}=\mathbf{A}^{\dagger}\mathbf{b}\), where \(\mathbf{A}^{\dagger}\) denotes the pseudo-inverse of \(\mathbf{A}\) satisfying the following Penrose equations:  
\[\mathbf{AA}^{\dagger}\mathbf{A}=\mathbf{A}, \tag{3}\] \[\mathbf{A}^{\dagger}\mathbf{AA}^{\dagger}=\mathbf{A}^{\dagger},\] (4) \[(\mathbf{AA}^{\dagger})^{\mathrm{T}}=\mathbf{AA}^{\dagger},\] (5) \[(\mathbf{A}^{\dagger}\mathbf{A})^{\mathrm{T}}=\mathbf{A}^{\dagger }\mathbf{A}. \tag{6}\]  
The pseudo-inverse has been investigated and applied to solve linear least-squares problems, and some computational methods have been developed to compute the Moore-Penrose inverse matrix [19; 20; 21; 22; 23; 24]. In terms of the Moore-Penrose inverse matrix \(\mathbf{A}^{\dagger}\), the general solution of the least-squares problem is given by  
\[\mathbf{x}=\mathbf{A}^{\dagger}\mathbf{b}+(\mathbf{I}-\mathbf{R} _{A})\mathbf{z}, \tag{7}\] \[\mathbf{R}_{A}=\mathbf{A}^{\dagger}\mathbf{A}, \tag{8}\]  
where \(\mathbf{z}\) is arbitrary. Because \(\mathbf{A}^{\dagger}\mathbf{b}\) and \((\mathbf{I}-\mathbf{R}_{A})\mathbf{z}\) are orthogonal, it follows that \(\mathbf{x}=\mathbf{A}^{\dagger}\mathbf{b}\) is the unique solution of Equation (1). QR factorization is a good method to solve Equation (2). The normalized Gram-Schmidt process is a modification of the original Gram-Schmidt orthogonalization by normalizing each vector after orthogonalization. Let \(\mathbf{Q}\) denote a \(q\times n\) matrix obtained by applying the normalized Gram-Schmidt process on the column vectors of \(\mathbf{A}\):  
\[\mathbf{Q}:=[\mathbf{q}_{1},\dots,\mathbf{q}_{n}], \tag{9}\]  
which possesses the following orthonormal property:  
\[\mathbf{Q}^{\mathrm{T}}\mathbf{Q}=\mathbf{I}_{n}. \tag{10}\]  
Hence, by \(\mathbf{QRx}=\mathbf{b}\), it is easy to solve \(\mathbf{x}\) by applying the backward substitution method to  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]0  
However, the QR method is only applicable to an over-determined system with \(q>n\). For an under-determined system with \(q<n\), the QR method is not applicable because \(\mathbf{R}\) is no longer an upper triangular matrix. Recently, many scholars have proposed many algorithms to tackle least-squares problems, like relaxed greedy randomized coordinate descent methods [25], the parallel approximate pseudo-inverse matrix technique in conjunction with the explicit preconditioned conjugate gradient least squares scheme [26], the new randomized Gauss-Seidel method [27], the QR-Cholesky method, the implicit QR method, the bidiagonal implicit QR method [28], the randomized block coordinate descent algorithm [29], an iterative pre-conditioning technique [30], splitting-based randomized iterative methods [31], the projection method [32], the index-search-method-based inner-outer iterative algorithm [33], the greedy double subspace coordinate descent method [34], a distributed algorithm [35], and multivalued projections [36].