最小二乘法

scipy.sparse.linalg實(shí)現(xiàn)了兩種稀疏矩陣最小二乘法lsqr和lsmr，前者是經(jīng)典算法，后者來自斯坦福優(yōu)化實(shí)驗(yàn)室，據(jù)稱可以比lsqr更快收斂。

這兩個(gè)函數(shù)可以求解Ax=b，或arg?min_x?∥Ax?b∥²，或arg?min_x?∥Ax?b∥²?+d²∥x?x0∥²，其中A必須是方陣或三角陣，可以有任意秩。

通過設(shè)置容忍度a_t ,b_t,可以控制算法精度，記r=b-A_x?為殘差向量，如果A_x=b是相容的，lsqr在∥r∥?a_t?∥A∥?∥x∥+b_t∥b∥時(shí)終止；否則將在∥A^T_r∥?a_t∥A∥?∥r∥。

如果兩個(gè)容忍度都是10^?6?，最終的∥r∥將有6位精度。

lsmr的參數(shù)如下

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

參數(shù)解釋：

• A 可謂稀疏矩陣、數(shù)組以及線性算子
• b 為數(shù)組
• damp 阻尼系數(shù)，默認(rèn)為0
• atol, btol 截止容忍度，是lsqr迭代的停止條件，即a_t ,b_t?。
• conlim 另一個(gè)截止條件，對(duì)于最小二乘問題，conlim應(yīng)該小于10⁸，如果A_x=b是相容的，則conlim最大可以設(shè)到10¹²
• iter_limint 迭代次數(shù)
• show 如果為True，則打印運(yùn)算過程
• calc_var 是否估計(jì)(A.T@A + damp**2*I)^{-1}的對(duì)角線
• x0 阻尼系數(shù)相關(guān)

lsqr和lsmr相比，沒有maxiter參數(shù)，但多了iter_lim, calc_va參數(shù)。

上述參數(shù)中，damp為阻尼系數(shù)，當(dāng)其不為0時(shí)，記作δ，待解決的最小二乘問題變?yōu)?/p>

返回值

lsmr的返回值依次為：

• x 即A_x=b中的x
• istop 程序結(jié)束運(yùn)行的原因
• itn 迭代次數(shù)
• normr ∥b?A_x∥
• normar ∥A^T?(b?Ax)∥
• norma ∥A∥
• conda A的條件數(shù)
• normx ∥x∥

lsqr的返回值為

1. x 即A_x=b中的x
2. istop 程序結(jié)束運(yùn)行的原因
3. itn 迭代次數(shù)
4. r1norm
5. anorm 估計(jì)的Frobenius范數(shù)Aˉ
6. acond Aˉ的條件數(shù)
7. arnorm ∥A^T_r?δ²(x?x₀)∥
8. xnorm ∥x∥
9. var (A^TA)^?1

二者的返回值較多，而且除了前四個(gè)之外，剩下的意義不同，調(diào)用時(shí)且須注意。

測(cè)試

下面對(duì)這兩種算法進(jìn)行驗(yàn)證，第一步就得先有一個(gè)稀疏矩陣

import numpy as np
from scipy.sparse import csr_array

np.random.seed(42)  # 設(shè)置隨機(jī)數(shù)狀態(tài)
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)

然后用這個(gè)稀疏矩陣乘以一個(gè)x，得到b

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

接下來對(duì)這兩個(gè)最小二乘函數(shù)進(jìn)行測(cè)試

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()

為了對(duì)比清晰，對(duì)圖像進(jìn)行放大，可以說二者不分勝負(fù)

接下來比較二者的效率，500 × 500 500\times500500×500這個(gè)尺寸顯然已經(jīng)不合適了，用2000×2000

from timeit import timeit

np.random.seed(42)  # 設(shè)置隨機(jī)數(shù)狀態(tài)
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)

測(cè)試結(jié)果如下

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

看來lsmr并沒有更快，看來斯坦福也不靠譜(滑稽)。