Jutho / KrylovKit.jl

Krylov methods for linear problems, eigenvalues, singular values and matrix functions
Other
284 stars 37 forks source link

re-export VectorInterface? #88

Closed rveltz closed 4 months ago

rveltz commented 4 months ago

Hi,

I rely on KK for BifurcationKit.jl. In particular, I implemented a type BorderedArray

It seems the only method I need to add is add. First, I thought having axpby! defined would not require to implement add since I have similar and copy. Can you confirm this ?

In the case I really need to add the method add, I need to import VectorInferface: add which requires to add VectorInterface as a dependency to BifurcationKit. Is this really necessary?

lkdvos commented 4 months ago

In principle, the required interface to be fully compatible with KrylovKit is the full interface specified by VectorInterface. I think for now though, there are still fallback methods available that will automatically dispatch to LinearAlgebra's axpy! and axpby!, as well as rmul!, whenever these are defined.

I think the easiest is to indeed add VectorInterface as a dependency, and just implement these methods for your custom type. Note that VectorInterface is quite lightweight, specifically for this reason, and as KrylovKit already depends on this it would already be in the implicit dependencies anyways.

An alternative solution is to make use of the fact that KrylovKit already depends on VectorInterface, and do something like:

import KrylovKit: VectorInterface
function VectorInterface.add(x::MyType, y::MyType, a::Number, b::Number)
    # my implementation
end

At the very least in the REPL this seems to work, and does not require an explicit dependency.

Note that I do not know which of these is considered best practice, but I hope to get people's opinion on this matter soon, see this discourse post

Jutho commented 4 months ago

Up till version 0.6, we required axpby! and a few other methods, as you can find in the old docs: https://jutho.github.io/KrylovKit.jl/v0.6/

Since v0.7, we switched to depending on VectorInterface, by which we can support a much larger range of types (such as nested arrays, tuples, …) out of the box. In particular, from VectorInterface we use both the non-mutating functions (scale, add, ...) as well as the possibly mutating methods (scale!!, add!!).

As @lkdvos mentioned, VectorInterface is a dependency of KrylovKit. We did not reexport it in order not to pollute the namespace, but you can indeed import VectorInterface from KrylovKit as mentioned there.

As for which methods to extend, it also depends a bit. For AbstractArray, there are default implementations that use broadcasting, or use BLAS methods (such as axpby!) for suitable types. However, for specific types, more optimised implementations might be preferred. For generic types, there are fallback that use + and *, or test (dynamically!) wheter axpby! and rmul! and friends can be applied. These fallbacks probably have some overhead, and I am not sure if we want to keep those in the future.

We might think about some way to automate generating these method definitions (since there are quite a few) for new types, e.g. using a macro where you just specify your type, whether the data can be mutated, and whether the implementation should be lowering to e.g. LinearAlgebra.axpby! and friends, or to broadcasting, ...). Any input from actual use cases like yours would be very welcome for that.

rveltz commented 4 months ago

Thank you for your helpful comments. In the end I implemented the interface using import KrylovKit: VectorInterface.

rveltz commented 4 months ago

Can you help me a bit with it please? It gives me NaN that were not there before and this block tagging a new version of BifurcationKit

I am doing this

VectorInterface.add!(y::BorderedArray, x::BorderedArray, α::Number = 1, β::Number = 1) = axpby!(β, x, α, y)
VectorInterface.add!!(y::BorderedArray, x::BorderedArray, α::Number = 1, β::Number = 1) = VectorInterface.add!(y, x, α, β)
@inline VectorInterface.scalartype(W::BorderedArray{vectype, Tv1}) where {vectype, Tv1} = eltype(BorderedArray{vectype, Tv1})
VectorInterface.zerovector(b::BorderedArray, S::Type{<:Number} = VectorInterface.scalartype(b)) = 0 * similar(b, S)
VectorInterface.scale(x::BorderedArray, α::Number) = mul!(similar(x), x, α)# α * x
VectorInterface.scale!(y::BorderedArray, x::BorderedArray, α::Number) = throw("")#mul!(y, x, α)
VectorInterface.scale!!(x::BorderedArray, α::Number) = α * x
VectorInterface.scale!!(y::BorderedArray, x::BorderedArray, α::Number) = mul!(y, x, α)
VectorInterface.inner(x::BorderedArray, y::BorderedArray) = dot(x, y)
Jutho commented 4 months ago

Just guessing, but if similar has nans in its unitialized memory, then 0 * … will preserve those. That’s a difference between Julia semantics and BLAS. Maybe use some version of fill!

rveltz commented 4 months ago

OK will try, thank you

lkdvos commented 4 months ago

I think the cleanest way is something like: VectorInterface.zerovector(x::BorderedArray, ::Type{S}) where {S} = BorderedArray(zerovector(x.u, S), zerovector(x.p, S))

In principle you can use a similar approach for all other functions, by defining them recursively, as vectorinterface should work on most AbstractArray, Tuple, Number, ... out of the box

rveltz commented 4 months ago

Ah interesting

lkdvos commented 4 months ago

Did everything work out for you? I think we will not be re-exporting VectorInterface, as it really is meant to be a lightweight dependency and there should be no cost associated to depending on it. According to the discussion on discourse, best practice would be to explicitly add VectorInterface to your dependencies, because now you risk that a breaking vectorinterface update would not be caught by your compat settings. I do admit that this is rather unlikely, and in the end it's a bit a matter of preference too. If you are happy with this, can I close this issue?

rveltz commented 4 months ago

Yeah, I agree, I should add VI to my project.toml which is not what I do at the moment