Python SciPy linalg.lobpcg用法及代码示例

本文简要介绍 python 语言中 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，可调用对象}

问题的埃尔米特线性算子，通常由稀疏矩阵给出。通常称为“stiffness matrix”。

X： ndarray、float32 或 float64

初始近似值k特征向量(非稀疏)。如果A已shape=(n,n)然后X一定有shape=(n,k).

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

可选的。默认为 B = None ，相当于身份。广义特征值问题中的右侧运算符(如果存在)。通常称为“mass matrix”。必须是 Hermitian 正定的。

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

可选的。默认为 M = None ，相当于身份。预处理器旨在加速收敛。

Y： ndarray、float32 或 float64，默认值：无

一个n-by-sizeYndarray 约束条件sizeY < n。迭代将在B- column-space 的正交补集Y.Y如果存在，则必须是满级。

tol： 标量，可选

默认为 tol=n*sqrt(eps) 。求解器对停止标准的容差。

maxiter： 整数，默认：20

最大迭代次数。

largest： 布尔值，默认值：真

如果为 True，则求解最大的特征值，否则求解最小的特征值。

verbosityLevel： 整数，可选

默认情况下verbosityLevel=0没有输出。控制解算器标准/屏幕输出。

retLambdaHistory： 布尔值，默认值：假

是否返回迭代特征值历史。

retResidualNormsHistory： 布尔值，默认值：假

是否返回剩余规范的迭代历史。

restartControl： 整数，可选。

如果残差与 retResidualNormsHistory 中记录的最小值相比跳跃 2**restartControl 倍，则迭代重新开始。默认值为 restartControl=20 ，因此很少需要重新启动以实现向后兼容性。

返回 ：：

lambda： 形状为 (k, ) 的 ndarray 。

k 近似特征值数组。

v： 与 X.shape 形状相同的 ndarray 。

k 近似特征向量数组。

lambdaHistory： ndarray，可选。

特征值历史，如果retLambda历史记录是True.

ResidualNormsHistory： ndarray，可选。

剩余范数的历史，如果ret 剩余规范历史是True.

注意：

迭代循环最多运行 maxit=maxiter (如果 maxit=None 则为 20 次)迭代，如果满足容差，则提前完成。 LOBPCG 打破了与先前版本的向后兼容性，现在以最佳精度返回迭代向量块，而不是最后一个迭代向量，作为解决可能出现的分歧的方法。

如果X.dtype == np.float32和用户提供的运算/乘法A,B， 和M都保留了np.float32数据类型，所有计算和输出都在np.float32.

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

如果两者都retLambda历史记录和ret 剩余规范历史是True，返回元组具有以下格式(lambda, V, lambda history, residual norms history).

下面的n 表示矩阵大小，k 表示所需特征值的数量(最小或最大)。

LOBPCG 代码内部解决了大小的特征问题3k在每次迭代中调用密集特征求解器八， 因此，如果k相比之下还不够小n，调用 LOBPCG 代码是没有意义的。此外，如果调用 LOBPCG 算法5k > n，它可能会在内部崩溃，因此代码调用标准函数八反而。不是那样的nLOBPCG 应该很大才能发挥作用，但比率n / k应该很大。你用它来调用 LOBPCGk=1和n=10，它虽然有效n是小。该方法适用于非常大的n / k.

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

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

2. 所需特征值与其余特征值的相对分离。人们可以改变k来改善分离。

3. 适当的预处理可缩小频谱扩展。例如，对于大型杆振动测试问题(在测试目录下)是ill-conditionedn，所以收敛会很慢，除非使用有效的预处理。对于这个特定问题，一个好的简单预条件函数将是一个线性求解A，这很容易编码，因为A是三对角线。

参考：

[1]

A. V. Knyazev (2001)，走向最优预条件本征求解器：局部最优块预条件共轭梯度法。 SIAM 科学计算杂志 23, no. 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 实现：https://github.com/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特征值问题A x = lambda x来解决。

>>> 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强制所有迭代值都为 dtypenp.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在这个例子中是对角线，可以写出矩阵乘积的函数A @ X使用对角线值vals例如，仅通过在 lambda-function 中进行广播的逐元素乘法

>>> 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=80)
>>> print(eigenvalues)  
[1. 2. 3.]

下一个示例说明计算同一矩阵的 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 乘以 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.])