`eigs` for sparse tridiagonal matrix: Arpack.jl vs scipy


using LinearAlgebra, SparseArrays, Arpack, PyCall, BenchmarkTools
@pyimport scipy.sparse as pysparse
@pyimport scipy.sparse.linalg as pylinalg

N = 100_000
X = sparse(Tridiagonal(1.0:N-1, 1.0:N, 1.0:N-1))
X = X + X'
Xpy = pysparse.csc_matrix((X.nzval, X.rowval .- 1, X.colptr .- 1), shape=size(X))

@btime Arpack.eigs($X, nev=6, which=:LM); # julia
@btime pylinalg.eigs($Xpy, k=6, which=:LM); # python

# Output:
#  3.475 s (4467 allocations: 23.83 MiB)
#  3.116 s (180 allocations: 9.16 MiB)

I’m setting up a sparse (about 1e-5 nz entries) tridiagonal real matrix. I want to know the, say, 6 largest eigenvalues.

Why is Julia slower here?

AFAIU, both scipy.sparse.linalg.eigs and Arpack.jl’s eigs are using ARPACK. Correct?

Any explanation is highly appreciated!

One more thing: what about scipy.sparse.linalg.eigsh? I thought this is a ARPACK function as well but can’t find a wrapper in Arpack.jl.

What’s the performance like if you just use a regular sparse matrix instead of tridiagonal?

Basically the same. In fact, I tried it with random sparse matrices first.