v0.5 issue

[matthieugomez]

trafficstars

With v0.5, there is a method error if the initial vector is Abstract, e.g.

using KrylovKit
n = 10
KrylovKit.eigsolve(rand(n, n), 1:n, 1, :LM)
# ERROR: MethodError: no method matching KrylovKit.ArnoldiFactorization(::Int64, ::KrylovKit.OrthonormalBasis{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Array{Float64,1}, ::Array{Float64,1})
 KrylovKit.eigsolve(rand(n, n), collect(1:n), 1, :LM)
# works

[rveltz]

I am not sure ranges follow the requirements of the package (axpy,...): https://jutho.github.io/KrylovKit.jl/latest/#Package-features-and-alternatives

[matthieugomez]

I believe it should work for any AbstractVector.

[Jutho]

As indicated by @rveltz , ranges do not support the minimal required interface. All AbstractVectors will probably work, provided they are at least mutable. As KrylovKit tries to assume as little as possible about the type of vectors that the user works with, I think that for this case the fix is much easier on the user side, by simply collecting the range into a vector. At the package side, I could at best provide a better error message.

[matthieugomez]

As KrylovKit tries to assume as little as possible about the type of vectors that the user works with,

No but, come on, in this instance, KrylovKit assumes that the second argument is a Vector. That is, the same error message would appear for a custom AbstractVector type that would satisfy the exact same interface as a Vector (incl. setindex!, axpy etc). So this is a bug.

This happens because, in https://github.com/Jutho/KrylovKit.jl/blob/master/src/factorizations/arnoldi.jl#L136, v is an AbstractVector while r is a Vector. ArnoldiFactorization expects them to have the same types

It could be fixed by replacing

    v = (one(T) / β₀) * x₀

by the commented version

v =  mul!(similar(x₀, T), x₀, 1/β₀)

[Jutho]

The problem is that I no longer depend on similar(x₀, T) to create a similar vector with a different scalar type, because I found that interface somewhat contrived. The problem, in my opinion, being that Julia conflates with eltype the type of what you get upon indexing/iteration, versus the scalar type over which this forms a vector space. Nested vectors would be an example where these two are different things.

Hence, multiplying with a scalar of a particular type is now the way I create similar vectors with a possibly different scalar type. It is somewhat unfortunate that a range decides to remain a range after being multiplied by a floating point number. At best, I could do something like

mul!(one(T)*similar(x₀), x₀, 1/β₀)

but then I will be allocating a similar(x₀), the memory of which is immediately lost by multiplying it out of place with another number. Maybe one additional allocation is not too bad, but it is not nice either.

[matthieugomez]

Yes it’s not obvious. How about relaxing the ArnoldiFactorization type so that v and r can have different types? Would that solve the issue?

[Jutho]

No I don't think this helps. After the first few lines in the initialize routine, v should have the type of all Krylov vectors that will be generated, r is the residual. The residual becomes the next Krylov vector (after orthonormalizing it with respect to the current Krylov subspace). So they should have the same type. Furthermore, similar(v) should create vectors of the same type of v. This behaviour not respected by ranges. I really think that ranges just don't satisfy the required minimal interface, and that it is much easier to just collect them as you did. Is there actually anything important about the fact that it is a range that you want to preserve in the computation?

[Jutho]

Ok, I found a replacement that works and does not result in additional allocations. If it passes my tests, it will be the next commit.

[Jutho]

This is no longer an issue so I will close this.