Диагонализация плохо обусловленной симметричной матрицы. Хорошие собственные значения, но плохие собственные векторы - PullRequest
Новая
       22

Диагонализация плохо обусловленной симметричной матрицы. Хорошие собственные значения, но плохие собственные векторы

0 голосов
/ 07 сентября 2018

Мне нужно диагонализировать плохо обусловленную разреженную матрицу с очень маленькими значениями. Я должен сказать, что C ++ с LAPACK может делать правильно, поэтому я надеюсь, что Джулия может сделать это также. Проблема проста: 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]

Заранее большое спасибо.

Ответы [ 2 ]

0 голосов
/ 07 сентября 2018

Причина вашей проблемы не в Джулии, а в ЛАПАКЕ. Это вызовы 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])

Более подробная информация об используемом алгоритме дана здесь и, очевидно, в вашей проблеме различия между опциями '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.

0 голосов
/ 07 сентября 2018

Возможно, я неправильно понял, но собственные значения и собственные векторы (с точностью до плавающей запятой) идентичны: 

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
...