Methods to diagonalize a Matrix{BigFloat}?

Specific Domains Numerics
1

Does anyone know about methods (or packages) that allow to diagonalize a Matrix{BigFloat}?

So far, I have only found (in IterativeSolvers.jl) that powm and invpowm do work with Matrix{BigFloat}, but using them I get only the largest and smallest eigenvectors, and I need all of them; svdl and eiglancz throw MethodErrors.

4 Likes
2

Would GenericSVD.jl or LinearAlgebra.jl do the trick?

3 Likes
3

Thanks @ChrisRackauckas, both seem to do what I need.

4

Which methods of LinearAlgebra and GenericSVD did you use to get all eigenvalues in BigFloat precision?

5

Don’t know which ones he ended up using, but LinearAlgebra.jl has a generic eigvals. Stuff isn’t documented or exported, but if you poke around you get to it:

julia> using LinearAlgebra

julia> A = big.(rand(2,2))
2×2 Array{BigFloat,2}:
 9.425471326934482529935621641925536096096038818359375000000000000000000000000000e-01  …  7.802185845048341672480773922870866954326629638671875000000000000000000000000000e-01
 5.767098416694038665042398861260153353214263916015625000000000000000000000000000e-01     1.806536534479177280587691711843945086002349853515625000000000000000000000000000e-01

julia> LinearAlgebra.EigenGeneral.eigvals!(A)
2-element Array{Complex{BigFloat},1}:
      1.333014778966178347989456524738500823065332531245815992476150525785817496429229+0.000000000000000000000000000000000000000000000000000000000000000000000000000000im
 -2.098139928248123669371251893615527048554936640583159924761505257858174964292295e-01-0.000000000000000000000000000000000000000000000000000000000000000000000000000000im
1 Like
6

Both packages work as expected. If you only need the eigenvalues, I think LinearAlgebra.jl is what you need; if you are also interested in the eigenvectors, GenericSVD.jl provides them.

1 Like
7

There does not seem to be a LinearAlgebra.EigenGeneral.eigvals! any more in Julia 1.0.
Is anyone aware of any other package implementing generic linear algebra routines? (besides GenericSVD)
(Sorry for the necrobump, but I figured out that the package ecosystem might have evolved since this question was first asked)

1 Like
8

What used to be LinAlg became the the stdlib LinearAlgebra so I renamed my LinearAlgebra to GenericLinearAlgebra. You no longer have to qualify the call to the BigFloat solver but can just do

using GenericLinearAlgebra

eigvals(big.(randn(4,4)))
1 Like
9

Thanks!

Do you know if there is a way to compute the eigenvectors as well? (either in this package or somewhere else)

What I am trying to do is simply to compute the square root of a matrix of Complex{BigFloat}, but the default method from Base fails because of the missing Schur decomposition. So I was trying to diagonalize the matrix as a fallback.

The reason behind all this is that I very recently translated some Fortran routine for the Takagi factorization to Julia (package here), and although it now “works” with arbitrary floating point types, it seems to have numerical stability issues when the input matrix contains elements varying by many orders of magnitude (even when using BigFloats with 1000’s of bits).

So I wanted to compare it to another algebraic method (link to paper) which relies on the SVD and is supposedly exact. I used the SVD from GenericSVD.jl, but I still need to find a way to compute the square root of an arbitrary complex matrix.

As a last resort, I can check in some linear algebra book if there is a way to get the eigendecomposition (or directly the square root) from the SVD of a hermitian matrix.

(For reference, it seems that eigvals is not exported, so I had to use GenericLinearAlgebra.eigvals instead.)

10

Wikipedia answered my question :slight_smile:

11

If your matrix is Hermitian, you can do

eigen(Hermitian(big.(randn(4,4))))

and get the vectors as well. I haven’t put all the pieces together for the non-Hermitian solver yet. I think all the individual components necessary for computer the eigenvectors are there but they’ll have to be combined and tested.

2 Likes
12

Thanks! Your example seems to work. I’ll have a quick look at the code as well.

I will compare the results to the SVD-based method, so this might help a bit with the testing.

13

FYI, I did a few tests with real symmetric matrices, and your implementation matches the SVD-based method with very high accuracy.

Both methods are capable of diagonalizing outrageously hierarchical matrices. I tried diagonalizing a 5×5 matrix with eigenvalues spanning 35 orders of magnitude, and both methods were accurate to within only 10^{35}\epsilon (using 1024 bit BigFloats, so \epsilon \sim 10^{-308}). The SVD-based method was slightly (10^3 times) more accurate, but in these extreme situations it does not really matter.

However, I am having issues diagonalizing matrices of Complex{BigFloat} using your package. The following code:

julia> using LinearAlgebra, GenericLinearAlgebra
julia> m = Hermitian(rand(Complex{BigFloat}, 5, 5));
julia> eigen(m)

leads to a method ambiguity with LinearAlgebra:

ERROR: MethodError: LinearAlgebra.eigen(::Hermitian{Complex{BigFloat},Array{Complex{BigFloat},2}}) is ambiguous. Candidates:
  eigen(A::Hermitian) in GenericLinearAlgebra at /home/jl/.julia/packages/GenericLinearAlgebra/s9cSL/src/eigenSelfAdjoint.jl:630
  eigen(A::Union{Hermitian{T,S}, Hermitian{Complex{T},S}, Symmetric{T,S}} where S where T<:Real) in LinearAlgebra at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/symmetric.jl:490
Possible fix, define
  eigen(::Union{Hermitian{T,S} where S<:(AbstractArray{#s549,2} where #s549<:T) where T<:Real, Hermitian{Complex{T},S} where S<:(AbstractArray{#s549,2} where #s549<:Complex{T}) where T<:Real})
Stacktrace:
 [1] top-level scope at none:0

Am I doing something wrong or is this a bug?

14

It is a bug. I thought I had fixed but apparently I didn’t. I’ll take a look.

15

Thanks!

16

Thinking about it again, the matrix which I need to diagonalize is not Hermitian, but unitary (so it is still normal and should have an eigendecomposition, but with complex eigenvalues).

Wikipedia claims that there is nonetheless a way to get the SVD from the eigendecomposition in this case. I will see if I can figure it out how to do the reverse operation, i.e. obtain the eigendecomposition from the SVD for an arbitrary normal matrix.

17

Did you work this out? For a unitary matrix, the Schur decomposition (which is mostly implemented in @andreasnoack’s package) should be a diagonalization. I wanted a full generic Schur at some point, but got lost trying to combine the Householder and Givens sequences in Andreas’ code into a unitary factor (what you would need for eigenvectors), and put it off.

18

For the problem at hand, I actually only need the singular values. I wanted to compute the full Takagi factorization for validation purpose, but I ended up checking the formula using Mathematica.

I may have a look at it at a later time, since I would rather use the exact formula in the TakagiFactorization package, and this requires me to diagonalize a unitary matrix. But I have not had the time so far, and my knowledge of numerical linear algebra is quite limited.

19

What is the status of the problem? Has anything changed after a year?

I am asking because I need to find both eigenvalues and eigenvectors of general complex matrices in arbitrary precision. It seems that GeneralLinearAlgebra has some routines for eigenvalues and Schur decomposition. May be it is possible to get eigenvectors somehow from Schur decomposition?

20

The GenericSchur package provides eigenvectors etc. for arbitrary precision complex matrices.

4 Likes