Unexpectedly slow performance of "eigs" with a Hermitian view

General Usage Performance
1

Solving quantum mechanics problems numerically (via the so-called “exact diagonalization” method) requires finding the lowest several eigenvalue and eigenvector of a large sparse Hermitian matrix. Since one wants to do this for matrix of size 2^L for L as large as possible, I want to understand the optimal way to do this in Julia.

The basic command for this is the eigs function from the Arpack package. (Which earlier posts such as https://discourse.julialang.org/t/special-method-for-eigenvectors-of-sparse-hermitian-complex-matrix/6233 point out does not have a special method for Hermitian matrices and defaults to the fully complex case.) Since we know the matrix is Hermitian, it makes sense to use a Hermitian view to store it, which takes half the memory since we only need the upper triangle.

Surprisingly, when using the Hermitian view for the matrix it takes several orders of magnitude (!) longer for the Hermitian view of the matrix than the full matrix. Is this due to a bug, some subtle way that views are implemented that makes this especially slow, or something else?

Example:

using LinearAlgebra
using SparseArrays
using Arpack
using BenchmarkTools

const N = 500 #size of the matrix
const fill = 100 #number of eigenvalues to fill

#make a random matrix
rand_cols = rand(1:N,fill);
rand_row = [rand(rc:N) for rc in rand_cols];
rand_values = rand(Complex{Float64},fill);

upper_tri = sparse(rand_cols,rand_row,rand_values,N,N)
upper_tri = upper_tri+adjoint(Diagonal(upper_tri)) #make sure the diagonal is real

#view of the matrix using the only upper half
hermitian_view = Hermitian(upper_tri,:U) 
#flip across the diagonal
hermitian_full = upper_tri + adjoint(upper_tri) - Diagonal(upper_tri) 

@show ishermitian(hermitian_full) #true
@show ishermitian(hermitian_view) #true
@show isapprox(Matrix(hermitian_full),Matrix(hermitian_view)) #they're the same matrix

evs1=eigs(hermitian_full; which=:SR, nev=2)[1]
evs2=eigs(hermitian_view; which=:SR, nev=2)[1]
@show isapprox(evs1,evs2) #same matrix => same eigenvalues

@show @btime eigs($hermitian_full; which=:SR, nev=2)[1]
@show @btime eigs($hermitian_view; which=:SR, nev=2)[1]

The last two lines gives
hermitian full → 3.853 ms (482 allocations: 243.83 KiB)
hermitian_view → 446.538 ms (217346 allocations: 10.17 MiB).

This is unchanged under any conditions I can think of, including wrapping everything in functions. This has been present since at least Julia 0.4, and persists under 1.0.

Thank you.

2

I think what you are seeing is Slow sparse matrix-vector product with symmetric matrices

julia> @btime mul!($c, $A, $b);
  120.341 μs (0 allocations: 0 bytes)

julia> _A = Hermitian(A);

julia> @btime mul!($c, $_A, $b);
  37.587 ms (4008 allocations: 187.88 KiB)
1 Like
3

Thanks for the link and quick reply. Let me see if I understand this correctly: this is a (somewhat) known issue that matrix multiplication with Symmetric or Hermitian views is very slow because 1. the views wrap full matrices instead of sparse ones and 2. there’s no specific implementation for matrix-vector multiplication which takes advantage of the extra structure. Is that right?

Unfortunately, the thread indicates there is no bug report or issue related to this. Do you know if this is still the case? It seems to me that if these views exist that indicate your matrix has extra structure, then the most basic operations of matrix-matrix and matrix-vector multiplication should take advantage of them.

4

That’s my understanding from the thread.

I don’t know if an issue was opened or not, but certainly the problem was not solved. But then this is open source, so the first person who really needs this feature is likely to be the one making the PR. Memory was not a bottleneck for me so I ended up using the underlying whole matrix. @dpo wrote some implementations in bmark_v1.jl · GitHub from that thread so he may have something to say.

5

As far as I know, there is still no resolution to the symmetric/hermitian matrix-vector product. I opened a pull request over a year ago but I’m told the ramifications of such changes are intricate: https://github.com/JuliaLang/julia/pull/22200. In the mean time, you should be able to use one of the implementations from my gist or from the pull request.

1 Like
6

I also have opened a somewhat related issue:

https://github.com/JuliaLang/julia/issues/28432

My take on this and similar issues is as follows: the standard library currently defines many matrix types that are sparse matrices (in the mathematical sense) but do not belong to an abstract sparse type. In particular, there is no AbstractSparseMatrix type in Julia. So the compiler falls back to AbstractArray (i.e., dense) matrix operations when it sees a function that hasn’t been implemented for the specific matrix type.

EDIT: Sorry, I misspoke: there is an AbstractSparseMatrix type in the standard library, but the necessary subtype structure and associated methods aren’t available.

EDIT: By the way, a work-around for related problems (not sure about yours) is to invoke the SparseMatrixCSC converter on your sparse-matrix object: use SparseMatrixCSC(A) in place of A in the argument list for the eigs call.

7

Thank you everyone for the careful responses. It’s somewhat frustrating that progress on this issue seems to have stalled, though I certainly recognize that its relatively low priority since there are clear and straightforward workarounds. Perhaps I’ll try to dig into what’s involved and see if I can make any headway.

The suggestion from @Stephen_Vavasis of wrapping everything in a SparseMatrixCSC call in eigs does “solve” the problem, but has to store both the lower and upper triangles of the matrix. It seems the best thing to do for the time being is live with double the memory usage. Of course, I’ll explore some of the PR’s that people linked here and see if those can do a better job.