Sparse matrix `eigsolve` defaults to non-sparse, dense `x₀`

[BenjaminRemez]

https://github.com/Jutho/KrylovKit.jl/blob/68ced469e833613e00c70871f82846b525a9d7c8/src/eigsolve/eigsolve.jl#L189

For sparse A, above line leads to a dense SparseVector x₀. Being dense, multiplications A*x₀ are slow and lead to excessive allocations. In my benchmarking, runtime for A with density ~5% exceeds that of calling eigsolve(collect(A)...). This also leads to inconsistent runtime/memory usage between calling eigsolve(A,...) and eigsolve(A, size(A,1), ...), as the latter correctly initializes a x₀::Vector{Float64}.

[lkdvos]

Hi Benjamin,

Thanks for the report, and also for the PR, this is definitely appreciated. I think it is an interesting problem, but there are several things to consider here:

• redefining the initialization to always allocate rand(T, size...) effectively means we always use Vector for AbstractMatrix. This is a rather opinionated choice, for example for CuMatrix this would immediately get us into trouble.
• Ultimately, it is very hard to write a solution which will work for every single case that someone might come up with.
• You always have the option to provide your own initial vector

In some sense, I think a cleaner solution might be to make this an explicit part of the interface, ie have some function that may be overloaded, something like:

initialize_vector(A, T) = rand!(similar(A, T))

which allows for more customization. Then, we might consider a package extension to implement this?

[BenjaminRemez]

I agree that right now there is a simple workaround for this problem. However, I think my focus is on the fact that this behavior went against my expectation of "all else being equal, calling eigsolve on A::SparseMatrixCSC should be faster than on A::Matrix, if A is large and very sparse". That right now you get the expected speedup (in my use case now, a respectable 30x ) for eigsolve(A,n,...), but instead a 10x slowdown for eigsolve(A,...) is the problem, which is undocumented and might be a pitfall for other new users.

I don't have particularly strong opinions on what the solution ought to be, so your suggestion certainly sounds like a good way to go :-) Here you mean that KrylovKit would have extensions for SparseArrays, GPUArrays, etc., or the other way around?

[BenjaminRemez]

Actually, it is not "undocumented", it contradicts the documentation, which would behave like my suggested commit: https://github.com/Jutho/KrylovKit.jl/blob/68ced469e833613e00c70871f82846b525a9d7c8/src/eigsolve/eigsolve.jl#L15