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

问题的厄米特线性算子，通常由稀疏矩阵给出。通常被称为“刚度矩阵”。

Xndarray，float32 或 float64

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

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

可选。默认情况下 B = None，相当于恒等式。如果存在广义特征值问题中的右手侧算子。通常被称为“质量矩阵”。必须是厄米特正定。

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，为了向后兼容，使重新启动很少发生。

返回:

lambdandarray，形状为 (k, )。

k 个近似特征值的数组。

vndarray，与 X.shape 相同的形状。

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 history, residual norms history)。

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

LOBPCG 代码在每次迭代中通过调用密集特征值求解器 eigh 内部解决大小为 3k 的特征值问题，因此如果 k 与 n 相比不够小，则调用 LOBPCG 代码没有意义。此外，如果调用 5k > n 的 LOBPCG 算法，则可能会在内部中断，因此代码会改为调用标准函数 eigh。这并不是说 n 必须很大才能使 LOBPCG 工作，而是说 n / k 的比率应该很大。如果使用 k=1 和 n=10 调用 LOBPCG，则它可以工作，尽管 n 很小。该方法适用于极大的 n / k。

收敛速度主要取决于三个因素

1. 对寻求特征向量 X 的初始近似值的质量。如果不知道更好的选择，则在原点周围随机分布的向量效果很好。

2. 所需特征值与其余特征值之间的相对分离度。可以更改 k 以改善分离度。

3. 适当的预处理以缩小光谱范围。例如，杆振动测试问题（在测试目录下）对于较大的 n 来说是病态的，因此除非使用有效的预处理，否则收敛速度会很慢。对于此特定问题，一个好的简单预处理函数将是对 A 进行线性求解，这很容易编码，因为 A 是三对角的。

参考文献

[1]

A. V. Knyazev (2001)，走向最佳预处理特征值求解器：局部最优块预处理共轭梯度法。SIAM 科学计算杂志 23，第 2 期，第 517-541 页。 DOI:10.1137/S1064827500366124

[2]

A. V. Knyazev、I. Lashuk、M. E. Argentati 和 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 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 的 3 个最小特征值，该矩阵由函数句柄 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 个最小或最大特征值。

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

返回在新标签页中打开