# How to do Arpack.eigs(v-\>A\*v, ...)?

**URL:** <https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192>\
**Category:** Numerics\
**Created:** [January 28, 2019, 6:05pm UTC](https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192 "2019-01-28T18:05:25Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![bsxfan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bsxfan/32/6914_2.png) [@bsxfan](https://discourse.julialang.org/u/bsxfan)\
**Post date:** [January 28, 2019, 6:05pm UTC](https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192/1 "2019-01-28T18:05:25Z")

</div>

I need to find a few of the principal eigenvalues and eigenvectors of a very large (symmetric, real) matrix that is implicitly defined via its matrix-vector product. I know Arpack can do this, because it only needs to do a few matrix-vector products, but I am not managing to use the Julia interface (Arpack.jl) to do this. I found some discussion in various Julia forums and I tried my own solution where I declared a new AbstractArray type, which implements mul!() and a few other methods, like size() and issymmetric(), but so far I have been unsuccesful in making it work.

I tried this in Julia 1.1 on Windows. Any help or advice would be appreciated.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [January 28, 2019, 8:04pm UTC](https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192/2 "2019-01-28T20:04:36Z")

</div>

Please provide the example that didn’t work.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [January 28, 2019, 8:26pm UTC](https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192/3 "2019-01-28T20:26:47Z")

</div>

I think you might be interested in [LinearMaps.jl](https://github.com/Jutho/LinearMaps.jl), and perhaps also Arpack alternatives such as [ArnoldiMethod.jl](https://github.com/haampie/ArnoldiMethod.jl) or [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl) .

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 28, 2019, 8:42pm UTC](https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192/4 "2019-01-28T20:42:27Z")

</div>

Regarding LinearMaps.jl, there is an example [here](https://github.com/Jutho/LinearMaps.jl#example).

---

<div class="post-metadata">

**Author:** ![bsxfan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bsxfan/32/6914_2.png) [@bsxfan](https://discourse.julialang.org/u/bsxfan)\
**Post date:** [January 29, 2019, 1:05pm UTC](https://discourse.julialang.org/t/how-to-do-arpack-eigs-v-a-v/20192/5 "2019-01-29T13:05:47Z")

</div>

Thank you for the pointer to LinearMaps. I was able to get

`Arpack.eigs(f::LinearMaps.FunctionMap,...)`

working to my satisfaction. Below is a code example comparing implicit vs explicit matrices passed to `eigs()`. Comparing the two cases, the results agree and the timings (after warm-up) are very similar.

Notes:

- To signal to `eigs()` whether the matrix is symmetric, I use the `Symmetric` type for the explicit case; and set a flag in the `FunctionMap` constructor for the implicit case.
- In the explicit case, when passing a general, Array{Float64,2} matrix to eigs(), I think it somehow examines the matrix for symmetry, because the type of the returned eigenvalues is real if the matrix is symmetric and complex if assymmetric. In the implicit case, it does not do this and returns real/complex results depending on the flag `issymmetric`. According to how I understand type stability in Julia, I would say the latter (implicit case) behaviour is better. In the explicit case, the return type depends on the value (not type) of the input argument. Is it worth raising an issue in `Arpack.jl` about this?

I tested the code below in Julia 1.1 on windows 10:

```julia
using Arpack
using LinearMaps: FunctionMap
using LinearAlgebra: checksquare, Symmetric

# construct a real, symmetric (positive definite) example matrix
posdef = false
if posdef
    A = randn(500,600)
    A = A*A'
else
    A = randn(500,500)
    A = A + A'
end
println("A is: ",typeof(A))
S = Symmetric(A)

# and also an assymetric matrix (with complex eigenvalues)
B = randn(500,500)

#warm up
d,V,nconv,niter,nmult,resid = eigs(S,nev=3)
d,V,nconv,niter,nmult,resid = eigs(A,nev=3)
d,V,nconv,niter,nmult,resid = eigs(B,nev=3)

println(" **Symmetric example matrix**")
println("explicit eigs(), using symmetry:")
@time d,V,nconv,niter,nmult,resid = eigs(S,nev=3)
println("principal eigenvalues: ",d)
println("-------")

println("explicit eigs(), ignoring symmetry:")
@time d,V,nconv,niter,nmult,resid = eigs(A,nev=3)
println("principal eigenvalues: ",d)
println("-------")

# wrap this code in a function to help compiler optimize (very slow otherwise)
function imp_eigs(f,m,sym::Bool,n_ev::Int)
    #A_imp = FunctionMap{Float64}(f,nothing, m, ismutating=false, issymmetric=sym)
    A_imp = FunctionMap{Float64}(f,nothing, m, issymmetric=sym)
    d,V,nconv,niter,nmult,resid = eigs(A_imp,nev=n_ev)
    return d
end

#functions for implicit matrix representation
f(v) = A*v
g(v) = B*v
m1 = checksquare(A)
m2 = checksquare(B)

#warm up
d = imp_eigs(f,m1,true,3)
d = imp_eigs(f,m2,false,3)
d = imp_eigs(g,m2,false,3)

println("implicit eigs(), using symmetry:")
@time d = imp_eigs(f,m1,true,3)
println("principal eigenvalues: ",d)
println("-------")

println("implicit eigs(), ignoring symmetry:")
@time d = imp_eigs(f,m1,false,3)
println("principal eigenvalues: ",d)
println("-------\n\n")

println(" **Asymmetric example matrix**")
println("explicit eigs():")
@time d,V,nconv,niter,nmult,resid = eigs(B,nev=3)
println("principal eigenvalues: ",d)
println("-------")

println("implicit eigs():")
@time d = imp_eigs(g,m2,false,3)
println("principal eigenvalues: ",d)
println("-------\n\n")

```