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, 线性算子, 可调用对象}: 问题的厄米线性算子，通常由一个稀疏矩阵给出。通常称为“刚度矩阵”。
Xndarray, float32 或 float64: 对 k 个特征向量的初始近似（非稀疏）。如果 A 的 shape=(n,n)，那么 X 必须有 shape=(n,k)。
B{稀疏矩阵, ndarray, 线性算子, 可调用对象}: 可选。默认情况下 B = None，这等同于单位矩阵。如果存在，这是广义特征问题中的右手边算子。通常称为“质量矩阵”。必须是厄米特正定矩阵。
M{稀疏矩阵, ndarray, 线性算子, 可调用对象}: 可选。默认情况下 M = None，这等同于恒等变换。预处理器旨在加速收敛。
Yndarray, float32 或 float64, 默认: None: 一个 n-by-sizeY 的约束 ndarray，其中 sizeY < n。迭代将在 Y 的列空间的 B-正交补中进行。如果存在，Y 必须是满秩的。
tol标量，可选: 默认值为 tol=n*sqrt(eps)。这是用于停止准则的求解器容差。
maxiterint, 默认值: 20: 最大迭代次数。
最大bool, 默认: True: 当为 True 时，求解最大的特征值，否则求解最小的特征值。
verbosityLevelint, 可选: 默认情况下 verbosityLevel=0 无输出。控制求解器的标准/屏幕输出。
retLambdaHistorybool, 默认值: False: 是否返回迭代特征值历史记录。
retResidualNormsHistorybool, 默认值: False: 是否返回残差范数的迭代历史。
重启控制int, 可选。: 如果残差的跳跃次数比 retResidualNormsHistory 中记录的最小值多 2**restartControl 次，迭代将重新开始。默认值为 restartControl=20，这使得重新启动对于向后兼容性来说是罕见的。

返回:

lambda : 形状为 (k, ) 的 ndarray。形状的 ndarray: k 个近似特征值的数组。
v : 与 X.shape 形状相同的 ndarray。与输入形状相同的 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 历史, 残差范数历史)。

在下面的内容中，n 表示矩阵的大小，k 表示所需特征值的数量（最小或最大）。

LOBPCG 代码在每次迭代中通过调用稠密特征求解器 eigh 内部求解大小为 3k 的特征问题，因此如果 k 相对于 n 不够小，调用 LOBPCG 代码就没有意义。此外，如果对 5k > n 调用 LOBPCG 算法，它可能会在内部中断，因此代码会调用标准函数 eigh 代替。LOBPCG 的工作并不要求 n 必须很大，而是比值 n / k 应该很大。如果你用 k=1 和 n=10 调用 LOBPCG，它虽然 n 很小但仍能工作。该方法旨在用于非常大的 n / k 。

收敛速度基本上取决于三个因素：

初始近似值 X 的质量对于寻找特征向量至关重要。如果不知道更好的选择，围绕原点随机分布的向量效果良好。
期望的特征值与剩余特征值之间的相对分离。可以通过改变 k 来改善分离效果。
适当的预处理以缩小谱宽。例如，杆振动测试问题（在测试目录下）对于较大的 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’s C and MATLAB implementations: 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]], dtype=int16)

第二个必需的输入参数 X 是一个二维数组，其行维度决定了请求的特征值的数量。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.])
>>> eigenvalues.dtype
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.])
>>> 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=90)
>>> print(eigenvalues)
[1. 2. 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 之后，我们仍然可以计算出最小的 3 个或最大的 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.])