Python實現兩種稀疏矩陣的最小二乘法

Posted on by

最小二乘法

scipy.sparse.linalg實現瞭兩種稀疏矩陣最小二乘法lsqr和lsmr，前者是經典算法，後者來自斯坦福優化實驗室，據稱可以比lsqr更快收斂。

這兩個函數可以求解Ax=b，或arg minx ∥Ax−b∥2，或arg minx ∥Ax−b∥2 +d2∥x−x0∥2，其中A必須是方陣或三角陣，可以有任意秩。

通過設置容忍度at ,bt,可以控制算法精度，記r=b-Ax 為殘差向量，如果Ax=b是相容的，lsqr在∥r∥⩽at∗∥A∥⋅∥x∥+bt∥b∥時終止；否則將在∥ATr∥⩽at∥A∥⋅∥r∥。

如果兩個容忍度都是10−6 ，最終的∥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 ,bt 。
  • conlim 另一個截止條件，對於最小二乘問題，conlim應該小於108，如果Ax=b是相容的，則conlim最大可以設到1012
  • iter_limint 迭代次數
  • show 如果為True，則打印運算過程
  • calc_var 是否估計(A.T@A + damp**2*I)^{-1}的對角線
  • x0 阻尼系數相關

lsqr和lsmr相比，沒有maxiter參數，但多瞭iter_lim, calc_va參數。

上述參數中，damp為阻尼系數，當其不為0時，記作δ，待解決的最小二乘問題變為

返回值

lsmr的返回值依次為：

  • x 即Ax=b中的x
  • istop 程序結束運行的原因
  • itn 迭代次數
  • normr ∥b−Ax∥
  • normar ∥AT (b−Ax)∥
  • norma ∥A∥
  • conda A的條件數
  • normx ∥x∥

lsqr的返回值為

  1. x 即Ax=b中的x
  2. istop 程序結束運行的原因
  3. itn 迭代次數
  4. r1norm
  5. anorm 估計的Frobenius范數Aˉ
  6. acond Aˉ的條件數
  7. arnorm ∥ATr−δ2(x−x0)∥
  8. xnorm ∥x∥
  9. var (ATA)−1

二者的返回值較多，而且除瞭前四個之外，剩下的意義不同，調用時且須註意。

測試

下面對這兩種算法進行驗證，第一步就得先有一個稀疏矩陣

import numpy as np
from scipy.sparse import csr_array

np.random.seed(42)  # 設置隨機數狀態
mat = np.random.rand(500,500)
mat[mat<0.9] = 0
csr = csr_array(mat)

然後用這個稀疏矩陣乘以一個x，得到b

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 × 500 500\times500500×500這個尺寸顯然已經不合適瞭，用2000×2000

from timeit import timeit

np.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並沒有更快，看來斯坦福也不靠譜(滑稽)。

以上就是Python實現兩種稀疏矩陣的最小二乘法的詳細內容，更多關於Python稀疏矩陣最小二乘法的資料請關註WalkonNet其它相關文章！

推薦閱讀：