# Eigenvalues and extended precision

**URL:** <https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119>\
**Category:** Numerics\
**Tags:** linearalgebra, precision, eigenvalues\
**Created:** [September 13, 2021, 7:48pm UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119 "2021-09-13T19:48:32Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)\
**Post date:** [September 13, 2021, 7:48pm UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/1 "2021-09-13T19:48:32Z")

</div>

I am interested in computing a small proportion (between 3 and 10% of the matrix size) of the eigenvalues of a matrix with extended precision. As a toy example, let’s take the matrix `A` define by:

```julia
T = Complex{Float128}
N = 64
A = Tridiagonal(rand(T, N - 1), rand(T, N), rand(T, N - 1))

```

The following eigensolver does not works:

```julia
julia> using Arpack
julia> eigs(A; nev=round(Int, 0.1 * N), which=:LM, ritzvec=false)
ERROR: LoadError: StackOverflowError:
...

```

from the `Arpack.jl` package;

```julia
julia> using KrylovKit
julia> eigsolve(A, round(Int, 0.1 * N), :LM)
ERROR: LoadError: MethodError: no method matching hschur!(...

```

from the `KrylovKit.jl` package;

```julia
julia> using ArnoldiMethod
julia> partialschur(A; nev=round(Int, 0.1 * N), which=LM())
ERROR: LoadError: MethodError: no method matching gemv!(...

```

from the `ArnoldiMethod.jl` package.

Is there a way/package to compute a small proportion of the eigenvalues with extended precision ?

* * *

Additional information:

- I have used `Float128` from the package `Quadmath.jl` but I could have use `Double64` from the `DoubleFloats.jl` or `Float64x2` from `MultiFloats.jl`, I get the same errors.
- `eigvals(Matrix(A))` from the standard library `LinearAlgebra.jl` works with extended precision but here the size `N` is small for exposition purpose, in my application it might be as big as 10 000 and I do not need all the eigenvalues only a small proportion of them. So if I can avoid the huge cost of computing all the eigenvalues, that would be great.
- All the example shown above works with `T = Complex{Float64}`.
- I know about the eigensolver `powm` in the `IterativeSolvers.jl` but it only give one eigenvalue and I am interested in a small proportion \> 1.

Edit: Formatting and add details.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [September 13, 2021, 8:32pm UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/2 "2021-09-13T20:32:47Z")

</div>

You can get the job done with GenericLinearAlgebra.jl and ArnoldiMethod.jl:

```julia
julia> partialschur(A; nev=round(Int, 0.1 * N), which=LM())[1]
PartialSchur decomposition (Complex{Float128}) of dimension 6
eigenvalues:
6-element Vector{Complex{Float128}}:
 1.76732917570178181966857786046259500e+00 + 1.88566438666511933882975520031760772e+00im
 1.65431839639080928590829614759631355e+00 + 1.97529966730726403527537198010151792e+00im
 1.29189845962540662486402509909018053e+00 + 1.85730981985326261230642237981537783e+00im
 1.29393589581771217265696071928253209e+00 + 1.85391512578430311844047859807080534e+00im
 1.03405214172759353767334445901983764e+00 + 1.94335216348051256609276579410600025e+00im
 1.66370260978943467402009678221541037e+00 + 1.42582140642411450137051978349554524e+00im

```

It’s dog-slow, though: `N = 1024` took ~3 minutes, versus 0.3 seconds for `Complex{Float64}`.

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [September 13, 2021, 8:36pm UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/3 "2021-09-13T20:36:34Z")

</div>

If you convert your Tridiagonal matrix to a general Matrix your code should work using `GenericLinearAlgebra.eigvals`:

``  
using LinearAlgebra, GenericLinearAlgebra, Quadmath, DoubleFloats

complex\_matrix(T,n) = rand(Complex{T}, n) \* rand(Complex{T}, n)’  
generic\_tridiag(m) = Matrix(Tridiagonal(m))

n = 4  
m128 = generic\_tridiag(complex\_matrix(Float128, n));  
md64 = Complex{Double64}.(m128);

m128eigvals = eigvals(m128)  
md64eigvals = eigvals(md64)

maximum(abs.(m128eigvals .- md64eigvals))

```julia

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 14, 2021, 6:35am UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/4 "2021-09-14T06:35:06Z")

</div>

Have you tried KrylovKit?

---

<div class="post-metadata">

**Author:** ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)\
**Post date:** [September 14, 2021, 6:51am UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/5 "2021-09-14T06:51:29Z")

</div>

Can you elaborate on how you use the `GenericLinearAlgebra.jl` package from their GitHub I do not understand how to use it in my case.

---

<div class="post-metadata">

**Author:** ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)\
**Post date:** [September 14, 2021, 7:20am UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/6 "2021-09-14T07:20:29Z")

</div>

Yes, I have already see that `eigvals` from the `LinearAlgebra.jl` works with extended precision (and you do not need the `GenericLinearAlgebra.jl` package for `eigenvals` to works) but the two mains issue with this method:

1. It transform the tridiagonal matrix into a dense matrix which become prohibitive as the matrix size grows.
2. It compute all the eigenvalues but I am interested in a small proportion of the eigenvalues (between 3 and 10% of the matrix size). So if I could avoid the cost of computing all the eigenvalues that would be great.

---

<div class="post-metadata">

**Author:** ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)\
**Post date:** [September 14, 2021, 7:35am UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/7 "2021-09-14T07:35:40Z")

</div>

Yes, I tried `KrylovKit.jl` in my second example (I have edited my post to be hopefully more clear).

It is possible that I am not using it properly for what I want to do. I have seen that you can give it the product matrix-vector instead of the matrix but I did not manage to make it works with extended precision. I get the same error.

```julia
julia> using KrylovKit
julia> Amap = LinearMap{T}(x -> A * x, N)
julia> eigsolve(Amap, N, nb_eig, :LM)[1]
ERROR: LoadError: MethodError: no method matching hschur!(...

```

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [September 14, 2021, 8:39am UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/8 "2021-09-14T08:39:09Z")

</div>

I think what you need in addition is `using GenericSchur`. At least this is what allows one to use `ArnoldiMethod.jl` with arbitrary float types.

```julia
julia> using ArnoldiMethod, GenericSchur, DoubleFloats

julia> T = Complex{Double64}
Complex{Double64}

julia> N = 64
64

julia> A = Tridiagonal(rand(T, N - 1), rand(T, N), rand(T, N - 1));

julia> pschur,hist = partialschur(A; nev=round(Int, 0.1 * N), which=LM());

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 14, 2021, 8:53am UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/9 "2021-09-14T08:53:23Z")

</div>

You should open an issue because it should be type agnostic

---

<div class="post-metadata">

**Author:** ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)\
**Post date:** [September 14, 2021, 1:24pm UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/10 "2021-09-14T13:24:22Z")

</div>

Thank you @fgerick for the answers.

I have found the main issue which was that the `ArnoldiMethod.jl` package was pinned at version `0.1.0` because of `BifurcationKit.jl`. When I remove `BifurcationKit.jl` and update `ArnoldiMethod.jl` to `0.2.0` it worked.

Actually, it even worked without `GenericLinearAlgebra.jl` and `GenericSchur.jl`. The following works:

```julia
using LinearAlgebra, ArnoldiMethod, Quadmath

T = Complex{Float128}
N = 64

A = Tridiagonal(rand(T, N - 1), rand(T, N), rand(T, N - 1))

pschur, hist = partialschur(A; nev=round(Int, 0.1 * N), which=LM())
display(pschur)

```

---

<div class="post-metadata">

**Author:** ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)\
**Post date:** [September 14, 2021, 1:31pm UTC](https://discourse.julialang.org/t/eigenvalues-and-extended-precision/68119/11 "2021-09-14T13:31:15Z")

</div>

Ok thanks, I will

edit:  
[https://github.com/Jutho/KrylovKit.jl/issues/51](https://github.com/Jutho/KrylovKit.jl/issues/51)