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, LinearOperator, 可调用对象}

问题的 Hermitian 线性算子,通常由稀疏矩阵给出。通常称为“刚度矩阵”。

Xndarray, float32 或 float64

k 个特征向量的初始近似值(非稀疏)。如果 Ashape=(n,n),则 X 必须具有 shape=(n,k)

B{稀疏矩阵, ndarray, LinearOperator, 可调用对象}

可选。默认情况下,B = None,这等效于单位矩阵。如果存在,则为广义特征值问题中的右手侧算子。通常称为“质量矩阵”。必须是 Hermitian 正定矩阵。

M{稀疏矩阵, ndarray, LinearOperator, 可调用对象}

可选。默认情况下,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 个近似特征值的数组。

vX.shape 形状相同的 ndarray。

k 个近似特征向量的数组。

lambdaHistoryndarray,可选。

如果 retLambdaHistoryTrue,则为特征值历史记录。

ResidualNormsHistoryndarray,可选。

如果 retResidualNormsHistoryTrue,则为残差范数历史记录。

注释

迭代循环最多运行 maxit=maxiter (如果 maxit=None 则为 20) 次迭代,如果满足容差则会提前结束。打破与之前版本的向后兼容性,LOBPCG 现在返回具有最佳精度的迭代向量块,而不是最后一次迭代的向量块,作为对可能发散的修复。

如果 X.dtype == np.float32 并且用户提供的 ABM 的操作/乘法都保留了 np.float32 数据类型,则所有计算和输出均采用 np.float32

迭代历史记录输出的大小等于最佳(受 maxit 限制)迭代次数加 3:初始、最终和后处理。

如果 retLambdaHistoryretResidualNormsHistory 均为 True,则返回元组的格式如下 (lambda, V, lambda history, residual norms history)

在以下内容中,n 表示矩阵大小,k 表示所需的特征值数量(最小或最大)。

LOBPCG 代码在每次迭代时通过调用密集特征值求解器 eigh 在内部求解大小为 3k 的特征值问题,因此如果 kn 相比不够小,则调用 LOBPCG 代码没有意义。此外,如果对 5k > n 调用 LOBPCG 算法,则可能会在内部中断,因此代码会改为调用标准函数 eigh。并不是说 n 应该很大才能使 LOBPCG 工作,而是比率 n / k 应该很大。如果你用 k=1n=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 是一个二维数组,其行维度决定了请求的特征值的数量。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)

传统的 callable 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_matmat 的情况下,计算同一矩阵 A 的 3 个最小特征值,但带有约束和预处理。

约束 - 可选的输入参数是一个二维数组,其中包含特征向量必须与之正交的列向量

>>> 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.])