Ctrl+K
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 是一種用於大型實對稱和複 Hermitian 定義廣義特徵值問題的預處理特徵值求解器。

參數:
A{稀疏矩陣, ndarray, 線性算子, 可呼叫物件}

問題的 Hermitian 線性算子，通常由稀疏矩陣給出。常稱為「剛度矩陣」。

Xndarray，float32 或 float64

特徵向量 k 的初始近似值（非稀疏）。如果 A 的 shape=(n,n)，則 X 必須具有 shape=(n,k)。

B{稀疏矩陣, ndarray, 線性算子, 可呼叫物件}

可選。預設為 B = None，相當於單位矩陣。廣義特徵值問題中，如果存在，則為右手邊算子。常稱為「質量矩陣」。必須是 Hermitian 正定矩陣。

M{稀疏矩陣, ndarray, 線性算子, 可呼叫物件}

可選。預設為 M = None，相當於單位矩陣。旨在加速收斂的預處理器。

Yndarray，float32 或 float64，預設值：None

一個 n-by-sizeY 的 ndarray 約束，其中 sizeY < n。迭代將在 Y 的列空間的 B-正交補集中執行。如果存在，Y 必須是滿秩的。

tol純量，可選

預設值為 tol=n*sqrt(eps)。求解器的停止準則容忍度。

maxiterint，預設值：20

最大迭代次數。

largestbool，預設值：True

為 True 時，求解最大特徵值，否則求解最小特徵值。

verbosityLevelint，可選

預設值 verbosityLevel=0 無輸出。控制求解器的標準/螢幕輸出。

retLambdaHistorybool，預設值：False

是否返回迭代特徵值歷史記錄。

retResidualNormsHistorybool，預設值：False

是否返回殘差範數的迭代歷史記錄。

restartControlint，可選。

如果殘差跳躍到 retResidualNormsHistory 中記錄的最小值的 2**restartControl 倍，則迭代重新啟動。預設值為 restartControl=20，為了向後相容性，重新啟動很少發生。

返回:
lambda形狀為 (k, ) 的 ndarray。

近似特徵值 k 的陣列。

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

近似特徵向量 k 的陣列。

lambdaHistoryndarray，可選。

特徵值歷史記錄，如果 retLambdaHistory 為 True。

ResidualNormsHistoryndarray，可選。

殘差範數的歷史記錄，如果 retResidualNormsHistory 為 True。

註解

迭代迴圈最多執行 maxit=maxiter 次迭代（如果 maxit=None 則為 20 次），如果達到容忍度則提前結束。為了打破與先前版本的向後相容性，LOBPCG 現在返回具有最佳精度的迭代向量區塊，而不是最後一個迭代向量，作為可能的發散的補救措施。

如果 X.dtype == np.float32 且使用者提供的 A、B 和 M 的運算/乘法都保留 np.float32 資料類型，則所有計算和輸出都以 np.float32 進行。

迭代歷史記錄輸出的尺寸等於最佳迭代次數（受 maxit 限制）加 3：初始、最終和後處理。

如果 retLambdaHistory 和 retResidualNormsHistory 均為 True，則返回元組的格式如下 (lambda, V, lambda history, residual norms history)。

在以下內容中，n 表示矩陣大小，k 表示所需特徵值的數量（最小或最大）。

LOBPCG 程式碼在每次迭代中透過呼叫密集特徵值求解器 eigh 在內部求解大小為 3k 的特徵值問題，因此如果 k 相對於 n 不夠小，則呼叫 LOBPCG 程式碼是沒有意義的。此外，如果為 5k > n 呼叫 LOBPCG 演算法，則可能會在內部中斷，因此程式碼會改為呼叫標準函數 eigh。並不是說 n 應該很大 LOBPCG 才能運作，而是比率 n / k 應該很大。如果您使用 k=1 和 n=10 呼叫 LOBPCG，即使 n 很小，它也能運作。該方法適用於極大的 n / k。

收斂速度基本上取決於三個因素

  1. 尋找特徵向量的初始近似值 X 的品質。如果沒有更好的選擇，則隨機分佈在原點附近的向量效果良好。

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

  3. 適當的預處理以縮小頻譜擴展。例如，桿振動測試問題（在 tests 目錄下）對於大的 n 是病態的，因此收斂速度會很慢，除非使用有效的預處理。對於這個特定問題，一個好的簡單預處理器函數是對 A 進行線性求解，這很容易編碼，因為 A 是三對角矩陣。

參考文獻

[1]

A. V. Knyazev (2001), Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23, no. 2, pp. 517-541. DOI:10.1137/S1064827500366124

[2]

A. V. Knyazev, I. Lashuk, M. E. Argentati, and E. Ovchinnikov (2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre and PETSc. arXiv:0705.2626

[3]

A. V. Knyazev 的 C 和 MATLAB 實作：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 x = lambda x 的稀疏對角矩陣 A。

>>> 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]], shape=(100, 100), dtype=int16)

第二個強制輸入參數 X 是一個 2D 陣列，其列維度決定了請求的特徵值的數量。X 是目標特徵向量的初始猜測。X 必須具有線性獨立的列。如果沒有可用的初始近似值，隨機定向向量通常效果最佳，例如，元件通常分佈在零附近或均勻分佈在間隔 [-1 1] 上。將初始近似值設定為 dtype np.float32 會強制所有迭代值為 dtype np.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.], dtype=float32)

lobpcg 只需要存取與 A 的矩陣乘積，而不是矩陣本身。由於矩陣 A 在此範例中是對角矩陣，因此可以使用對角線值 vals 僅撰寫矩陣乘積 A @ X 的函數，例如，透過與 lambda 函數中的廣播進行元素級乘法

>>> 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.], dtype=float32)
>>> eigenvalues, _ = lobpcg(A_matmat, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)

傳統的可呼叫物件 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.], dtype=float32)

效率最低的可呼叫物件選項是 aslinearoperator

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

我們現在切換到計算三個最小的特徵值，指定

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

和 largest=False 參數

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

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

約束 - 可選的輸入參數是一個 2D 陣列，包含特徵向量必須與之正交的列向量

>>> 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 平移 50 到 -49, …, 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.])
在本頁面上