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
,则为残差范数历史记录。
- lambda形状为
注释
迭代循环最多运行
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
。收敛速度基本上取决于三个因素
初始近似值 X 对所寻求的特征向量的质量。如果不知道更好的选择,则随机分布在原点周围的向量效果很好。
所需特征值与其余特征值的相对分离。可以更改
k
以改善分离。适当的预处理以缩小频谱扩展。例如,杆振动测试问题(在 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
会强制所有迭代值都为 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.], 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.])