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

[BenjaminRemez]

KrylovKit.jl/src/eigsolve/eigsolve.jl

Line 189 in 68ced46

x₀ = Random.rand!(similar(A, T, size(A, 1)))

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:

KrylovKit.jl/src/eigsolve/eigsolve.jl

Line 15 in 68ced46

provided, it is then chosen as `rand(T, size(A, 1))`. If the linear map is encoded more