## 4 A Numerical Algorithm  
In this section, we develop an iterative algorithm to solve the least-squares problem, starting from an initial guess \(\mathbf{x}_{0}\). We assume that the value of \(\mathbf{x}_{k}\) at the \(k\)th step is already known, and \(\mathbf{x}_{k+1}\) is computed at the next step via the iterative algorithm. According to the value \(\mathbf{x}_{k}\), the \(k\)th step residual \(\mathbf{r}_{k}=\mathbf{b}-\mathbf{A}\mathbf{x}_{k}\) can be computed. When the initial guess \(\mathbf{x}_{0}\) is given, an initial residual is written as follows:  
\[\mathbf{r}_{0}=\mathbf{b}-\mathbf{A}\mathbf{x}_{0}. \tag{109}\]  
Upon letting  
\[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]  
Equation (2) is equivalent to solving \(\mathbf{z}\) from  
\[\mathbf{A}\mathbf{z}=\mathbf{r}_{0}, \tag{111}\]  
which is deduced by  
\[\mathbf{A}\mathbf{z}=\mathbf{A}\mathbf{x}-\mathbf{A}\mathbf{x}_{0}=\mathbf{b}- \mathbf{A}\mathbf{x}_{0}=\mathbf{r}_{0}.\]  
For system (111), we seek \(\mathbf{z}\) in the VAKS by  
\[\mathbf{z}=\alpha_{0}\mathbf{u}_{0}+\sum_{k=1}^{m}\alpha_{k}\mathbf{u}_{k}= \alpha_{0}\mathbf{u}_{0}+\mathbf{U}\mathbf{\alpha}\in\mathcal{K}^{\prime}_{m}, \tag{112}\]  
where  
\[\mathbf{u}_{0}=\mathbf{A}^{\mathrm{T}}\mathbf{r}_{0}, \tag{113}\]  
and \(\mathbf{u}_{k}:=(\mathbf{A}^{\mathrm{T}}\mathbf{A})^{k}\mathbf{u}_{0},\ k=1, \ldots,m\). The constants \(\alpha_{k},\ k=1,\ldots,m\) and \(\alpha_{0}\) are determined by the following two minimizations:  
\[\min_{\mathbf{z}}\biggl{\{}f=\frac{\|\mathbf{A}\mathbf{z}\|^{2} }{(\mathbf{r}_{0}^{\mathrm{T}}\mathbf{A}\mathbf{z})^{2}}\biggr{\}}, \tag{114}\] \[\min_{\mathbf{z}}r=\|\mathbf{r}_{0}-\mathbf{A}\mathbf{z}\|^{2}. \tag{115}\]  
After inserting Equation (112) for \(\mathbf{z}\), the first minimization can derive \(\mathbf{\alpha}\), while the second minimization can derive \(\alpha_{0}\). Since two minimizations in Equations (114) and (115) are adopted to determine the descent vector \(\mathbf{z}\) in Equation (112), the resulting iterative algorithm to solve the least-square problem might be labeled as a double optimal algorithm (DOA). To treat the rank deficient least-squares problem, the dimension \(m\) of the VAKS must be \(m<m_{0}=\mathrm{rank}(\mathbf{A})\), such that \(\mathrm{rank}(\mathbf{J}_{k})=m\) and \(\mathbf{C}_{k}=\mathbf{J}_{k}^{\mathrm{T}}\mathbf{J}_{k}>0\) is an \(m\times m\) invertible matrix. Consequently, the DOA is depicted by (i), giving \(\mathbf{A}\), \(\mathbf{b}\), \(m<m_{0}\), an initial value \(\mathbf{x}_{0}\) and the convergence criterion \(\varepsilon\) and (ii) doing \(k=0,1,\ldots\),\[\mathbf{r}_{k}=\mathbf{b}-\mathbf{A}\mathbf{x}_{k},\] \[\mathbf{u}_{0}^{k}=\mathbf{A}^{\mathrm{T}}\mathbf{r}_{k},\] \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]0 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]1 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]2 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]3 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]4 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]5 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]6 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]7 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]8 \[\mathbf{z}=\mathbf{x}-\mathbf{x}_{0}, \tag{110}\]9  
until the convergence with \(\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|<\varepsilon\) or \(\|\mathbf{r}_{k+1}\|<\varepsilon\). We call \(\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|\) the relative residual. In the above, \(\mathbf{U}_{k}\) is \(n\times m\) matrix, \(\mathbf{J}_{k}\) is a \(q\times m\) matrix, \(\mathbf{C}_{k}\) and \(\mathbf{D}_{k}\) are \(m\times m\) matrices, \(\mathbf{V}_{k}\) is an \(n\times q\) matrix and \(\mathbf{P}_{k}\) is a \(q\times q\) matrix. The computational cost is expanded to compute \(\mathbf{D}_{k}\), which is however quite time-saving because \(m\) is a small number. \(\mathbf{A}^{\mathrm{T}}\mathbf{A}\) is an \(n\times n\) fixed matrix computed once and used in the construction of \(\mathbf{U}_{k}\), which requires \(mn\) operations.