diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml new file mode 100644 index 0000000..9e47014 --- /dev/null +++ b/.github/workflows/Documentation.yml @@ -0,0 +1,37 @@ +name: Documentation +on: + workflow_dispatch: + push: + branches: + - main + tags: '*' + pull_request: +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + docs: + name: Build and deploy documentation + runs-on: ubuntu-latest + permissions: + actions: write + contents: write + pull-requests: read + statuses: write + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 + with: + version: '1' + - uses: julia-actions/cache@v3 + - name: Install docs dependencies + run: | + julia --project=docs -e ' + using Pkg + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + run: julia --project=docs docs/make.jl diff --git a/.gitignore b/.gitignore index 0f84bed..95731a5 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ *.jl.cov *.jl.mem /Manifest.toml +/docs/Manifest.toml +/docs/build/ diff --git a/README.md b/README.md index b6337dd..ce589c2 100644 --- a/README.md +++ b/README.md @@ -6,114 +6,49 @@ | [![][docs-stable-img]][docs-stable-url] [![][docs-dev-img]][docs-dev-url] | [![][aqua-img]][aqua-url] [![CI][github-img]][github-url] [![][codecov-img]][codecov-url] | [![license][license-img]][license-url] | [docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg -[docs-dev-url]: https://jutho.github.io/VectorInterface.jl/latest +[docs-dev-url]: https://quantumkithub.github.io/VectorInterface.jl/dev [docs-stable-img]: https://img.shields.io/badge/docs-stable-blue.svg -[docs-stable-url]: https://jutho.github.io/VectorInterface.jl/stable +[docs-stable-url]: https://quantumkithub.github.io/VectorInterface.jl/stable -[github-img]: https://github.com/Jutho/VectorInterface.jl/workflows/CI/badge.svg -[github-url]: https://github.com/Jutho/VectorInterface.jl/actions?query=workflow%3ACI +[github-img]: https://github.com/QuantumKitHub/VectorInterface.jl/workflows/CI/badge.svg +[github-url]: https://github.com/QuantumKitHub/VectorInterface.jl/actions?query=workflow%3ACI [aqua-img]: https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg [aqua-url]: https://github.com/JuliaTesting/Aqua.jl -[codecov-img]: https://codecov.io/gh/Jutho/VectorInterface.jl/branch/main/graph/badge.svg -[codecov-url]: https://codecov.io/gh/Jutho/VectorInterface.jl +[codecov-img]: https://codecov.io/gh/QuantumKitHub/VectorInterface.jl/branch/main/graph/badge.svg +[codecov-url]: https://codecov.io/gh/QuantumKitHub/VectorInterface.jl [license-img]: http://img.shields.io/badge/license-MIT-brightgreen.svg?style=flat [license-url]: LICENSE.md --- -An Julia interface proposal for vector-like objects. +A small Julia interface for vector-like objects. +`VectorInterface.jl` proposes a minimal, unified set of methods that any type representing an element of a vector space should support, so that generic algorithms (Krylov methods, ODE integrators, gradient optimizers, ...) can be written once and reused across `Number`s, `AbstractArray`s, nested arrays, tuples, named tuples, and arbitrary user-defined types — without conflating vector-space operations with the iteration, container, or `LinearAlgebra` interfaces. -### Motivation +## Installation -This package proposes a Julia interface for vector-like objects, i.e. a small set of methods that should be supported by types representing objects that can be interpreted as living in a vector space. +```julia +julia> using Pkg; Pkg.add("VectorInterface") +``` -Having such a unified interface is useful, because this small collection of methods can be used to write a large collection of algorithms in a universal and agnostic manner. Unfortunately, the current interfaces existing in the Julia ecosystem are not adequate for this purpose. +## Quick example -### Context +```julia +using VectorInterface -Recall the basic properties of vectors in a given vector space. A vector space is a set of objects, called vectors, with the basic properties that: +v = [1.0, 2.0, 3.0] +w = [4.0, 5.0, 6.0] -* Vectors can be added (which is commutative, there is a neutral element and inverses) -* The neutral element is called the zero vector -* Vectors can be rescaled with a scalar coefficient, taking values in an underlying scalar field. (and there are all kind of relations with vector addition, such as distributivity, ...) -Together, these two operations give rise to the construction of linear combinations. On top of those properties, many vector spaces used in practical applications have more structure. -* Often there is an inner product, or at least a norm. -* Finally, many useful vector spaces admit a finite basis, such that any vector can be written as a linear combination of a finite number of basis vectors. +scalartype(v) # Float64 +zerovector(v) # [0.0, 0.0, 0.0] +scale!!(v, 2.0) # tries in place: [2.0, 4.0, 6.0] +add!!(v, w, 0.5) # v <- v + 0.5 * w +inner(v, w) # conjugate-linear inner product +norm(v) # re-exported from LinearAlgebra +``` -Many quantities in science and engineering (and thus in scientific computing) behave as vectors, typically with real or complex numbers as underlying field. Even for quantities that do not (but rather belong to a manifold), their derivatives (tangents and cotangents) do, which is important for automatic differentiation. +Every mutating operation comes in three variants: `f` (allocating), `f!` (strict in-place), and `f!!` (try in-place, fall back to allocating), following the convention from [`BangBang.jl`](https://github.com/JuliaFolds2/BangBang.jl). -More importantly, many algorithms can be formulated using the basic operations (linear combinations and inner products), in combination with some recipe for spitting out new vectors (e.g. applying a linear map): gradient optimization, ODE integrators, Krylov methods, ... - -### Current situation and problems in Julia - -The most elementary Julia type that acts as a vector is `Vector{T<:Number}`, and by extension (almost *) any subtype of `AbstractArray{T<:Number}`. However, many other types which are not subtypes of `AbstractArray` also are vectors conceptually, e.g. types for representing functions (ApproxFun ecoystem), `Tangent` types in the AD ecosystem, and many more. - -However, I have several gripes with the Julia interface (or lack thereof) to access the basic vector operations, which makes writing generic algorithms hard. In particular: - -1. Vector addition and vector rescaling almost always works simply as `v + w` and `v * α`. However, for efficiency, we often want to be able to do these operations in place. For instances of `AbstractArray`, we then need to be `using LinearAlgebra` (which also includes a lot of stuff that we might not need). For vector rescaling within the `AbstractArray` family, we can use `rmul!(v, α)` or `lmul!(α, v)` or `mul!(w, v, α)`, but that interface is conflated with the concept of matrix multiplication etc, whereas these are two very different concepts. There used to be a `scale!` method (in analogy with a corresponding method in BLAS level 1) in older Julia 0.x versions (though it was also not perfect, and was also used for matrix times diagonal matrix). For vector addition on the other hand, the options for in-place methods are the cryptic `axpy!` and `axpby!` methods (referring to their BLAS level 1 analogues), with thus a very un-Julian interface (the vector that is modified is the final rather than the first argument). - -2. In programming, rather than the scalar field (reals or complex), we of course want to know the specific scalar type, with which the vectors can be natively rescaled. For an instance `v` of `AbstractArray{T<:Number}`, this is the type `T` and it can be obtained as `eltype(v)`. However, because `eltype` is also used by the iteration interface, other types which might have an iteration behaviour that is distinct from their vector-like behaviour, cannot overload `eltype` for both purposes. An example from `Base` would be a nested array, e.g. the type `Vector{Vector{T<:Number}}` still constitutes a set of vectors with scalar type `T`, but `eltype` equal to `Vector{T}` - -3. To get the zero vector associated with some vector `v`, we can use `zero(v)`, which is fine as an interface, as it is defined to be the neutral element with respect to `+`. However, `zero` on e.g. a nested array will fail because of how it is implemented. Furthermore, to make a zero vector in place, you could use `fill!(v, 0)` for `v::AbstractArray{T<:Number}`, but that is a very `AbstractArray` specific interface. The only more general solution is to restort to scaling by zero, e.g. using `rmul!(v, 0)`, if available. - -4. Closely related to the previous two points, we often want to be able to create equivalent vectors but with a modified scalar type. For vectors in `AbstractArray{T<:Number}`, there is `similar(v, T′)` for this, but that is again a very array-specific method, and fails once more for e.g. nested arrays. - -5. Most (but not all) vector-like objects in programming belong to a finite-dimensional vector space. For `v::AbstractArray{T<:Number}`, this dimension is given by `length(v)`, but again this interface is also used for the iteration length, and so new types might face an incompatibility as with `eltype`. And for structured arrays, `length(v)` might also not express the vector space dimension, for e.g. `UpperTriangular{T<:Number}`, the natural vector space dimension is `n*(n+1)/2`, not `n*n`. - -6. The inner product and norm corresponds to the Julia methods `LinearAlgebra.dot` and `LinearAlgebra.norm`. Unlike in some of the previous points, `dot` and `norm` natively support nested arrays. However, `dot` is so loose in its implementation, that it also happily computes an inner product between things which are probably not vectors from the same vector space, such as `dot( (1, 2, [3, 4]), [[1], (2,), (3,4)])`. In particular, `dot` and `norm` also accept tuples, whereas tuples do not behave as vectors with respect to the previous methods (`+`, `*`, `zero`). - -In summary, the main problem is that there actually is no formal standardized vector interface in Julia, despite its broad potential applicability for writing very generic algorithms. There are standardized interfaces for containers (`AbstractArray`) and for iterators, which have become conflated with a hypothetical vector interface. - -### Existing solutions - -Different ecosystems have responded to this hiatus in different ways. Several Krylov and optimization packages merely restrict their applicability to instances of `(Abstract)Array{T<:Number}` or even simply `Vector{T<:Number}` (like their Fortran and C analogues would probably do). The DifferentialEquations.jl ecosystem does more or less the same, restricting to `AbstractArray` (if I remember correctly), but provides a bunch of packages such as `RecursiveArrayTools.jl` and `ArrayInterface.jl` to accommodate for more complex use cases. Finally, the AD ecosystem (Zygote.jl and ChainRules.jl) use custom `Tangent` types for which they define the necessary operations, using a lot of internal machinery to destructure custom types. - -Forcing everything to be a subtype of `AbstractArray` is an unsatisfactory solution. Some vector like objects might not be naturally represented with respect to a basis, and thus have no notion of indexing, and might not even be finite-dimensional. The `AbstractArray` or container interface is and should be distinct from a general vector (space) interface. - -### New solution -With VectorInterface.jl, I have tried to create a simple package to resolve my gripes. As I hope that I am not alone with those, I would like this to be useful for the community and could eventually evolve into a standardized interface. Therefore, I would very much value comments. Everything is up for bikeshedding. I tried to come up with a design which is compatible with `LinearAlgebra` (e.g. not stealing names) and does not commit type piracy. Currently, VectorInterface.jl provides the following minimalistic interface: - -* `scalartype(v)`: accesses the native scalar type `T<:Number` of a vector-like object; also works in the type domain (i.e. `scalartype(typeof(v))`). -* `zerovector(v)`, `zerovector!(v)` and `zerovector!!(v)`: produce a zero vector of the same type as `v`; the second method tries to do this in-place for mutable types, whereas the third method is inspired by BangBang.jl and tries to do it in-place when possible, and out of place otherwise. - - *Comment: Ideally, the `zerovector` functionality would be provided by `Base.zero`.* - -* `zerovector(v, S<:Number)` creates a zero vector of similar type, but with a modified scalar type that is now given by `S`. `zerovector!(v)` sets `v` equal to the zero vector; in that case the scalar type cannot be changed. Finally `zerovector!!(v, S)` works, but reduces to `zerovector(v, S)` if `scalartype(v)` is not `S`. - - *Comment: Given that there is a tendency to zero out uninitialized memory, I think it is fine to merge the concept of constructing a new vector with different scalar type with that of constructing the zero vector.* - -* `scale(v, α)`, `scale!(v, α)` and `scale!!(v, α)` rescale the vector `v` with the scalar coefficient `α`. The second method tries to do this in place, but will fail if `α` cannot be converted to `scalartype(v)` (or if `v` contains immutable contents), whereas the third method is the BangBang-style unified solution. There is also `scale!(w, v, α)` and `scale!!(w, v, α)` to rescale `v` out of place, but storing the result in `w`. - -* `add(w, v, [, α = 1, β = 1])`, `add!(w, v, [, α = 1, β = 1])` and `add!!(w, v, [, α = 1, β = 1]))` compute `w * β + v * α`, where (by now self-explanatory) the second method stores the result in `w` (and errors if not possible), and the third method tries to store in `w` but doesn't error. - -* `inner(v, w)` works almost equivalent to `dot(v, w)`, is sometimes a bit more efficient, but also more strict in what arguments it allows. - -* `norm(v)` simply reexports `LinearAlgebra.norm` - -These methods are implemented for instances `v` of type `<:Number` (scalars are also vectors over themselves) and `<:AbstractArray` (both `<:AbstractArray{<:Number}` and nested array). - -In addition, the interface is currently also defined for tuples and named tuples, again with arbitrary nesting. So instances of e.g. `Vector{NTuple{3,Matrix{Float64}}}` are currently also supported. - -Furthermore, in-place methods will work recursively so as to try to maximally work in place. What I mean by that is, if you have nested vectors, say `v = [[1., 2.], [3., 4.]]`, then `rmul!(v, 0.5)` will keep the outer array, but will work on the inner arrays using regular `*` multiplication, and will thus allocate two new inner arrays in this example. In contrast, `scale!(v, 0.5)` does work in-place on the inner arrays and is free of allocations. - -Similarly, for `v = ((1., 2.), [3., 4.])`, `scale!!(v, 0.5)` could of course not work in-place on the outer tuple or inner tuple, but would still work in-place on the inner array. Hence, the return value of `scale!!(v, 0.5)` is of course `((0.5, 1.), [1.5, 2.])`, but after his operation, `v` would be `((1., 2.), [1.5, 2.])`. - -### Open questions - -Upon development, there were several decisions to be made, about which I was not sure. -* Should tuples actually be supported? In Base, they are not treated as vectors, though they are supported by `dot` and `norm`. - -> Currently, tuples and even named tuples are supported (which was a request) -* Should there be some `vectordim` function (name up to debate) to probe the vector space dimension? - -> No such function is currently defined -* Should `LinearAlgebra.dot` be exported? Or should `inner` just become `dot` and should and is the looseness of `dot` of no concern? - -> `dot` is currently not exported; `inner` is the preferred method for inner products -* Should there be fallbacks in place for user types that did already implement `rmul!`, `mul!`, `axpy!`, `axpby!` and `dot` (the latter relating to the previous question)? - -> Currently, there is a fallback to these methods, but with the intention to remove these fallbacks in future versions. - -All thoughts and comments are very welcome. - -(*) There is one (actually two) subtype of `AbstractArray` in `LinearAlgebra` that does not behave as a vector, in the sense that its instances cannot represent the zero vector or cannot be rescaled without changing its type, namely `UnitUpperTriangular` and `UnitLowerTriangular`. The fixed unit diagonal prevents these types from constituting a vector space. It seems like the unit diagonal also poses issues for broadcasting, as operations that are preserving ones are much more rare than operations preserving zeros (which is necessary for any structured or unstructured sparseness). diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..feeb91a --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,7 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" + +[compat] +Documenter = "1" +julia = "1.10" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..ae3105c --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,29 @@ +using Documenter +using VectorInterface + +DocMeta.setdocmeta!(VectorInterface, :DocTestSetup, :(using VectorInterface); recursive = true) + +makedocs( + sitename = "VectorInterface.jl", + modules = [VectorInterface], + authors = "Jutho Haegeman and contributors", + pages = [ + "Home" => "index.md", + "Manual" => [ + "Motivation" => "man/motivation.md", + "Interface" => "man/interface.md", + ], + "API reference" => "api.md", + ], + format = Documenter.HTML( + prettyurls = get(ENV, "CI", nothing) == "true", + canonical = "https://quantumkithub.github.io/VectorInterface.jl", + ), + checkdocs = :exports, +) + +deploydocs( + repo = "github.com/QuantumKitHub/VectorInterface.jl.git", + devbranch = "main", + push_preview = true, +) diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 0000000..46cafba --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,14 @@ +```@meta +CurrentModule = VectorInterface +``` + +# API reference + +```@index +Pages = ["api.md"] +``` + +```@autodocs +Modules = [VectorInterface] +Order = [:type, :function] +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..e592930 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,36 @@ +```@meta +CurrentModule = VectorInterface +``` + +# VectorInterface.jl + +A small Julia interface for vector-like objects. + +`VectorInterface.jl` proposes a minimal set of methods that any type representing an element of a vector space should support. +Code written against this interface can then operate uniformly on `Number`s, `AbstractArray`s, nested arrays, tuples, named tuples, and arbitrary user-defined types — without conflating vector-space operations with the iteration, container, or `LinearAlgebra` interfaces. + +## Installation + +```julia +julia> using Pkg; Pkg.add("VectorInterface") +``` + +## What's provided + +- A scalar-type query: [`scalartype`](@ref). +- Construction of zero vectors: [`zerovector`](@ref), [`zerovector!`](@ref), [`zerovector!!`](@ref). +- Scalar multiplication: [`scale`](@ref), [`scale!`](@ref), [`scale!!`](@ref). +- Linear combinations: [`add`](@ref), [`add!`](@ref), [`add!!`](@ref). +- Inner product and norm: [`inner`](@ref), and `norm` is re-exported from `LinearAlgebra`. +- Singleton helpers [`One`](@ref) and [`Zero`](@ref) for hard-coded coefficients in linear combinations. + +Each `!`-suffixed method modifies its first argument in place. +Each `!!`-suffixed method tries to do so but falls back to allocating when in-place updates are not possible (e.g. for immutable types or incompatible scalar types). +This convention follows [`BangBang.jl`](https://github.com/JuliaFolds2/BangBang.jl). + +## Contents + +```@contents +Pages = ["man/motivation.md", "man/interface.md", "api.md"] +Depth = 2 +``` diff --git a/docs/src/man/interface.md b/docs/src/man/interface.md new file mode 100644 index 0000000..bd7927e --- /dev/null +++ b/docs/src/man/interface.md @@ -0,0 +1,188 @@ +```@meta +CurrentModule = VectorInterface +DocTestSetup = quote + using VectorInterface +end +``` + +# Interface + +`VectorInterface.jl` exports a small, fixed set of methods that together constitute the vector-space interface. +This page walks through each method group, with examples. + +For every mutating operation, three variants are provided: + +- `f(args...)` — pure, allocating, always returns a new object. +- `f!(args...)` — strict in-place, errors if the in-place update is not possible (e.g. the target is immutable or the scalar types are incompatible). +- `f!!(args...)` — tries in-place, falls back to allocating when needed. + This follows the convention of [`BangBang.jl`](https://github.com/JuliaFolds2/BangBang.jl) and is the right choice when writing generic algorithms that should be both efficient on mutable inputs and correct on immutable ones. + +## Scalar type + +```@docs; canonical=false +scalartype +``` + +`scalartype` returns the type of scalars over which the vector-like object behaves as a vector. +For an `AbstractArray{T <: Number}` this is `T`, and for a nested array `Vector{Vector{Float64}}` it is still `Float64` (whereas `eltype` would return `Vector{Float64}`). + +```jldoctest +julia> scalartype([1.0, 2.0, 3.0]) +Float64 + +julia> scalartype(Vector{Vector{Float64}}) +Float64 +``` + +Custom types should implement the type-domain method `scalartype(::Type{MyType})`. +The value-domain method `scalartype(x)` is defined generically as `scalartype(typeof(x))`. + +## Zero vector + +```@docs; canonical=false +zerovector +zerovector! +zerovector!! +``` + +`zerovector(v)` returns a new zero vector of the same type as `v`. +An optional second argument `S <: Number` overrides the scalar type of the result: + +```jldoctest +julia> zerovector([1.0, 2.0, 3.0]) +3-element Vector{Float64}: + 0.0 + 0.0 + 0.0 + +julia> zerovector([1.0, 2.0, 3.0], ComplexF64) +3-element Vector{ComplexF64}: + 0.0 + 0.0im + 0.0 + 0.0im + 0.0 + 0.0im +``` + +`zerovector!(v)` zeros `v` in place, and `zerovector!!(v)` does the same when possible and otherwise allocates. +For nested vectors, the in-place variants recurse, so `zerovector!(v)` on a `Vector{Vector{Float64}}` does not allocate any new inner arrays. + +## Scaling + +```@docs; canonical=false +scale +scale! +scale!! +``` + +`scale(v, α)` rescales the vector `v` by the scalar `α` and returns a new object. +`scale!(v, α)` and `scale!(w, v, α)` perform the rescaling in place (into `v` or into `w`, respectively), and require the scalar types to be compatible. +`scale!!` tries to do the operation in place and falls back to allocating otherwise. + +```jldoctest +julia> v = [1.0, 2.0, 3.0]; + +julia> scale(v, 2.0) +3-element Vector{Float64}: + 2.0 + 4.0 + 6.0 + +julia> scale!(v, 0.5); v +3-element Vector{Float64}: + 0.5 + 1.0 + 1.5 +``` + +A key design point: `scale!` on a **nested** vector recurses into the inner vectors and operates fully in place, in contrast to `LinearAlgebra.rmul!` which allocates new inner arrays. +For `v = [[1.0, 2.0], [3.0, 4.0]]`, `scale!(v, 0.5)` does not allocate. + +## Linear combinations + +```@docs; canonical=false +add +add! +add!! +``` + +`add(y, x, α, β)` computes `y * β + x * α` and returns a new object. +The in-place variants `add!`/`add!!` overwrite `y` with the result. +The default coefficients are [`One()`](@ref), so the short forms `add(y, x)` and `add(y, x, α)` compute `y + x` and `y * 1 + x * α` respectively. + +```jldoctest +julia> add([1.0, 2.0], [10.0, 20.0]) +2-element Vector{Float64}: + 11.0 + 22.0 + +julia> add([1.0, 2.0], [10.0, 20.0], 0.1) +2-element Vector{Float64}: + 2.0 + 4.0 + +julia> add([1.0, 2.0], [10.0, 20.0], 0.5, 2.0) +2-element Vector{Float64}: + 7.0 + 14.0 +``` + +Custom types should implement only the four-argument form. +When desired, the `One()`-coefficient case can be specialized for efficiency by dispatching on `α::One` and/or `β::One`. + +## Inner product and norm + +```@docs; canonical=false +inner +``` + +`inner(x, y)` is similar to `LinearAlgebra.dot` but is stricter about the arguments it accepts: both `x` and `y` must be elements of the same vector space, with compatible scalar types and shapes. +The norm is re-exported directly from `LinearAlgebra`: + +```jldoctest +julia> inner([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) +32.0 + +julia> norm([3.0, 4.0]) +5.0 +``` + +For complex vectors, `inner(x, y)` is conjugate-linear in its first argument, matching the mathematical convention. + +## Hard-coded coefficients: `One` and `Zero` + +```@docs; canonical=false +One +Zero +``` + +`One` and `Zero` are singleton subtypes of `Number` used to represent hard-coded constant coefficients in linear combinations. +They allow methods like [`add`](@ref) to dispatch on a unit coefficient at compile time, avoiding unnecessary multiplications. +They are the default values for the `α` and `β` coefficients in [`add`](@ref), [`add!`](@ref), and [`add!!`](@ref). + +## Supported types + +Out of the box, `VectorInterface.jl` provides implementations for: + +- `<:Number` — scalars are also vectors over themselves. +- `<:AbstractArray` — both `AbstractArray{<:Number}` and nested arrays, with recursive in-place behaviour for the latter. +- `Tuple` and `NamedTuple` with a homogeneous element scalar type, again with arbitrary nesting. + For example, `Vector{NTuple{3, Matrix{Float64}}}` is supported. + +A general fallback exists for types that already implement `LinearAlgebra.rmul!`, `axpy!`, `axpby!`, `mul!`, and `LinearAlgebra.dot`, but this fallback emits a warning and is intended only as a transitional measure that may be removed in a future release. +New types should implement the `VectorInterface.jl` methods directly. + +## Implementing the interface for a new type + +A new vector-like type `T` should implement, at a minimum: + +- `scalartype(::Type{T})` +- `zerovector(x::T, ::Type{S})` where `S <: Number` +- `scale(x::T, α::Number)` +- `add(y::T, x::T, α::Number, β::Number)` +- `inner(x::T, y::T)` +- `LinearAlgebra.norm(x::T)` + +If the type is mutable, it should additionally implement the `!`-variants (`zerovector!`, `scale!`, `add!`). +The `!!`-variants then have generic fallbacks: they delegate to the `!`-variant when the type is mutable and to the non-mutating variant otherwise. + +The `MinimalVec` type included in the package's source (`src/minimalvec.jl`) is a worked example: it wraps an `AbstractVector` and implements exactly the `VectorInterface.jl` API, deliberately excluding all other `AbstractArray` methods. +It is useful both as a reference and as a test harness for checking that an algorithm depends only on the minimal interface. diff --git a/docs/src/man/motivation.md b/docs/src/man/motivation.md new file mode 100644 index 0000000..4ca009e --- /dev/null +++ b/docs/src/man/motivation.md @@ -0,0 +1,88 @@ +# Motivation + +This page explains why `VectorInterface.jl` exists, what gap it fills in the Julia ecosystem, and why existing interfaces (`AbstractArray`, `LinearAlgebra`, iteration) are not adequate for writing generic algorithms that operate on vector-like objects. + +## What is a vector? + +Recall the basic properties of vectors in a given vector space. +A vector space is a set of objects — called vectors — with the following structure: + +- Vectors can be **added**, and this addition is commutative, has a neutral element (the **zero vector**), and admits inverses. +- Vectors can be **rescaled** by a scalar coefficient, taken from an underlying scalar field. + Scalar multiplication interacts with addition through the usual distributivity laws. + +Together, these two operations give rise to the construction of **linear combinations**. +Beyond this minimum structure, many vector spaces used in practice carry additional structure: + +- An **inner product**, or at least a **norm**. +- A finite **basis**, so that any vector can be written as a finite linear combination of basis vectors. + +Many quantities in science and engineering — and therefore in scientific computing — behave as vectors, typically over the real or complex numbers. +Even for quantities that do not (but rather belong to a manifold), their derivatives (tangents and cotangents) do, which makes a vector interface important for automatic differentiation. + +More importantly, **many algorithms can be formulated using just these basic operations**: gradient optimization, ODE integrators, Krylov methods, and so on. +A standardized vector interface enables a single implementation of such an algorithm to be applied to any vector-like type. + +## Problems with the current Julia situation + +The most elementary Julia type that acts as a vector is `Vector{T<:Number}`, and by extension almost any subtype of `AbstractArray{T<:Number}`. +However, many other types which are not subtypes of `AbstractArray` are also vectors conceptually — types representing functions (the `ApproxFun` ecosystem), `Tangent` types in the automatic-differentiation ecosystem, and many more. + +The Julia interface to access the basic vector operations has several rough edges that make writing generic algorithms hard: + +1. **Vector addition and rescaling don't have a clean in-place interface.** + `v + w` and `v * α` work, but for efficiency we often want in-place operations. + For `AbstractArray`, this requires `using LinearAlgebra`, which pulls in a lot of additional functionality. + Within `LinearAlgebra`, scaling is done via `rmul!(v, α)` / `lmul!(α, v)` / `mul!(w, v, α)` — an interface conflated with matrix multiplication, although scaling a vector and multiplying by a matrix are conceptually very different operations. + For in-place addition, the only options are the cryptic `axpy!` / `axpby!`, whose argument order (the modified vector is last) is not Julian. + +2. **`eltype` is overloaded between vectors and iterators.** + For an instance `v` of `AbstractArray{T<:Number}`, the scalar type is `T = eltype(v)`. + But `eltype` is also part of the iteration interface, so types whose iteration behaviour differs from their vector-like behaviour cannot consistently overload `eltype`. + A simple example from `Base` is a nested array: the type `Vector{Vector{T<:Number}}` still represents vectors with scalar type `T`, but `eltype` reports `Vector{T}`. + +3. **`zero` and `fill!` don't compose well.** + `zero(v)` is fine in principle — it returns the additive neutral element — but the default implementation fails on e.g. nested arrays. + For an in-place zero, `fill!(v, 0)` works for `v::AbstractArray{T<:Number}` but is a very `AbstractArray`-specific interface. + The most general workaround is to scale by zero with `rmul!`, when available. + +4. **`similar` is array-specific.** + Constructing an equivalent vector with a modified scalar type is done via `similar(v, T′)` for `AbstractArray`s, but again this fails for nested arrays and is array-shaped in spirit. + +5. **`length` is overloaded.** + For `v::AbstractArray{T<:Number}`, the vector-space dimension is `length(v)`, but `length` is also part of the iteration interface. + New types may face an `eltype`-style incompatibility. + For structured arrays, `length(v)` may not even reflect the vector-space dimension — for `UpperTriangular{T<:Number}`, the natural dimension is `n*(n+1)/2`, not `n*n`. + +6. **`dot` is too permissive and `norm` lives in `LinearAlgebra`.** + Inner products and norms are `LinearAlgebra.dot` and `LinearAlgebra.norm`. + Unlike some of the previous issues, these natively support nested arrays. + But `dot` happily computes an inner product between things which are probably not vectors from the same vector space, e.g. `dot((1, 2, [3, 4]), [[1], (2,), (3, 4)])`. + In particular `dot` and `norm` accept tuples, even though tuples do not behave as vectors with respect to `+`, `*`, or `zero`. + +In summary: there is no formal standardized vector interface in Julia, even though one has broad applicability for writing generic algorithms. +The interfaces for **containers** (`AbstractArray`) and **iterators** are well defined, but they have become conflated with a hypothetical vector interface. + +## Existing approaches in the ecosystem + +Different ecosystems have responded to this gap in different ways: + +- Many Krylov and optimization packages restrict their applicability to `(Abstract)Array{T<:Number}` or even simply `Vector{T<:Number}`, much like their Fortran and C analogues. +- The DifferentialEquations.jl ecosystem also restricts to `AbstractArray`, with helpers such as `RecursiveArrayTools.jl` and `ArrayInterface.jl` to accommodate more complex use cases. +- The automatic-differentiation ecosystem (Zygote.jl and ChainRules.jl) uses custom `Tangent` types for which the necessary operations are defined, relying on internal machinery to destructure custom types. + +Forcing everything to be a subtype of `AbstractArray` is an unsatisfactory solution. +Some vector-like objects are not naturally represented with respect to a basis, have no notion of indexing, and may not even be finite-dimensional. +The `AbstractArray` (container) interface is — and should be — distinct from a vector-space interface. + +## How `VectorInterface.jl` addresses this + +`VectorInterface.jl` provides a small set of methods that together constitute a minimal vector-space interface. +The design goals are: + +- **Compatibility with `LinearAlgebra`:** no name clashes, no type piracy. +- **Separation of concerns:** the vector interface is independent of the container, iteration, and matrix-multiplication interfaces. +- **In-place, out-of-place, and "try-in-place" variants** for each mutating operation, distinguished by the `!` and `!!` suffixes. +- **Recursive in-place behaviour** for nested vectors: a `scale!(v, α)` on a nested array works in-place all the way down, freeing the caller of allocations that `rmul!` would still incur. + +See the [Interface](@ref) page for the full set of methods, with examples. diff --git a/test/issues.jl b/test/issues.jl index a3ff3e2..959f85a 100644 --- a/test/issues.jl +++ b/test/issues.jl @@ -10,7 +10,7 @@ using TestExtras @test @constinferred add(fill(0.5), fill(0.8), 0.3f0, 0.25im) isa Array{ComplexF64, 0} end -# came up here: https://github.com/Jutho/VectorInterface.jl/issues/20 +# came up here: https://github.com/QuantumKitHub/VectorInterface.jl/issues/20 @testset begin @test promote_type(Bool, One) === promote_type(Bool, Zero) === promote_type(One, Bool) === promote_type(Zero, Bool) === Bool