根据 Julia 的 Arpack 包的文档(参见https://julialinearalgebra.github.io/Arpack.jl/stable/eigs/),我计算了以所需 CSC 格式编码的稀疏矩阵的一些最大和最小幅度特征值,并注意到主要关于最小幅度特征值的差异给我残差在的数量级上,尽管所有这些特征值都被声称是收敛的.
我将结果与 Python 的 scipy.sparse.linalg 包(参见https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html)的计算进行了比较,它给出了适当的结果,即数量级的残差。
稀疏矩阵来自 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 值。
看起来很马车。还是有其他解释?