## 1 Introduction
As just mentioned, there are many algorithms for solving linear least-squares problems. For small- and medium-sized problems, one can use QR-based methods, and, if the problem is rank-deficient or ill conditioned, singular value decomposition (SVD) can be adopted [37; 38]. For large-sized least-squares problems, Krylov subspace methods are more efficient than direct methods if the matrix is sparse or otherwise structured. In the present paper, we plan to develop an efficient and accurate Krylov subspace method for large-sized least-squares problems with dense and maybe ill-conditioned matrices, not limited to sparse matrices. It is also interesting that the proposed iterative algorithm can solve least-squares problems with extremal cases of \(q\ll n\) or \(q\gg n\). In fact, for under-determined systems, many methods are not applicable, for instance, QR and conjugate gradient (CG)-type methods, for which the basis vectors of the Krylov subspace are orthogonalized implicitly. In addition to Equation (1), [PERSON] [39] considered the following minimization  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]1  
for solving ill-posed linear systems of Equation (2) with \(q=n\). [PERSON] [40; 41] employed the double optimization technique to iteratively solve inverse problems. [PERSON] [42] replaced the Krylov subspace by the column space of \(\mathbf{A}\) and derived an effective iterative algorithm based on the maximal projection in Equation (12). These results are quite promising. However, for a non-square system (2), there exists no similar method of double optimization in the literature. It is well known that the least-squares solution of the minimization in Equation (1) is obtained by  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]2  
When the derivative of \(f\) in Equation (12) vs. \(\mathbf{x}\) is equal to zero, we can obtain  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]3  
Taking the transpose of Equation (2) to \(\mathbf{x}^{\mathsf{T}}\mathbf{A}^{\mathsf{T}}=\mathbf{b}^{\mathsf{T}}\) and inserting it into Equation (14), we can immediately obtain  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]4  
which implies Equation (13) by \(\mathbf{b}^{\mathsf{T}}\mathbf{A}\mathbf{x}/\mathbf{b}^{\mathsf{T}}\mathbf{A} \mathbf{x}=1\). Many numerical methods for least-squares problems are based on Equation (1); no numerical methods for least-squares problems are based on Equation (12). To solve the least-squares problem (2), the usual approach is the Karush-Kuhn-Tucker equation method, which sets Equation (2) to a square linear system but with a larger dimension \(q+n\). Let  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]5  
be the residual vector. By considering the enforcing constraint:  
\[\mathbf{Ax}=\mathbf{b}, \tag{2}\]6  
which is equivalent to the normal form of Equation (2), Equations (16) and (17) constitute a \(q+n\)-dimensional linear system:  
\[\left[ \begin{array}{cc}\mathbf{I}_{q}&\mathbf{A}\\ \mathbf{A}^{\mathsf{T}}&\mathbf{0}\end{array} \right]\left[ \begin{array}{c} \mathbf{r}\\ \mathbf{x}\end{array} \right]=\left[ \begin{array}{c}\mathbf{b}\\ \mathbf{0}\end{array} \right]. \tag{18}\]  
The advantage of the so-called Karush-Kuhn-Tucker equation [15] is that one can employ a linear solver to solve Equation (18) to obtain the least-squares solution of Equation (2); however, a disadvantage is that the dimension is enlarged to \(q+n\). A major novelty of the present paper is that we employ both Equations (1) and (12) to develop a novel numerical method in a varying affine Krylov subspace. The innovative contributions are as follows: this is the first time a new concept of varying affine Krylov subspaces has been introduced, and a double-optimal iterative algorithm has been derived based on the double optimization method to solve least-squares problems. A new variant of the Karush-Kuhn-Tucker equation can significantly improve the accuracy of solving least-squares problems by adopting the partial pivoting Gaussian elimination method.