LSQR算法的程序设计（带阻尼与无阻尼）

LSQR算法是计算大型稀疏线性方程组的算法，由Paige和Saunders于1982年提出。该方法主要为求解以下线性方程组（A为m*n的矩阵，m>n），且保持二阶残差范数最小：
$$\begin{gathered} Ax=b \\ Min\parallel Ax-b{\parallel }_2 \end{gathered}$$
对于传统的最小二乘问题，可以直接求出方程的解，但是对于大型稀疏矩阵，往往是病态问题（超定方程组或欠定方程组），因此，LSQR方法不是直接求解x，而是遵守最小二乘的原则，在法方程限定的空间中寻求最佳解。

寻找最佳解的过程如下，可以基于该过程编程：

1. 初始条件
   （a）${\beta }_1{u}_1=b{\alpha }_1{v}_1=A^T{u}_1$
   （b）$w_1=v_1{x}_0=0$
   （c）${\bar{\varphi }}_1={\beta }_1{\rho }_1={\alpha }_1$

2. 对i=1,2,3，…重复步骤（3）-（6）

3. 对角化系数矩阵
   对角化的目的是将系数矩阵双对角化，以便进行QR分解，具体原理我会在另一篇博客说明。分为两个步骤，第一是把 $A{v}_i-{\alpha }_i{u}_i$ 进行标准化，系数是${\beta }_{i+1}$。第二是把$A^T{u}_{i+1}-{\beta }_{i+1}{v}_i$进行标准化，系数是${\alpha }_{i+1}$。
   （a）${\beta }_{i+1}{u}_{i+1}=A{v}_i-{\alpha }_i{u}_i$
   （b）${\alpha }_{i+1}{v}_{i+1}=A^T{u}_{i+1}-{\beta }_{i+1}{v}_i$

4. 计算QR分解的中间变量
   带阻尼的LSQR算法和无阻尼的LSQR算法的区别主要在该步骤中，关于QR分解的原理我会在另一篇博客详细说明。
   无阻尼LSQR算法的QR分解中间变量赋值过程：
   $$\begin{array}{l}(a){\rho }_i=\sqrt{{\bar{\rho }}_i^2+{\beta }_{i+1}^2}\\ (b)c_i={\bar{\rho }}_i/{\rho }_i\\ (c)s_i={\beta }_{i+1}/{\rho }_i\\ (d){\theta }_{i+1}=s_i{\alpha }_{i+1}\\ (e){\bar{\rho }}_{i+1}=-c_i{\alpha }_{i+1}\\ (f){\varphi }_i=c_i{\bar{\varphi }}_i\\ (g){\bar{\varphi }}_{i+1}=s_i{\bar{\varphi }}_i\end{array}$$

带阻尼LSQR算法的QR分解中间变量的赋值过程：
$$\begin{array}{ll}(a)& {\bar{\rho }}_{i+1}=\sqrt{{\bar{\rho }}_i^2+{\lambda }^2}\\ (b)& c_i={\bar{\rho }}_i/{\bar{\rho }}_{i+1}\\ (c)& s_i=\lambda /{\bar{\rho }}_{i+1}\\ (d)& {\psi }_i=s_i{\bar{\varphi }}_i\\ (e)& {\bar{\varphi }}_{i+1}=c_i{\varphi }_i\\ (f)& {\rho }_i=\sqrt{{\bar{\rho }}_i^2+{\lambda }^2+{\beta }_{i+1}^2}\\ (g)& c_{i+1}={\bar{\rho }}_{i+1}/{\rho }_i\\ (h)& s_{i+1}={\beta }_{i+1}/{\rho }_i\\ (i)& {\theta }_{i+1}=s_{i+1}{\alpha }_{i+1}\\ (j)& {\bar{\rho }}_{i+1}^o=-c_{i+1}{\alpha }_{i+1}\\ (k)& {\varphi }_{i+1}=c_{i+1}{\bar{\varphi }}_{i+1}\\ (l)& {\bar{\varphi }}_{i+1}=s_{i+1}{\varphi }_{i+1}\end{array}$$

5. 更新x和w，到此步骤为止，$x_i$和$w_i$能够计算出来了
   （a）$x_{i+1}=x_i+w_i$
   （b）$w_{i+1}=v_{i+1}-w_i$

6. 判断迭代条件，即停机准则
   停机准则有以下三条，涉及到范数的计算，关于范数的计算可以参考我的另一篇博客。
   （a）$\parallel r_k\parallel \le BTOL\parallel b\parallel +ATOL\parallel A\parallel \parallel x\parallel$
   （b）$\frac{\parallel A^T{r}_k\parallel }{\parallel A\parallel \parallel r_k\parallel }\le ATOL$
   （c）$cond(A)>CONLIM$