Generalised eigenvalue problem with BigFloat

Specific Domains Numerics
,
1

I am trying to solve the generalised eigenvalue problem in arbitrary precision, using the BigFloat data type.

using LinearAlgebra
using GenericLinearAlgebra

setprecision(128)

println(eigen(Hermitian(rand(BigFloat, 2, 2)), Hermitian(rand(BigFloat, 2, 2))))

It produces this error:

Exception has occurred: MethodError

MethodError: no method matching eigen!(::Hermitian{BigFloat, Matrix{BigFloat}}, ::Hermitian{BigFloat, Matrix{BigFloat}}) 
The function `eigen!` exists, but no method is defined for this combination of argument types.

However, using BigFloat for the standard eigenvalue problem works

println(eigen(Hermitian(rand(BigFloat, 2, 2))))

Does anyone know why this is or how I could fix it? Thanks so much for any help, I’ve been asked to use Julia for my masters project and have no experience using it!

1 Like
2

The symmetric eigenvalue problem is easier to implement

Arb.jl has a general eigenvalue solver with its own arbitrary precision ball arithmetic, getting rigorous bounds too. This might suit your needs

1 Like
3

I don’t think such a package exists. I guess at least one of Arblib.jl or ArbNumerics.jl should suffice, though.

4

There are methods for generalized eigensystems of complex matrices with arbitrary precision in GenericSchur.jl. I haven’t gotten around to writing all the conventional wrappers but you can do something like

using LinearAlgebra, GenericSchur
A = rand(Complex{BigFloat},n,n)
B = rand(Complex{BigFloat},n,n)
S = schur(A,B)
V = eigvecs(S)
# check result:
resid = norm(A * V - V * B * Diagonal(S.values))

I haven’t implemented methods for the plain real or Hermitian cases.

2 Likes
5

In the Hermitian/real symmetric case, you usually transform the problem by Cholesky factorizing the B matrix. I.e. something like

julia> using LinearAlgebra

julia> using GenericLinearAlgebra

julia> setprecision(128)
128

julia> A = Hermitian(rand(BigFloat, 2, 2))
2×2 Hermitian{BigFloat, Matrix{BigFloat}}:
 0.750726  0.800863
 0.800863  0.146439

julia> B = Hermitian(rand(BigFloat, 2, 2))
2×2 Hermitian{BigFloat, Matrix{BigFloat}}:
 0.497388   0.0396879
 0.0396879  0.298677

julia> Fchol = cholesky(B)
Cholesky{BigFloat, Matrix{BigFloat}}
U factor:
2×2 UpperTriangular{BigFloat, Matrix{BigFloat}}:
 0.705257  0.0562743
  ⋅        0.543609

julia> Ã = Fchol.U'\A/Fchol.U
2×2 Matrix{BigFloat}:
 1.50934  1.93268
 1.93268  0.0792299

julia> Feigen = eigen(Ã)
Eigen{BigFloat, BigFloat, Matrix{BigFloat}, Vector{BigFloat}}
values:
2-element Vector{BigFloat}:
 -1.266437199907713533384656554486557838751
  2.855004257222210224062038053659169376169
vectors:
2×2 Matrix{BigFloat}:
  0.571405  -0.820668
 -0.820668  -0.571405

julia> V = Fchol.U\Feigen.vectors
2×2 Matrix{BigFloat}:
  0.930669  -1.07977
 -1.50967   -1.05113

julia> norm(A*V - B*V*Diagonal(Feigen.values))
4.156000133564589919546802616566101790705e-39
2 Likes
6

Thankyou! Are you aware of any additional time constraints the cholsesky factorisation would have on the calculation? I am using matrix sizes of > (30)^3. If not don’t worry I can do some analysis :slight_smile:

7

Thanks so much !! :))

8

The Cholesky factorization wouldn’t be the limiting factor here. The eigenvalue computation would be much more costly than the Cholesky. Using a symmetric eigen solver, as would be available when transforming the problem with the Cholesky factor, is also typically much faster than using a general solver. However, if the dimension is ~30^3 then I think it will be almost prohibitively costly to solve with BigFloat elements. Do you need all the eigenvalues? If not, then you should probably consider an iterative solver.

1 Like
9

I am aware of the computational cost required, this code has already been used in python to solve the non-BO Schrodinger equation (w/ Float64). Implementation into Julia is in order to increase precision, as it has been extraordinarily difficult in python. The python code also uses sparse matrices, but as far as I am aware, there is no package that allows for solving the eigenvalue problem with precision > 64 bits for sparse matrices in Julia. For this purpose, only the minimum eigenvalue(s) is required. Do you know an iterative solver that would work using the BigFloat type?

10

Do you know an iterative solver that would work using the BigFloat type?

I believe that should be possible with either Home · ArnoldiMethod.jl or Home · KrylovKit.jl but I haven’t tested it.