稀疏特征值问题与ARPACK#

简介#

ARPACK [1] 是一个Fortran包，它提供了快速找到大型稀疏矩阵的几个特征值/特征向量的例程。为了找到这些解，它只需要通过矩阵进行左乘。这个操作是通过一个*反向通信*接口执行的。这种结构的结果是，ARPACK能够找到任何将向量映射到向量的线性函数的特征值和特征向量。

ARPACK提供的所有功能都包含在两个高级接口中：scipy.sparse.linalg.eigs 和 scipy.sparse.linalg.eigsh。eigs 提供了求解实数或复数非对称方阵的特征值/特征向量的接口，而 eigsh 提供了求解实对称或复数厄米矩阵的接口。

基本功能#

ARPACK可以求解以下形式的标准特征值问题

\[A \mathbf{x} = \lambda \mathbf{x}\]

或以下形式的一般特征值问题

\[A \mathbf{x} = \lambda M \mathbf{x}.\]

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' 可能会导致执行时间缓慢和/或异常结果。 更好的方法是使用*移位-反演模式*。

移位-反演模式#

移位-反演模式依赖于以下观察结果。对于广义特征值问题

\[A \mathbf{x} = \lambda M \mathbf{x},\]

可以证明

\[(A - \sigma M)^{-1} M \mathbf{x} = \nu \mathbf{x},\]

其中

\[\nu = \frac{1}{\lambda - \sigma}.\]

示例#

假设您想找到一个大矩阵的最小和最大特征值及其对应的特征向量。 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)  # 创建一个对称矩阵

我们现在有一个对称矩阵 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.],       # 可能会有所不同（符号）
       [ 0.  1.  0.],
       [-0.  0. -1.]])

结果如预期。ARPACK 恢复了所需的特征值，并且它们与之前已知的结果相匹配。此外，特征向量是正交的，正如我们所期望的那样。现在，让我们尝试求解最小量级的特征值：

>>> evals_small, evecs_small = eigsh(X, 3, which='SM')
Traceback (most recent call last):       # 可能会有所不同（收敛性）
...
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],    # 可能会有所不同（符号）
       [-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.],           # 可能会有所不同（符号）
       [-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.],    # 可能会有所不同（符号）
       [ 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.],     # 可能会有所不同（符号）
       [-0. -0.  1.],
       [ 1.  0.  0.]]

特征值以不同的顺序出现，但它们都在这里。 请注意，移位-反演模式需要内部求解矩阵的逆。 这由 eigsh 和 eigs 自动处理，但用户也可以指定该操作。 有关详细信息，请参阅 scipy.sparse.linalg.eigsh 和 scipy.sparse.linalg.eigs 的文档字符串。

使用 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.],     # 可能会有所不同（符号）
       [-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 # 消除边缘效应
>>> 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]) # 可能会有所不同

需要注意的是，该算子的特征值均为虚数。此外，scipy.sparse.linalg.eigs 中的关键字 which='LI' 生成了具有最大绝对虚部的特征值（包括正负）。同样，在 PyLops 库中有一个更高级的一阶导数算子实现，名为 FirstDerivative 算子。

参考文献#

On this page