Julia 中 Arpack 对最小量级特征值的准确性问题

计算科学 特征值 收敛 朱莉娅
2021-12-17 07:50:48

根据 Julia 的 Arpack 包的文档(参见https://julialinearalgebra.github.io/Arpack.jl/stable/eigs/),我计算了以所需 CSC 格式编码的稀疏矩阵的一些最大和最小幅度特征值,并注意到主要关于最小幅度特征值的差异给我残差的数量级上,尽管所有这些特征值都被声称是收敛的.θAvθv∣∣/∣∣v107

我将结果与 Python 的 scipy.sparse.linalg 包(参见https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html)的计算进行了比较,它给出了适当的结果,即数量级的残差。105

稀疏矩阵来自 Groebner 基础理论中的麦考利矩阵,因此解释起来有点复杂,但我能够用随机生成的相同大小和相似条目的稀疏矩阵重现奇怪的行为。

这个朱莉娅代码

using Arpack, LinearAlgebra, SparseArrays

rf(n) = rand(-2^15+1:2^15, n)
S = sprand(229375, 16384, 0.0005, rf)
S = convert(SparseMatrixCSC{Float64,Int64},S)
A = transpose(S)*S
l = eigs(A, nev=4, which=:LM)
s = eigs(A, nev=4, which=:SM)

println("Number of converged eigenvalues of largest magnitude: ", l[3])
println("Number of converged eigenvalues of smallest magnitude: ", s[3])

# Tests

println( "Epsilon according to documentation: ", eps(real(eltype(A)))/2 )

# Largest Eigenvalue

x = 1
theta = l[1][x]
v = l[2][:,x]

println( "Residual of largest eigenvalue: ", norm(A * v - v * theta) )
println( "Convergence criterion for largest eigenvalue: ", (residual <= epsilon * max(cbrt(epsilon^2),abs(theta))) )
println( "Upper bound in convergence criterion for largest eigenvalue: ", epsilon * max(cbrt(epsilon^2),abs(theta)))

# Smallest Eigenvalue

y = 1
iota = s[1][x]
w = s[2][:,x]

println( "Residual of smallest eigenvalue: ", norm(A * w - w * iota) )
println( "Convergence criterion test for smallest eigenvalue: ", (residual <= epsilon * max(cbrt(epsilon^2),abs(iota))) )
println( "Upper bound in convergence criterion for smallest eigenvalue: ", epsilon * max(cbrt(epsilon^2),abs(iota)))

给出以下输出

Number of converged eigenvalues of largest magnitude: 4
Number of converged eigenvalues of smallest magnitude: 4
Epsilon according to documentation: 1.1102230246251565e-16
Residual of largest eigenvalue: 0.00018192639656053935
Convergence criterion for largest eigenvalue: false
Upper bound in convergence criterion for largest eigenvalue: 8.018644716606565e-6
Residual of smallest eigenvalue: 0.00011757473085532065
Convergence criterion test for smallest eigenvalue: false
Upper bound in convergence criterion for smallest eigenvalue: 2.056712217644155e-6

特征值看起来不错,但上述 Julia 的 Arpack 包文档中描述的收敛标准似乎与实际实现不同。

我正在使用 Julia 1.5.3

Installed Arpack_jll ─── v3.5.0+3
Installed OpenBLAS_jll ─ v0.3.9+5
Installed Arpack ─────── v0.5.1

在我上面提到的特定用例中,计算出的最小幅度特征值远没有任何用处。调整参数tol也没有帮助。有人可以解释这种奇怪的行为或帮助我评估正确的收敛标准并将计算调整为有用的残差吗?

更新

以下 Julia 代码提供了上述 BAD 残差的示例。

using Arpack, LinearAlgebra, SparseArrays

rr = [1,2,3,5,1,2,3,4,5,6,1,2,3,4,5,7,2,3,4,6,7,8,1,2,3,5,6,7,2,4,5,6,7,8,3,4,5,6,7,4,6,8]
cc = [1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,7,8,8,8]
dd = [4.02018713e8,1.42895766e8,1.91997182e8,-5.94015895e8,1.42895766e8,9.88669461e8,9.1763724e7,-3.932278e6,-2.3465949e8,-1.198600777e9,1.91997182e8,9.1763724e7,3.156467582e9,-6.25317095e8,-1.03130053e9,-2.317164234e9,-3.932278e6,-6.25317095e8,4.565270498e9,-4.93246362e8,6.25317095e8,-1.691847139e9,-5.94015895e8,-2.3465949e8,-1.03130053e9,2.883067709e9,-7.60572644e8,-1.257751284e9,-1.198600777e9,-4.93246362e8,-7.60572644e8,4.2947594e9,7.60572644e8,-4.9717864e8,-2.317164234e9,6.25317095e8,-1.257751284e9,7.60572644e8,3.574915518e9,-1.691847139e9,-4.9717864e8,5.815904688e9]

A = sparse(rr,cc,dd)

l = Arpack.eigs(A, nev=6, which=:LM)
s = Arpack.eigs(A, nev=6, which=:SM)

println("Number of converged largest magnitude Ritz values: ", l[3])
println("Number of converged smallest magnitude Ritz values: ", s[3])

epsilon = eps(real(eltype(A)))/2 # LAPACK's epsilon
tol = epsilon # tol default

println("Machine epsilon and tolerance default: ", epsilon)

l_res = [norm( A * l[2][:,i] - l[1][i] * l[2][:,i]) / norm( l[2][:,i] ) for i in 1:size(l[1])[1] ]
println("Residues of largest magnitude Ritz values: ", l_res)

s_res = [norm( A * s[2][:,i] - s[1][i] * s[2][:,i]) / norm( s[2][:,i] ) for i in 1:size(s[1])[1] ]
println("Residues of smallest magnitude Ritz values: ", s_res)

l_con = [(l_res[i] <= tol * max(cbrt(epsilon^2),abs(l[1][i]))) for i in 1:size(l[1])[1]]
println("Convergence criteria of largest magnitude Ritz values: ", l_con)

s_con = [(s_res[i] <= tol * max(cbrt(epsilon^2),abs(s[1][i]))) for i in 1:size(s[1])[1]]
println("Convergence criteria of smallest magnitude Ritz values: ", s_con)

输出是:

Number of converged largest magnitude Ritz values: 6
Number of converged smallest magnitude Ritz values: 6
Machine epsilon and tolerance default: 1.1102230246251565e-16
Residues of largest magnitude Ritz values: [2.5183513490019574e-6, 2.0857933315848784e-6, 2.3053979403750557e-6, 1.0323827311807143e-6, 4.182955279236495e-7, 1.3601799608960461e-6]
Residues of smallest magnitude Ritz values: [1.4473263034070623e-6, 3.3449964080033436e9, 2.50443508742856e9, 4.42611459230314e9, 1.1814710823584437e9, 6.968580232683431e8]
Convergence criteria of largest magnitude Ritz values: Bool[0, 0, 0, 0, 0, 0]
Convergence criteria of smallest magnitude Ritz values: Bool[0, 0, 0, 0, 0, 0]

如果您分析(声称的)最小幅度 Ritz 值,您会注意到相应的 Ritz 向量是零向量的估计值。此外,声称已经收敛的所有其他最小量级 Ritz 值都是虚假的。可以使用上面指定索引和数据的数组来用 scipy 实例化一个稀疏矩阵并计算可信赖的最小幅度 Ritz 值。

看起来很马车。还是有其他解释?

3个回答

看起来一切工作都比较顺利?您的矩阵的阶数为 1e10,因此 1e-4 的残差实际上接近机器精度。收敛准则确实被违反了,但不是很多;不知道那里发生了什么,你确定 Arpack 真的保证它还是尽力而为?令我惊讶的是,您在 python 中得到了不同的结果,因为它们都调用了同一个底层库(尽管这样做可能会有所不同)。eigs更直接地与 numpy's 相媲美eigs如果你愿意eigsh,你应该eigs(Symmetric(A))在朱莉娅做。

好的,所以在编辑新版本:这不是你的错,它是https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/87您可以自己手动调用 Arpack,或者更好地使用纯 julia 求解器(KrylovKit.jl 和 IterativeSolvers.jl 是不错的选择)

eigsh( ... sigma=0 )看起来也像 python-scipy arpack 中的错误。这是 arpack 与 numpy eigh(Lapack syevd,黄金标准)的比较,在您的矩阵/1e6 上运行:

sigma=0 --
evals numpy eigh:      [-8.23e-13 370 515 2.82e3 3.64e3 5.21e3 5.9e3 7.22e3]
evals arpack:     [-146 -5.68e-14 96.1 844 3.43e3 5.06e3 5.18e3] ❓

eigcheck AV - VΛ numpy: max 2.5e-12 at λ 3.64e3  [1.1e-12 1e-12 2.4e-12 1.9e-12 2.5e-12 1.4e-12 1.4e-12 1.8e-12]
eigcheck AV - VΛ arpack: max 1.7e3 at λ -146  [1.7e3 7.1e-13 9.3e2 1.3e3 3.1e2 7.5e2 1.9e2]

--------------------------------------------------------------------------------
arpack sigma=-1 agrees pretty well with numpy --
evals numpy eigh: [-8.23e-13 370 515 2.82e3 3.64e3 5.21e3 5.9e3 7.22e3]
evals arpack:     [7.73e-14 370 515 2.82e3 3.64e3 5.21e3 5.9e3]  # neig = n - 1

eigcheck AV - VΛ numpy: max 2.5e-12 at λ 3.64e3  [1.1e-12 1e-12 2.4e-12 1.9e-12 2.5e-12 1.4e-12 1.4e-12 1.8e-12]
eigcheck AV - VΛ arpack: max 4.2e-10 at λ 2.82e3  [6.4e-13 4.6e-11 7.5e-11 4.2e-10 4.8e-11 4.3e-11 3.8e-10]

任何人都可以尝试 julia arpacksigma=-1吗?

关于 Arpack 的注意事项

  • 使用eigsh( ... v0=np.ones ),而不是默认猜测特征向量v0= 随机。更好的是,v0从以前的运行或更粗的网格中获取。

  • 对于给定 sigma 附近的特征值,scipy arpack 调用密集/稀疏/线性运算符的不同内部求解器,默认为 LU/稀疏 LU/GMRES Ax = bA但是复杂的来回 arpack求解器是隐藏的,完全模糊。(可以制作一个带有回调的 LinOp 来跟踪这一点,或者使用例如 LGMRES 而不是 GMRES。)

  • 如果可以的话,将特征值接近 0 的云移到右半平面: sigma=-1在 上进行移位反转A + I,条件更好,然后移回。在这种情况下,sigma=-1有效。(也换LLT: 如果L是数值奇异的,LLT很可能有微小的负特征值,这对求解器来说很难。)

我在pybinding代码中偶然发现了这个:

   if sigma == 0:
        # eigsh can cause problems when sigma is exactly zero
        sigma = np.finfo(model.hamiltonian.dtype).eps