SciML / SciMLOperators.jl

SciMLOperators.jl: Matrix-Free Operators for the SciML Scientific Machine Learning Common Interface in Julia
https://docs.sciml.ai/SciMLOperators/stable
MIT License
44 stars 10 forks source link
differential-equations julia linear-algebra ode scientific-machine-learning sciml

SciMLOperators.jl

Unified operator interface for SciML.ai and beyond

Join the chat at https://julialang.zulipchat.com #sciml-bridged Global Docs

codecov Build Status

ColPrac: Contributor's Guide on Collaborative Practices for Community Packages SciML Code Style

SciMLOperators is a package for managing linear, nonlinear, time-dependent, and parameter dependent operators acting on vectors, (or column-vectors of matrices). We provide wrappers for matrix-free operators, fast tensor-product evaluations, pre-cached mutating evaluations, as well as Zygote-compatible non-mutating evaluations.

The lazily implemented operator algebra allows the user to update the operator state by passing in an update function that accepts arbitrary parameter objects. Further, our operators behave like AbstractMatrix types thanks to overloads defined for methods in Base, and LinearAlgebra.

Therefore, an AbstractSciMLOperator can be passed to LinearSolve.jl, or NonlinearSolve.jl as a linear/nonlinear operator, or to OrdinaryDiffEq.jl as an ODEFunction. Examples of usage within the SciML ecosystem are provided in the documentation.

Installation

SciMLOperators.jl is a registered package and can be installed via

julia> import Pkg
julia> Pkg.add("SciMLOperators")

Examples

Let M, D, F be matrix-based, diagonal-matrix-based, and function-based SciMLOperators respectively.

N = 4
f(u, p, t) = u .* u
f(v, u, p, t) = v .= u .* u

M = MatrixOperator(rand(N, N))
D = DiagonalOperator(rand(N))
F = FunctionOperator(f, zeros(N), zeros(N))

Then, the following codes just work.

L1 = 2M + 3F + LinearAlgebra.I + rand(N, N)
L2 = D * F * M'
L3 = kron(M, D, F)
L4 = M \ D
L5 = [M; D]' * [M F; F D] * [F; D]

Each L# can be applied to AbstractVectors of appropriate sizes:

p = nothing # parameter struct
t = 0.0     # time

u = rand(N)
v = L1(u, p, t) # == L1 * u

u_kron = rand(N^3)
v_kron = L3(u_kron, p, t) # == L3 * u_kron

For mutating operator evaluations, call cache_operator to generate in-place cache so the operation is nonallocating.

α, β = rand(2)

# allocate cache
L2 = cache_operator(L2, u)
L4 = cache_operator(L4, u)

# allocation-free evaluation
L2(v, u, p, t) # == mul!(v, L2, u)
L4(v, u, p, t, α, β) # == mul!(v, L4, u, α, β)

Roadmap