Symmetric eigen memory usage for eigenvectors.

[m-a-stewart]

Commit 57e9a0d to reduce allocations in syevr! seems to have a potential negative impact on the method

eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)

when the interval is such that a small number of eigenvalues/vectors are returned. The resizing and reshaping of the storage passed to LAPACK in syevr! means that space for an $n\times n$ eigenvector matrix is retained even if only a few eigenvectors are returned. This can be painful for large matrices.

With 1.12.0-rc1 the following code illustrates this:

using LinearAlgebra

n = 2000
A = Hermitian(randn(n,n))
e, V = eigen(A)
println("All eigenvectors:")
@show size(V)
@show Base.summarysize(V)

println("\nEigenvectors in [-1,1]:")
e1, V1 = eigen(A, -1, 1)
@show size(V1)
@show Base.summarysize(V1)

println("\nEigenvectors in [-1,1], copied:")
V1 = copy(V1)
@show size(V1)
@show Base.summarysize(V1)

The output is:

All eigenvectors:
size(V) = (2000, 2000)
Base.summarysize(V) = 32000048

Eigenvectors in [-1,1]:
size(V1) = (2000, 30)
Base.summarysize(V1) = 32000048

Eigenvectors in [-1,1], copied:
size(V1) = (2000, 30)
Base.summarysize(V1) = 480048

[dkarrasch]

IIRC then that commit doesn't change the behavior, it just avoids creating copies of arrays. Do you observe a difference in array size compared to, for instance, v1.11?

[m-a-stewart]

IIRC then that commit doesn't change the behavior, it just avoids creating copies of arrays. Do you observe a difference in array size compared to, for instance, v1.11?

Yes. A run of the same code with 1.11.6 gives:

All eigenvectors:
size(V) = (2000, 2000)
Base.summarysize(V) = 32000048

Eigenvectors in [-1,1]:
size(V1) = (2000, 29)
Base.summarysize(V1) = 464048

Eigenvectors in [-1,1], copied:
size(V1) = (2000, 29)
Base.summarysize(V1) = 464048

The copying that was previously in syevr!, while undesirable when computing $n$ eigenvectors, seems beneficial in this case. It allows garbage collection of the much larger array that needs to be passed to LAPACK to potentially hold $n$ eigenvectors.

[dkarrasch]

I'm not sure I'm following. This is a random matrix with a random spectrum, so it could be that there is one less eigenvalue in the unit disk. Or, maybe I need help: where exactly do I need to look, what behavior do you expect?

[m-a-stewart]

I'm not sure I'm following. This is a random matrix with a random spectrum, so it could be that there is one less eigenvalue in the unit disk. Or, maybe I need help: where exactly do I need to look, what behavior do you expect?

Sorry. I did think to mention something about the exact number of eigenvalues in the interval being not significant, but then didn't follow through on that. The random matrix was just easy to generate and the number of eigenvalues in the interval don't vary that much each time it is run.

What I was trying to illustrate with the example is that in 1.11, 30 eigenvectors were returned in a $2000\times 30$ matrix that used less than 0.5MB. In 1.12, about the same number of eigenvectors were returned in a $2000\times 29$ matrix that used about 31MB, which is the same memory that would be used by a $2000\times 2000$ matrix.

[dkarrasch]

Ah, now I spotted the difference. Thanks!

[dkarrasch]

How does that express itself in terms of performance? I guess the summarysize reflects the fact that the array was larger initially, but it was allocated in either case. I think we need a compelling use case to weigh the pros and cons of the "new" behavior.

[m-a-stewart]

My use case is making a project I'm working on have a smaller memory footprint when working with large matrices. I have some code for the symmetric definite generalized eigenvalue problem. The current code was proof of concept for a paper I wrote last year. I was pretty relaxed about allocation of potentially large matrices and I have a branch where I'm working on doing more things in place and making it more practical for larger matrices. Some of that is avoiding allocations where possible, but I'm also concerned about memory used by what my code returns. For the generalized eigenvalue problem, you start with two matrices. A $20,000\times 20,000$ matrix is certainly within reach of dense matrix algorithms. Such a matrix requires about 3GB. If you also use 3GB used to store a $20,000\times 30$ matrix, that has a pretty significant impact on overall memory use. Although I know you can't avoid allocating the $20,000 \times 20,000$ to set up the call to syevr!, the way things are now, on my larger tests, I get a 3GB $20,000\times 30$ matrix and the unused storage won't be garbage collected. That's a lot of unused memory to hold onto just to avoid an allocation.

I would also argue that this defies reasonable user expectations for the specific method for computing eigenvalues/eigenvectors in an interval. The main reason to use this method is to save on time and storage. Time savings are somewhat limited by the fact that a full tridiagonal reduction is required either way. The time savings are really in the computation of a limited number of eigenvectors. Just playing around with it with some non-rigorous benchmarks, if you don't need many eigenvectors, maybe it runs in 1/3 the time, which is nice, but not earth-shattering. However, when it comes to storage required for the eigenvector matrix it can potentially be orders of magnitude better. Returning a small $20,000\times 30$ matrix that still requires 3GB seems inconsistent with what a user expecting reduced memory usage might want from this method. Further, this behavior is not documented. Anyone using this method on large matrices is going to have to do some debugging to figure out why their code is holding onto so much memory. A 3GB $20,000\times 30$ matrix is not going be an obvious suspect.

Note that I'm not arguing more generally for other methods. I think the new changes are great for other methods that compute all eigenvectors. And having figured out where all my memory went, I have an obvious solution of just returning a copy of the eigenvector matrix so the small 3GB matrix can be garbage collected. So I'm not personally in a difficult spot with this. But I do think that an allocation for the final eigenvector matrix makes sense for this specific method.

[jaemolihm]

It appears that the issue is that resize! does not modify the size of the underlying MemoryRef. Thus, the resulting array may have a small size, but its ref is much larger.

julia> Z = zeros(100 * 100);

# Without sizehint! (current behavior)
julia> A = reshape(resize!(Z, 100 * 1), 100, 1);
julia> Base.summarysize(A)
80048
julia> length(A.ref.mem)
10000

# With sizehint!
julia> B = reshape(sizehint!(resize!(Z, 100 * 1), 100 * 1), 100, 1);
julia> Base.summarysize(B)
848
julia> length(B.ref.mem)
100

So the fix might be reshape(sizehint!(resize!(Z, ldz * zm), ldz * zm), ldz, zm). I have not tested it.

[m-a-stewart]

sizehint! looks like the behavior I would prefer. It would not introduce an allocation if all eigenvalues/eigenvectors are computed. I suppose there might be some objection that an allocation might not be helpful if you get a large number of eigenvalues in the interval, but I wouldn't have thought that would be the most important use case for a method that allows you to isolate just the eigenvalues in an interval. Certainly the possibility of holding onto enough memory for a full matrix of eigenvectors seems like something that will trip up people who really only want a small number of eigenvectors. It caught me by surprise to the extent that I ran out of memory while thinking I was saving only a small number of vectors.