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 个特征向量（非稀疏）。如果 A 的 shape=(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 个近似特征值的数组。

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 历史记录, 残差 范数 历史记录)。

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

LOBPCG 代码在每次迭代时通过调用密集特征值求解器 eigh 在内部求解大小为 3k 的特征值问题，因此如果 k 相对于 n 不够小，则调用 LOBPCG 代码是没有意义的。此外，如果一个人调用 LOBPCG 算法求解 5k > n，则很可能会在内部中断，因此代码会调用标准函数 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)，迈向最佳预处理特征值求解器：局部最优块预处理共轭梯度法。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)，hypre 和 PETSc 中的块局部最优预处理特征值求解器 (BLOPEX)。arXiv:0705.2626

[3]

A. V. Knyazev 的 C 和 MATLAB 实现：lobpcg/blopex

示例

在您的浏览器中尝试一下！

我们的第一个例子是最小的 - 通过求解非广义特征值问题 A x = lambda x，在没有约束或预处理的情况下找到对角矩阵的最大特征值。

>>> import numpy as np
>>> from scipy.sparse import diags_array
>>> 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 = diags_array(vals, offsets=0, shape=(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)

传统的 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)

效率最低的 callable 选项是 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 个最小特征值，但具有约束和预处理。

约束 - 可选的输入参数是一个 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.])

返回在新标签页中打开