文章目录

• • 最小二乘法
  • 返回值
  • 测试

最小二乘法

scipy.sparse.linalg实现了两种稀疏矩阵最小二乘法lsqr和lsmr，前者是经典算法，后者来自斯坦福优化实验室，据称可以比lsqr更快收敛。

这两个函数可以求解Ax=bAx=bAx=b，或arg min⁡x∥Ax−b∥2\argmin_x\Vert Ax-b\Vert^2argminx​∥Ax−b∥2，或arg min⁡x∥Ax−b∥2+d2∥x−x0∥2\argmin_x\Vert Ax-b\Vert^2+d^2\Vert x-x_0\Vert^2argminx​∥Ax−b∥2+d2∥x−x0​∥2，其中AAA必须是方阵或三角阵，可以有任意秩。

通过设置容忍度at,bta_t, b_tat​,bt​，可以控制算法精度，记r=b−Axr=b-Axr=b−Ax为残差向量，如果Ax=bAx=bAx=b是相容的，lsqr在∥r∥⩽at∗∥A∥⋅∥x∥+bt∥b∥\Vert r\Vert\leqslant a_t*\Vert A\Vert\cdot\Vert x\Vert + b_t\Vert b\Vert∥r∥⩽at​∗∥A∥⋅∥x∥+bt​∥b∥时终止；否则将在∥ATr∥⩽at∥A∥⋅∥r∥\Vert A^T r\Vert\leqslant a_t\Vert A\Vert \cdot\Vert r\Vert∥ATr∥⩽at​∥A∥⋅∥r∥。
如果两个容忍度都是10−610^{-6}10−6，最终的∥r∥\Vert r\Vert∥r∥将有6位精度。

lsmr的参数如下

lsmr(A, b, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, maxiter=None, show=False, x0=None)

参数解释：

• A 可谓稀疏矩阵、数组以及线性算子
• b 为数组
• damp 阻尼系数，默认为0
• atol, btol 截止容忍度，是lsqr迭代的停止条件，即at,bta_t, b_tat​,bt​。
• conlim 另一个截止条件，对于最小二乘问题，conlim应该小于10810^8108，如果Ax=bAx=bAx=b是相容的，则conlim最大可以设到101210^{12}1012
• iter_limint 迭代次数
• show 如果为True，则打印运算过程
• calc_var 是否估计(A.T@A + damp**2*I)^{-1}的对角线
• x0 阻尼系数相关

lsqr和lsmr相比，没有maxiter参数，但多了iter_lim, calc_va参数。

上述参数中，damp为阻尼系数，当其不为0时，记作δ\deltaδ，待解决的最小二乘问题变为

[AδI]x=[bδx0]\begin{bmatrix}A\\\delta I\end{bmatrix} x=\begin{bmatrix}b\\\delta x_0 \end{bmatrix} [AδI​]x=[bδx0​​]

返回值

lsmr的返回值依次为：

• x 即Ax=bAx=bAx=b中的xxx
• istop 程序结束运行的原因
• itn 迭代次数
• normr ∥b−Ax∥\Vert b-Ax\Vert∥b−Ax∥
• normar ∥AT(b−Ax)∥\Vert A^T(b-Ax)\Vert∥AT(b−Ax)∥
• norma ∥A∥\Vert A\Vert∥A∥
• conda A的条件数
• normx ∥x∥\Vert x\Vert∥x∥

lsqr的返回值为

• x 即Ax=bAx=bAx=b中的xxx
• istop 程序结束运行的原因
• itn 迭代次数
• r1norm ∥b−Ax∥\Vert b-Ax\Vert∥b−Ax∥
• r2norm ∥b−Ax∥2+δ2∥x−x0∥2\sqrt{\Vert b-Ax\Vert^2+\delta^2\Vert x-x_0\Vert^2}∥b−Ax∥2+δ2∥x−x0​∥2​
• anorm 估计的Frobenius范数Aˉ\bar AAˉ
• acond Aˉ\bar AAˉ的条件数
• arnorm ∥ATr−δ2(x−x0)∥\Vert A^Tr-\delta^2(x-x_0)\Vert∥ATr−δ2(x−x0​)∥
• xnorm ∥x∥\Vert x\Vert∥x∥
• var (ATA)−1(A^TA)^{-1}(ATA)−1

二者的返回值较多，而且除了前四个之外，剩下的意义不同，调用时且须注意。

测试

下面对这两种算法进行验证，第一步就得先有一个稀疏矩阵

import numpy as np
from scipy.sparse import csr_arraynp.random.seed(42)  # 设置随机数状态
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)

然后用这个稀疏矩阵乘以一个xxx，得到bbb

xs = np.arange(500)
b = mat @ xs

接下来对这两个最小二乘函数进行测试

from scipy.sparse.linalg import lsmr, lsqr
import matplotlib.pyplot as plt
mx = lsmr(csr, b)[0]
qx = lsqr(csr, b)[0]
plt.plot(xs, lw=0.5)
plt.plot(mx, lw=0, marker='*', label="lsmr")
plt.plot(qx, lw=0, marker='.', label="lsqr")
plt.legend()
plt.show()

为了对比清晰，对图像进行放大，可以说二者不分胜负

接下来比较二者的效率，500×500500\times500500×500这个尺寸显然已经不合适了，用2000×20002000\times20002000×2000

from timeit import timeitnp.random.seed(42)  # 设置随机数状态
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)
timeit(lambda : lsmr(csr, b), number=10)
timeit(lambda : lsqr(csr, b), number=10)

测试结果如下

>>> timeit(lambda : lsqr(csr, b), number=10)
0.5240591000001587
>>> timeit(lambda : lsmr(csr, b), number=10)
0.6156221000019286

看来lsmr并没有更快，看来斯坦福也不靠谱(滑稽)。