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

Specific Domains Numerics
1

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.

1 Like
2

Please provide the example that didn’t work.

1 Like
3

I think you might be interested in LinearMaps.jl, and perhaps also Arpack alternatives such as ArnoldiMethod.jl or KrylovKit.jl .

4

Regarding LinearMaps.jl, there is an example here.

5

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:

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")