Unified operator interface for SciML.ai and beyond
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.
SciMLOperators.jl is a registered package and can be installed via
julia> import Pkg
julia> Pkg.add("SciMLOperators")
Let M, D, F be matrix-based, diagonal-matrix-based, and function-based
SciMLOperators respectively.
N = 4
f(v, u, p, t) = v .* u
f(w, v, u, p, t) = w .= v .* u
M = MatrixOperator(rand(N, N))
D = DiagonalOperator(rand(N))
# Fix: Specify that we're providing a parameter placeholder
F = FunctionOperator(f, zeros(N), zeros(N); p=nothing, isconstant=true)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 - must be nothing or compatible with the FunctionOperator
t = 0.0 # time
u = rand(N) # update vector
v = rand(N) # action vector
w = zeros(N) # output vector
# Out-of-place evaluation
result1 = L1(v, u, p, t) # L1 acting on v after updating with u
# In-place evaluation (after caching)
L1_cached = cache_operator(L1, v)
L1_cached(w, v, u, p, t) # Result stored in w
# For tensor product operators
v_kron = rand(N^3)
u_kron = rand(N^3) # Update vector for tensor product
w_kron = zeros(N^3) # Output vector for tensor product
# Cache the tensor product operator with the correct sized vector
L3_cached = cache_operator(L3, v_kron)
# Out-of-place evaluation (action vector v_kron, update vector u_kron)
result_kron = L3_cached(v_kron, u_kron, p, t)
# In-place evaluation (output to w_kron)
L3_cached(w_kron, v_kron, u_kron, p, t)For mutating operator evaluations, call cache_operator to generate
in-place cache so the operation is nonallocating.
α, β = rand(2)
# Allocate and cache operators first
L2 = cache_operator(L2, v)
L4 = cache_operator(L4, v)
# Allocation-free evaluation with separate update and action vectors
w = zeros(N)
L2(w, v, u, p, t) # w = L2 * v
L4(w, v, u, p, t) # w = L4 * v
# In-place evaluation with scaling: w = α*(L*v) + β*w
result_w = rand(N) # Start with random vector
L2(result_w, v, u, p, t, α, β) # result_w = α*(L2*v) + β*result_w
L4(result_w, v, u, p, t, α, β) # result_w = α*(L4*v) + β*result_w