对角化病态对称矩阵。良好的特征值，但不好的特征向量

Question

我需要将具有非常小值的病态稀疏矩阵对角化。我不得不说带有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.16228e-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]

非常感谢你提前。

Answer 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

这是接近eigvals为Float64。

Answer 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