使用 ARPACK 的稀疏特徵值問題#
簡介#
ARPACK [1] 是一個 Fortran 套件，提供快速尋找大型稀疏矩陣的一些特徵值/特徵向量的例程。為了找到這些解，它只需要與相關矩陣進行左乘。此操作透過反向通訊介面執行。這種結構的結果是，ARPACK 能夠找到將向量映射到向量的任何線性函數的特徵值和特徵向量。
ARPACK 中提供的所有功能都包含在兩個高階介面 scipy.sparse.linalg.eigs 和 scipy.sparse.linalg.eigsh 中。eigs 提供用於尋找實數或複數非對稱方陣的特徵值/向量的介面，而 eigsh 提供用於實數對稱或複數 Hermitian 矩陣的介面。
基本功能#
ARPACK 可以解決以下形式的標準特徵值問題
或以下形式的廣義特徵值問題
ARPACK 的強大之處在於它只能計算特徵值/特徵向量對的指定子集。這是透過關鍵字 which 完成的。which 的以下值可用
which = 'LM': 具有最大量值的特徵值 (eigs,eigsh)，也就是複數歐幾里得範數中最大的特徵值。which = 'SM': 具有最小量值的特徵值 (eigs,eigsh)，也就是複數歐幾里得範數中最小的特徵值。which = 'LR': 具有最大實部的特徵值 (eigs)。which = 'SR': 具有最小實部的特徵值 (eigs)。which = 'LI': 具有最大虛部的特徵值 (eigs)。which = 'SI': 具有最小虛部的特徵值 (eigs)。which = 'LA': 具有最大代數值的特徵值 (eigsh)，也就是包含任何負號的最大特徵值。which = 'SA': 具有最小代數值的特徵值 (eigsh)，也就是包含任何負號的最小特徵值。which = 'BE': 來自頻譜兩端的特徵值 (eigsh)。
請注意，ARPACK 通常更擅長尋找極端特徵值，也就是具有大量值的特徵值。特別是，使用 which = 'SM' 可能會導致執行時間緩慢和/或異常結果。更好的方法是使用移位-反轉模式。
移位-反轉模式#
移位-反轉模式依賴於以下觀察。對於廣義特徵值問題
可以證明
其中
範例#
假設您想找到大型矩陣的最小和最大特徵值以及相應的特徵向量。ARPACK 可以處理多種形式的輸入：稠密矩陣，例如 numpy.ndarray 實例、稀疏矩陣，例如 scipy.sparse.csr_matrix，或從 scipy.sparse.linalg.LinearOperator 衍生的通用線性運算子。對於此範例，為了簡單起見，我們將建構一個對稱正定矩陣。
>>> import numpy as np
>>> from scipy.linalg import eig, eigh
>>> from scipy.sparse.linalg import eigs, eigsh
>>> np.set_printoptions(suppress=True)
>>> rng = np.random.default_rng()
>>>
>>> X = rng.random((100, 100)) - 0.5
>>> X = np.dot(X, X.T) # create a symmetric matrix
現在我們有一個對稱矩陣 X，用於測試例程。首先，使用 eigh 計算標準特徵值分解
>>> evals_all, evecs_all = eigh(X)
隨著 X 的維度增長，此例程會變得非常慢。特別是，如果只需要幾個特徵向量和特徵值，ARPACK 可能是一個更好的選擇。首先，讓我們計算 X 的最大特徵值 (which = 'LM') 並將它們與已知結果進行比較
>>> evals_large, evecs_large = eigsh(X, 3, which='LM')
>>> print(evals_all[-3:])
[29.22435321 30.05590784 30.58591252]
>>> print(evals_large)
[29.22435321 30.05590784 30.58591252]
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1. 0. 0.], # may vary (signs)
[ 0. 1. 0.],
[-0. 0. -1.]])
結果如預期。ARPACK 恢復了所需的特徵值，並且它們與先前已知的結果相符。此外，特徵向量是正交的，正如我們所預期的。現在，讓我們嘗試求解最小量值的特徵值
>>> evals_small, evecs_small = eigsh(X, 3, which='SM')
Traceback (most recent call last): # may vary (convergence)
...
scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence:
ARPACK error -1: No convergence (1001 iterations, 0/3 eigenvectors converged)
糟糕。我們看到，如上所述，ARPACK 並不完全擅長尋找小特徵值。有幾種方法可以解決此問題。我們可以增加容差 (tol) 以加快收斂速度
>>> evals_small, evecs_small = eigsh(X, 3, which='SM', tol=1E-2)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 0.99999999 0.00000024 -0.00000049], # may vary (signs)
[-0.00000023 0.99999999 0.00000056],
[ 0.00000031 -0.00000037 0.99999852]])
這可行，但我們會失去結果的精確度。另一個選項是將最大迭代次數 (maxiter) 從 1000 增加到 5000
>>> evals_small, evecs_small = eigsh(X, 3, which='SM', maxiter=5000)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1. 0. 0.], # may vary (signs)
[-0. 1. 0.],
[ 0. 0. -1.]])
我們得到了我們希望的結果，但計算時間要長得多。幸運的是，ARPACK 包含一種模式，可以快速確定非外部特徵值：移位-反轉模式。如上所述，此模式涉及將特徵值問題轉換為具有不同特徵值的等效問題。在這種情況下，我們希望找到接近零的特徵值，因此我們將選擇 sigma = 0。然後，轉換後的特徵值將滿足 \(\nu = 1/(\lambda - \sigma) = 1/\lambda\)，因此我們的小特徵值 \(\lambda\) 變成大特徵值 \(\nu\)。
>>> evals_small, evecs_small = eigsh(X, 3, sigma=0, which='LM')
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1. 0. 0.], # may vary (signs)
[ 0. -1. -0.],
[-0. -0. 1.]])
我們得到了我們希望的結果，並且計算時間少得多。請注意，從 \(\nu \to \lambda\) 的轉換完全在後台進行。使用者無需擔心細節。
移位-反轉模式不僅提供了一種快速獲得一些小特徵值的方法。假設您希望找到內部特徵值和特徵向量，例如，最接近 \(\lambda = 1\) 的那些。只需設定 sigma = 1，ARPACK 將處理其餘部分
>>> evals_mid, evecs_mid = eigsh(X, 3, sigma=1, which='LM')
>>> i_sort = np.argsort(abs(1. / (1 - evals_all)))[-3:]
>>> evals_all[i_sort]
array([0.94164107, 1.05464515, 0.99090277])
>>> evals_mid
array([0.94164107, 0.99090277, 1.05464515])
>>> print(np.dot(evecs_mid.T, evecs_all[:,i_sort]))
array([[-0. 1. 0.], # may vary (signs)
[-0. -0. 1.],
[ 1. 0. 0.]]
特徵值以不同的順序出現，但它們都在那裡。請注意，移位-反轉模式需要矩陣反轉的內部解。這由 eigsh 和 eigs 自動處理，但該操作也可以由使用者指定。有關詳細資訊，請參閱 scipy.sparse.linalg.eigsh 和 eigs 的 docstring。
使用 LinearOperator#
現在我們考慮您想要避免建立稠密矩陣並改用 scipy.sparse.linalg.LinearOperator 的情況。我們的第一個線性運算子應用輸入向量與使用者提供給運算子本身的向量 \(\mathbf{d}\) 之間的逐元素乘法。此運算子模擬對角矩陣，其對角線上的元素為 \(\mathbf{d}\)，並且它具有主要優點，即正向和伴隨運算是簡單的逐元素乘法，而不是矩陣-向量乘法。對於對角矩陣，我們預期特徵值等於對角線上的元素，在本例中為 \(\mathbf{d}\)。將使用 eigsh 獲得的特徵值和特徵向量與應用於稠密矩陣時使用 eigh 獲得的特徵值和特徵向量進行比較
>>> from scipy.sparse.linalg import LinearOperator
>>> class Diagonal(LinearOperator):
... def __init__(self, diag, dtype='float32'):
... self.diag = diag
... self.shape = (len(self.diag), len(self.diag))
... self.dtype = np.dtype(dtype)
... def _matvec(self, x):
... return self.diag*x
... def _rmatvec(self, x):
... return self.diag*x
>>> N = 100
>>> rng = np.random.default_rng()
>>> d = rng.normal(0, 1, N).astype(np.float64)
>>> D = np.diag(d)
>>> Dop = Diagonal(d, dtype=np.float64)
>>> evals_all, evecs_all = eigh(D)
>>> evals_large, evecs_large = eigsh(Dop, 3, which='LA', maxiter=1e3)
>>> evals_all[-3:]
array([1.53092498, 1.77243671, 2.00582508])
>>> evals_large
array([1.53092498, 1.77243671, 2.00582508])
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1. 0. 0.], # may vary (signs)
[-0. -1. 0.],
[ 0. 0. -1.]]
在本例中，我們建立了一個快速而簡單的 Diagonal 運算子。外部函式庫 PyLops 在 Diagonal 運算子以及其他幾個運算子中提供了類似的功能。
最後，我們考慮一個線性運算子，該運算子模擬一階導數模板的應用。在這種情況下，運算子等效於實數非對稱矩陣。我們再次將估計的特徵值和特徵向量與將相同一階導數應用於輸入訊號的稠密矩陣的特徵值和特徵向量進行比較
>>> class FirstDerivative(LinearOperator):
... def __init__(self, N, dtype='float32'):
... self.N = N
... self.shape = (self.N, self.N)
... self.dtype = np.dtype(dtype)
... def _matvec(self, x):
... y = np.zeros(self.N, self.dtype)
... y[1:-1] = (0.5*x[2:]-0.5*x[0:-2])
... return y
... def _rmatvec(self, x):
... y = np.zeros(self.N, self.dtype)
... y[0:-2] = y[0:-2] - (0.5*x[1:-1])
... y[2:] = y[2:] + (0.5*x[1:-1])
... return y
>>> N = 21
>>> D = np.diag(0.5*np.ones(N-1), k=1) - np.diag(0.5*np.ones(N-1), k=-1)
>>> D[0] = D[-1] = 0 # take away edge effects
>>> Dop = FirstDerivative(N, dtype=np.float64)
>>> evals_all, evecs_all = eig(D)
>>> evals_large, evecs_large = eigs(Dop, 4, which='LI')
>>> evals_all_imag = evals_all.imag
>>> isort_imag = np.argsort(np.abs(evals_all_imag))
>>> evals_all_imag = evals_all_imag[isort_imag]
>>> evals_large_imag = evals_large.imag
>>> isort_imag = np.argsort(np.abs(evals_large_imag))
>>> evals_large_imag = evals_large_imag[isort_imag]
>>> evals_all_imag[-4:]
array([-0.95105652, 0.95105652, -0.98768834, 0.98768834])
>>> evals_large_imag
array([0.95105652, -0.95105652, 0.98768834, -0.98768834]) # may vary
請注意，此運算子的特徵值都是虛數。此外，scipy.sparse.linalg.eigs 的關鍵字 which='LI' 產生具有最大絕對虛部的特徵值（正數和負數）。同樣，PyLops 函式庫中以 FirstDerivative 運算子的名稱提供了更進階的一階導數運算子實作。