对角化病态对称矩阵。好的特征值但坏的特征向量答案

【问题标题】：Diagonalization ill-conditioned symmetric matrix. Good eigenvalues but bad eigenvectors对角化病态对称矩阵。好的特征值但坏的特征向量
【发布时间】：2018-09-07 21:33:16
【问题描述】：

我需要对具有非常小的值的病态稀疏矩阵进行对角化。我不得不说带有 LAPACK 的 C++ 能够正确完成，所以我希望 Julia 也能做到。问题很简单eigvals(Matrix) 给出了正确的频谱，但eigen(Matrix) 给出了一个糟糕的频谱，因此，糟糕的特征向量。所以我的问题是：有什么方法可以正确计算特征向量？我在这里粘贴一个最小的完整示例：

using LinearAlgebra
using SparseArrays

hops=[-1.0e-60, -1.0e-55, -1.0e-50, -1.0e-45, -1.0e-40, -1.0e-35, -1.0e-30, -1.0e-25, -1.0e-20, -1.0e-15, -1.0e-10, -1.0e-5, -0.00316228, -1.0e-5, -1.0e-10, -1.0e-15, -1.0e-20, -1.0e-25, -1.0e-30, -1.0e-35, -1.0e-40, -1.0e-45, -1.0e-50, -1.0e-55, -1.0e-60]

ham=diagm(-1 => hops, 1=>hops)
ham_dense=Array(ham)
s1=eigvals(ham_dense)
s2,basis=eigen(ham_dense)
println(s1)
println(s2)

具体来说，我们有 eigvals 给：

[-0.00316231, -3.16228e-8, -3.16228e-13, -3.16228e-18, -3.16228e-23, -3.16228e-28, -3.16228e-33, -3.16228e-38, -3.16228e-43, -3.16228e-48, -3.16228e-53, -3.16228e-58, -3.16225e-63, 3.16225e -63、3.16228e-58、3.16228e-53、3.16228e-48、3.16228e-43、3.16228e-38、3.16228e-33、3.16228e-28、3.16228e-23、3.16228e-18、3.126 -13, 3.16228e-8, 0.00316231]

而eigen得到的频谱为：

[-0.00316231, -3.16228e-8, -3.16225e-13, -2.09351e-18, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.24468e-18, 3.1623e-13, 3.16228e-8, 0.00316231]

非常感谢您。

【问题讨论】：

标签： arrays julia linear-algebra

【解决方案1】：

问题的原因不是 Julia，而是 LAPACK。这些是最终在您的问题中发生的对 LAPACK 的调用：

julia> A = Symmetric(ham_dense);

julia> LAPACK.syevr!('N', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)
([-0.00316231, -3.16228e-8, -3.16228e-13, -3.16228e-18, -3.16228e-23, -3.16228e-28, -3.16228e-33, -3.16228e-38, -3.16228e-43, -3.16228e-48  …  3.16228e-48, 3.16228e-43, 3.16228e-38, 3.16228e-33, 3.16228e-28, 3.16228e-23, 3.16228e-18, 3.16228e-13, 3.16228e-8, 0.00316231], Array{Float64}(26,0))

julia> LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)
([-0.00316231, -3.16228e-8, -3.16226e-13, -2.01616e-18, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.33528e-18, 3.1623e-13, 3.16228e-8, 0.00316231], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

here 提供了有关所用算法的更多详细信息，显然在您的问题中，'V' 和 'N' 选项之间的差异很重要。最可能的原因是矩阵中最大和最小绝对非零值的相对差异为3.1622799999999996e57，因此可能会出现舍入误差。

要获得更准确的结果，您可以尝试 https://github.com/andreasnoack/GenericLinearAlgebra.jl 解决此类问题（它还没有 100% 的覆盖率并且是实验性的；您必须直接从 GitHub 安装它）：

julia> using GenericLinearAlgebra

julia> eigvals!(BigFloat.(ham))
26-element Array{Complex{BigFloat},1}:
  3.162311622436982597307210858940985344637820907597401918790008106376535395523792e-03 + 0.0im
 -3.162311622436982597307210858940985344637820907597401918790008106376535395523792e-03 - 0.0im
  3.162275320748298914345541008173998213886552122046290838685673282360067640110347e-08 + 0.0im
 -3.162275320748298914345541008173998213886552122046290838685673282360067640110347e-08 - 0.0im
  3.162279999906412335597731744982886550300231679921373378044484886374018562948466e-13 + 0.0im
 -3.162279999906412335597731744982886550300231679921373378044484886374018562948466e-13 - 0.0im
  3.162275320432081129857004982100630443653566381027866582713636819787519332736638e-18 + 0.0im
 -3.162275320432081129857004982100630443653566381027866582713636819787519332736638e-18 - 0.0im
   3.16227999990640833509930880688425213299040187766949574394524744644634876014558e-23 + 0.0im
  -3.16227999990640833509930880688425213299040187766949574394524744644634876014558e-23 - 0.0im
  3.162275320432081720151694936640991555004802493020434614166313056473370703548901e-28 + 0.0im
 -3.162275320432081720151694936640991555004802493020434614166313056473370703548901e-28 - 0.0im
  3.162279999906408618697844989393166731890284189322183903857450691187082378309407e-33 + 0.0im
 -3.162279999906408618697844989393166731890284189322183903857450691187082378309407e-33 - 0.0im
  3.162275320432081242789139430274227532672254451244152725835934892025795336479976e-38 + 0.0im
 -3.162275320432081242789139430274227532672254451244152725835934892025795336479976e-38 - 0.0im
  3.162279999906408121819532864964337223950046982981611768712594199335770669717062e-43 + 0.0im
 -3.162279999906408121819532864964337223950046982981611768712594199335770669717062e-43 - 0.0im
  3.162275320432081589442009096101640817970407333051878158079593749435340884536119e-48 + 0.0im
 -3.162275320432081589442009096101640817970407333051878158079593749435340884536119e-48 - 0.0im
   3.16227999990640827051609368034343266099945811896062048249442816603909306889103e-53 + 0.0im
  -3.16227999990640827051609368034343266099945811896062048249442816603909306889103e-53 - 0.0im
  3.162275320432081506938731037100931333563539294943227057058654130546144762012267e-58 + 0.0im
 -3.162275320432081506938731037100931333563539294943227057058654130546144762012267e-58 - 0.0im
  3.162248377469425472537552742972431066440507090082172984570144491979906956991795e-63 + 0.0im
 -3.162248377469425472537552742972431066440507090082172984570144491979906956991795e-63 - 0.0im

对于Float64，它接近于eigvals。

【讨论】：

你好。非常感谢您的回答。确实是lapack的问题。在 C++ 中，我使用 dsyevx。我需要使用那个特定的功能。 ccall 可以吗？我做不到
应该可以拨打ccall，但我没试过。

【解决方案2】：

也许我理解错了，但是特征值和特征向量是相同的（直到浮点精度）：

julia> ev = eigvals(ham_dense);

julia> evec = eigvecs(ham_dense);

julia> a,b = eigen(ham_dense);

julia> isapprox(a, ev)
true

julia> isapprox(b, evec)
true

julia> ev - a
26-element Array{Float64,1}:
 -7.806255641895632e-18
 -4.9900762330669964e-18
 -2.435853087382333e-18
 -1.146119995342351e-18
 -3.162279999906409e-23
 -3.1622753204320815e-28
 -3.1622799999064074e-33
 -3.1622753204320813e-38
 -3.1622799999064078e-43
 -3.162275320432081e-48
 -3.162279999906408e-53
 -3.1622753204320827e-58
 -3.162248377469362e-63
  3.1622483774694875e-63
  3.1622753204320827e-58
  3.162279999906408e-53
  3.162275320432081e-48
  3.1622799999064078e-43
  3.1622753204320813e-38
  3.1622799999064074e-33
  3.1622753204320815e-28
  3.1622799999064087e-23
 -1.1730014662417663e-18
 -2.4358773191193703e-18
 -4.9901225551812994e-18
 -4.336808689942018e-19

【讨论】：