Python SciPy linalg.lobpcg用法及代碼示例

本文簡要介紹 python 語言中 scipy.sparse.linalg.lobpcg 的用法。

用法:

scipy.sparse.linalg.lobpcg(A, X, B=None, M=None, Y=None, tol=None, maxiter=None, largest=True, verbosityLevel=0, retLambdaHistory=False, retResidualNormsHistory=False, restartControl=20)#

局部最優分塊預條件共軛梯度法 (LOBPCG)。

LOBPCG 是用於大型實對稱和複雜埃爾米特定廣義本征問題的預處理本征求解器。

參數 ：：

A： {稀疏矩陣，ndarray，LinearOperator，可調用對象}

問題的埃爾米特線性算子，通常由稀疏矩陣給出。通常稱為“stiffness matrix”。

X： ndarray、float32 或 float64

初始近似值k特征向量(非稀疏)。如果A已shape=(n,n)然後X一定有shape=(n,k).

B： {稀疏矩陣，ndarray，LinearOperator，可調用對象}

可選的。默認為 B = None ，相當於身份。廣義特征值問題中的右側運算符(如果存在)。通常稱為“mass matrix”。必須是 Hermitian 正定的。

M： {稀疏矩陣，ndarray，LinearOperator，可調用對象}

可選的。默認為 M = None ，相當於身份。預處理器旨在加速收斂。

Y： ndarray、float32 或 float64，默認值：無

一個n-by-sizeYndarray 約束條件sizeY < n。迭代將在B- column-space 的正交補集Y.Y如果存在，則必須是滿級。

tol： 標量，可選

默認為 tol=n*sqrt(eps) 。求解器對停止標準的容差。

maxiter： 整數，默認：20

最大迭代次數。

largest： 布爾值，默認值：真

如果為 True，則求解最大的特征值，否則求解最小的特征值。

verbosityLevel： 整數，可選

默認情況下verbosityLevel=0沒有輸出。控製解算器標準/屏幕輸出。

retLambdaHistory： 布爾值，默認值：假

是否返回迭代特征值曆史。

retResidualNormsHistory： 布爾值，默認值：假

是否返回剩餘規範的迭代曆史。

restartControl： 整數，可選。

如果殘差與 retResidualNormsHistory 中記錄的最小值相比跳躍 2**restartControl 倍，則迭代重新開始。默認值為 restartControl=20 ，因此很少需要重新啟動以實現向後兼容性。

返回 ：：

lambda： 形狀為 (k, ) 的 ndarray 。

k 近似特征值數組。

v： 與 X.shape 形狀相同的 ndarray 。

k 近似特征向量數組。

lambdaHistory： ndarray，可選。

特征值曆史，如果retLambda曆史記錄是True.

ResidualNormsHistory： ndarray，可選。

剩餘範數的曆史，如果ret 剩餘規範曆史是True.

注意：

迭代循環最多運行 maxit=maxiter (如果 maxit=None 則為 20 次)迭代，如果滿足容差，則提前完成。 LOBPCG 打破了與先前版本的向後兼容性，現在以最佳精度返回迭代向量塊，而不是最後一個迭代向量，作為解決可能出現的分歧的方法。

如果X.dtype == np.float32和用戶提供的運算/乘法A,B， 和M都保留了np.float32數據類型，所有計算和輸出都在np.float32.

迭代曆史輸出的大小等於最佳(受 maxit 限製)迭代的數量加 3：初始、最終和後處理。

如果兩者都retLambda曆史記錄和ret 剩餘規範曆史是True，返回元組具有以下格式(lambda, V, lambda history, residual norms history).

下麵的n 表示矩陣大小，k 表示所需特征值的數量(最小或最大)。

LOBPCG 代碼內部解決了大小的特征問題3k在每次迭代中調用密集特征求解器八， 因此，如果k相比之下還不夠小n，調用 LOBPCG 代碼是沒有意義的。此外，如果調用 LOBPCG 算法5k > n，它可能會在內部崩潰，因此代碼調用標準函數八反而。不是那樣的nLOBPCG 應該很大才能發揮作用，但比率n / k應該很大。你用它來調用 LOBPCGk=1和n=10，它雖然有效n是小。該方法適用於非常大的n / k.

收斂速度主要取決於三個因子：

1. 尋找特征向量的初始近似值 X 的質量。如果不知道更好的選擇，則圍繞原始向量隨機分布的效果很好。

2. 所需特征值與其餘特征值的相對分離。人們可以改變k來改善分離。

3. 適當的預處理可縮小頻譜擴展。例如，對於大型杆振動測試問題(在測試目錄下)是ill-conditionedn，所以收斂會很慢，除非使用有效的預處理。對於這個特定問題，一個好的簡單預條件函數將是一個線性求解A，這很容易編碼，因為A是三對角線。

參考：

[1]

A. V. Knyazev (2001)，走向最優預條件本征求解器：局部最優塊預條件共軛梯度法。 SIAM 科學計算雜誌 23, no. 2，第 517-541 頁。 DOI:10.1137/S1064827500366124

[2]

A. V. Knyazev、I. Lashuk、M. E. Argentati 和 E. Ovchinnikov (2007)，在 hypre 和 PETSc 中阻止局部最優預條件特征值求解器 (BLOPEX)。 arXiv:0705.2626

[3]

A. V. Knyazev 的 C 和 MATLAB 實現：https://github.com/lobpcg/blopex

例子：

我們的第一個示例是簡約的 - 通過在沒有約束或預處理的情況下解決非廣義特征值問題 A x = lambda x 來找到對角矩陣的最大特征值。

>>> import numpy as np
>>> from scipy.sparse import spdiags
>>> from scipy.sparse.linalg import LinearOperator, aslinearoperator
>>> from scipy.sparse.linalg import lobpcg

方陣大小為

>>> n = 100

它的對角線條目是 1, …, 100 定義為

>>> vals = np.arange(1, n + 1).astype(np.int16)

此測試中的第一個強製輸入參數是稀疏對角矩陣A特征值問題A x = lambda x來解決。

>>> A = spdiags(vals, 0, n, n)
>>> A = A.astype(np.int16)
>>> A.toarray()
array([[  1,   0,   0, ...,   0,   0,   0],
       [  0,   2,   0, ...,   0,   0,   0],
       [  0,   0,   3, ...,   0,   0,   0],
       ...,
       [  0,   0,   0, ...,  98,   0,   0],
       [  0,   0,   0, ...,   0,  99,   0],
       [  0,   0,   0, ...,   0,   0, 100]], dtype=int16)

第二個強製輸入參數X是一個二維數組，行維數決定了請求的特征值的數量。X是對目標特征向量的初步猜測。X必須具有線性獨立的列。如果沒有可用的初始近似值，則隨機定向的向量通常效果最好，例如，分量通常分布在零附近或均勻分布在區間 [-1 1] 上。設置 dtype 的初始近似值np.float32強製所有迭代值都為 dtypenp.float32加快運行速度，同時仍然允許精確的特征值計算。

>>> k = 1
>>> rng = np.random.default_rng()
>>> X = rng.normal(size=(n, k))
>>> X = X.astype(np.float32)

>>> eigenvalues, _ = lobpcg(A, X, maxiter=60)
>>> eigenvalues
array([100.])
>>> eigenvalues.dtype
dtype('float32')

lobpcg隻需要訪問矩陣乘積A而不是矩陣本身。由於矩陣A在這個例子中是對角線，可以寫出矩陣乘積的函數A @ X使用對角線值vals例如，僅通過在 lambda-function 中進行廣播的逐元素乘法

>>> A_lambda = lambda X: vals[:, np.newaxis] * X

或常規函數

>>> def A_matmat(X):
...     return vals[:, np.newaxis] * X

並使用這些可調用對象之一的句柄作為輸入

>>> eigenvalues, _ = lobpcg(A_lambda, X, maxiter=60)
>>> eigenvalues
array([100.])
>>> eigenvalues, _ = lobpcg(A_matmat, X, maxiter=60)
>>> eigenvalues
array([100.])

傳統的可調用 LinearOperator 不再是必需的，但仍然支持作為 lobpcg 的輸入。明確指定 matmat=A_matmat 可以提高性能。

>>> A_lo = LinearOperator((n, n), matvec=A_matmat, matmat=A_matmat, dtype=np.int16)
>>> eigenvalues, _ = lobpcg(A_lo, X, maxiter=80)
>>> eigenvalues
array([100.])

效率最低的可調用選項是 aslinearoperator ：

>>> eigenvalues, _ = lobpcg(aslinearoperator(A), X, maxiter=80)
>>> eigenvalues
array([100.])

我們現在轉而計算指定的三個最小特征值

>>> k = 3
>>> X = np.random.default_rng().normal(size=(n, k))

和largest=False參數

>>> eigenvalues, _ = lobpcg(A, X, largest=False, maxiter=80)
>>> print(eigenvalues)  
[1. 2. 3.]

下一個示例說明計算同一矩陣的 3 個最小特征值A由函數句柄給出A_matmat但有約束和預處理。

約束 - 可選輸入參數是由特征向量必須正交的列向量組成的二維數組

>>> Y = np.eye(n, 3)

預調節器的作用相當於A在此示例中，但精度降低了np.float32盡管最初的X因此所有迭代和輸出都是完整的np.float64.

>>> inv_vals = 1./vals
>>> inv_vals = inv_vals.astype(np.float32)
>>> M = lambda X: inv_vals[:, np.newaxis] * X

現在讓我們首先求解矩陣 A 的特征值問題，無需進行 80 次迭代的預處理

>>> eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, largest=False, maxiter=80)
>>> eigenvalues
array([4., 5., 6.])
>>> eigenvalues.dtype
dtype('float64')

通過預處理，我們隻需要對同一個 X 進行 20 次迭代

>>> eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, M=M, largest=False, maxiter=20)
>>> eigenvalues
array([4., 5., 6.])

注意傳入的向量Y是 3 個最小特征值的特征向量。上麵返回的結果與這些結果正交。

初級矩陣A可能是不確定的，例如在轉移之後vals從 1, …, 100 到 -49, …, 50 乘以 50，我們仍然可以計算 3 個最小或最大特征值。

>>> vals = vals - 50
>>> X = rng.normal(size=(n, k))
>>> eigenvalues, _ = lobpcg(A_matmat, X, largest=False, maxiter=99)
>>> eigenvalues
array([-49., -48., -47.])
>>> eigenvalues, _ = lobpcg(A_matmat, X, largest=True, maxiter=99)
>>> eigenvalues
array([50., 49., 48.])