The AbstractSciMLOperator Interface

abstract type AbstractSciMLOperator{T}

Subtypes of AbstractSciMLOperator represent linear, nonlinear, time-dependent operators acting on vectors, or matrix column-vectors. A lazy operator algebra is also defined for AbstractSciMLOperators.

Mathematical Notation

An AbstractSciMLOperator $L$ is an operator which is used to represent the following type of equation:

\[w = L(u,p,t)[v]\]

where L[v] is the operator application of $L$ on the vector $v$.

Construction and Extension Rules

AbstractSciMLOperator is an interface, not a concrete constructor. New operator types should subtype it with the scalar element type T and must preserve their mathematical action when they participate in lazy algebra.

Required Interface

An AbstractSciMLOperator can be called like a function in the following ways:

• L(v, u, p, t) - Out-of-place application where v is the action vector and u is the update vector
• L(w, v, u, p, t) - In-place application where w is the destination, v is the action vector, and u is the update vector
• L(w, v, u, p, t, α, β) - In-place application with scaling: w = α*(L*v) + β*w

Operator state can be updated separately from application:

• update_coefficients!(L, u, p, t) for in-place operator update
• L = update_coefficients(L, u, p, t) for out-of-place operator update

SciMLOperators also overloads Base.*, LinearAlgebra.mul!, LinearAlgebra.ldiv! for operator evaluation without updating operator state. An AbstractSciMLOperator behaves like a matrix in these methods. Allocation-free methods, suffixed with a ! often need cache arrays. To precache an AbstractSciMLOperator, call the function L = cache_operator(L, input_vector).

Required Interface For Subtypes

A concrete subtype must define Base.size(L) and one of the following application paths:

• Base.:*(L, v) for out-of-place matrix-like application.
• LinearAlgebra.mul!(w, L, v) and, when has_mul!(L) is true, LinearAlgebra.mul!(w, L, v, α, β) for in-place application.
• Base.convert(AbstractMatrix, L) when isconvertible(L) is true.

If the operator state depends on (u, p, t) or accepted keyword arguments, the subtype must implement update_coefficients for out-of-place state updates or update_coefficients! for in-place state updates. Composite operators assume these update methods may be called recursively on every operator returned by getops(L).

Subtypes that need preallocated work arrays for allocation-free application must implement cache_self(L, v) for their own caches, cache_internals(L, v) for child-operator caches, or both. cache_operator(L, v) calls these hooks and downstream solvers may call it before repeated mul! evaluations.

Caching Rules

cache_operator(L, v) may return either L or a cached replacement. A subtype that advertises has_mul!(L) == true must ensure the cached result is ready for repeated mul! calls with compatible vectors. Cache hooks may not change the mathematical action, dimensions, or trait values of the operator.

Trait Rules

Trait functions such as isconstant, islinear, isconvertible, has_concretization, has_mul, has_mul!, has_ldiv, and has_ldiv! are part of the public operator interface. A trait returning true is a promise that the corresponding operation is valid for inputs with compatible sizes. For example, has_mul!(L) means mul!(w, L, v) is available, and has_concretization(L) means either convert(AbstractMatrix, L) or convert(Number, L) can materialize the operator state without changing its mathematical action.

isconstant(L) means repeated calls to update_coefficients[!] are not required to keep L current. islinear(L) means the action is linear in the vector being multiplied; state dependence on (u, p, t) is still allowed for a linear operator.

Keyword Arguments

When an operator accepts keywords during updates, its constructor must record the accepted names with accepted_kwargs. Composite operators forward only those accepted keywords to each component. Extension authors must therefore accept (u, p, t; kwargs...) consistently in every update and application method they advertise.

Standard Actions

The behavior of a SciMLOperator is indistinguishable from an AbstractMatrix. These operators can be passed to linear solver packages, and even to ordinary differential equation solvers. The list of overloads to the AbstractMatrix interface includes, but is not limited to, the following:

• Base: size, zero, one, +, -, *, /, \, ∘, inv, adjoint, transpose, convert
• LinearAlgebra: mul!, ldiv!, lmul!, rmul!, factorize, issymmetric, ishermitian, isposdef
• SparseArrays: sparse, issparse

Multidimensional arrays and batching

SciMLOperator can also be applied to AbstractMatrix subtypes where operator-evaluation is done column-wise.

K = 10
u_mat = rand(N, K)

v_mat = F(u_mat, p, t) # == mul!(v_mat, F, u_mat)
size(v_mat) == (N, K) # true

L can also be applied to AbstractArrays that are not AbstractVecOrMats so long as their size in the first dimension is appropriate for matrix-multiplication. Internally, SciMLOperators reshapes an N-dimensional array to an AbstractMatrix, and applies the operator via matrix-multiplication.

Operator update

This package can also be used to write state-dependent, time-dependent, and parameter-dependent operators, whose state can be updated per a user-defined function. The updates can be done in-place, i.e. by mutating the object, or out-of-place, i.e. in a non-mutating, Zygote-compatible way.

For example,

u = rand(N)
p = rand(N)
t = rand()

# out-of-place update
mat_update_func = (A, u, p, t) -> t * (p * u')
sca_update_func = (a, u, p, t) -> t * sum(p)

M = MatrixOperator(zero(N, N); update_func = mat_update_func)
α = ScalarOperator(zero(Float64); update_func = sca_update_func)

L = α * M
L = cache_operator(L, v)

# L is initialized with zero state
L * v == zeros(N) # true

# update operator state with `(u, p, t)`
L = update_coefficients(L, u, p, t)
# and multiply
L * v != zeros(N) # true

# updates state and evaluates L*v at (u, p, t)
L(v, u, p, t) != zeros(N) # true

The out-of-place evaluation function L(v, u, p, t) calls update_coefficients under the hood, which recursively calls the update_func for each component SciMLOperator. Therefore, the out-of-place evaluation function is equivalent to calling update_coefficients followed by Base.*. Notice that the out-of-place evaluation does not return the updated operator.

On the other hand, the in-place evaluation function, L(w, v, u, p, t), mutates L, and is equivalent to calling update_coefficients! followed by mul!. The in-place update behavior works the same way, with a few <!>s appended here and there. For example,

w = rand(N)
v = rand(N)
u = rand(N)
p = rand(N)
t = rand()

# in-place update
_A = rand(N, N)
_d = rand(N)
mat_update_func! = (A, u, p, t) -> (copy!(A, _A); lmul!(t, A); nothing)
diag_update_func! = (diag, u, p, t) -> copy!(diag, N)

M = MatrixOperator(zero(N, N); update_func! = mat_update_func!)
D = DiagonalOperator(zero(N); update_func! = diag_update_func!)

L = D * M
L = cache_operator(L, v)

# L is initialized with zero state
L * v == zeros(N) # true

# update L in-place
update_coefficients!(L, v, p, t)
# and multiply
mul!(w, v, u, p, t) != zero(N) # true

# updates L in-place, and evaluates w=L*v at (u, p, t)
L(w, v, u, p, t) != zero(N) # true

The update behavior makes this package flexible enough to be used in OrdinaryDiffEq. As the parameter object p is often reserved for sensitivity computation via automatic-differentiation, a user may prefer to pass in state information via other arguments. For that reason, we allow update functions with arbitrary keyword arguments.

mat_update_func = (A, u, p, t; scale = 0.0) -> scale * (p * u')

M = MatrixOperator(zero(N, N); update_func = mat_update_func,
    accepted_kwargs = (:state,))

M(v, u, p, t) == zeros(N) # true
M(v, u, p, t; scale = 1.0) != zero(N)

AbstractSciMLScalarOperator{T} <: AbstractSciMLOperator{T}

Abstract interface for a scalar-valued linear scaling operator.

Interface Rules

Subtypes must provide convert(Number, operator) for their current scalar value. They may implement state updates through update_coefficients and update_coefficients! using the same (u, p, t; kwargs...) contract as AbstractSciMLOperator. Scalar operators act on numbers and arrays, and their addition, multiplication, division, and inversion remain lazy so that later updates affect the composed expression.

Use ScalarOperator to construct a concrete scalar operator. Custom scalar operator types should only claim traits such as has_ldiv when the converted scalar has the corresponding operation for the current state.

update_coefficients(L, u, p, t; kwargs...)

Update the state of L based on u, input vector, p parameter object, t, and keyword arguments. Internally, update_coefficients calls the user-provided update_func method for every component operator in L with the positional arguments (u, p, t) and keyword arguments corresponding to the symbols provided to the operator via kwarg accepted_kwargs.

This method is out-of-place, i.e. fully non-mutating and Zygote-compatible.

Example

using SciMLOperators

mat_update_func = (A, u, p, t; scale = 1.0) -> p * p' * scale * t

M = MatrixOperator(zeros(4,4); update_func = mat_update_func,
                   accepted_kwargs = Val((:scale,)))

L = M + IdentityOperator(4)

u = rand(4)
p = rand(4)
t = 1.0
v = rand(4)

# Update the operator to `(u,p,t)` and apply it to `v`
L = update_coefficients(L, u, p, t; scale = 2.0)
result = L * v

# Or use the interface which separates the update from the application
result = L(v, u, p, t; scale = 2.0)

update_coefficients!(L, u, p, t; kwargs...)

Update in-place the state of L based on u, input vector, p parameter object, t, and keyword arguments. Internally, update_coefficients! calls the user-provided mutating update_func! method for every component operator in L with the positional arguments (u, p, t) and keyword arguments corresponding to the symbols provided to the operator via kwarg accepted_kwargs.

Example

using SciMLOperators

_A = rand(4, 4)
mat_update_func! = (A, u, p, t; scale = 1.0) -> copy!(A, _A)

M = MatrixOperator(zeros(4,4); update_func! = mat_update_func!)

L = M + IdentityOperator(4)

u = rand(4)
p = rand(4)
t = 1.0
v = rand(4)

update_coefficients!(L, u, p, t)
L * v

cache_operator(L, u)

Allocate caches for L for in-place evaluation with u-like input vectors.

concretize(L) -> AbstractMatrix

concretize(L) -> Number

Convert SciMLOperator to a concrete type via eager fusion. This method is a no-op for types that are already concrete.

DEFAULT_UPDATE_FUNC(A, u, p, t)

The default update function for AbstractSciMLOperators, a no-op that leaves the operator state unchanged.

isconstant(_)

Checks if an L's state is constant or needs to be updated by calling update_coefficients.

iscached(L)

Checks whether L has preallocated caches for inplace evaluations.

Check if SciMLOperator L has preallocated cache-arrays for in-place computation.

Checks if size(L, 1) == size(L, 2).

islinear(_)

Checks if L is a linear operator.

isconvertible(L) -> Bool

Checks if L can be cheaply converted to an AbstractMatrix via eager fusion.

has_adjoint(L)

Check if adjoint(L) is lazily defined.

has_expmv(L)

Check if expmv(L, v, t), equivalent to exp(t * A) * v, is defined for Number t, and AbstractArray u of appropriate size.

has_expmv!(L)

Check if expmv!(w, L, v, t), equivalent to mul!(w, exp(t * A), v), is defined for Number t, and AbstractArrays w, v of appropriate sizes.

has_exp(L)

Check if exp(L) is defined lazily defined.

has_mul(L)

Check if L * v is defined for AbstractArray u of appropriate size.

has_mul!(L)

Check if mul!(w, L, v) is defined for AbstractArrays w, v of appropriate sizes.

has_ldiv(L)

Check if L \ v is defined for AbstractArray v of appropriate size.

has_ldiv!(L)

Check if ldiv!(w, L, v) is defined for AbstractArrays w, v of appropriate sizes.

has_concretization(_)

Return whether L can be materialized into a concrete scalar or matrix representation for fallback operations.

Arguments

• L: An operator-like object.

Returns

true when concretize(L) is expected to succeed by calling convert(AbstractMatrix, L) or convert(Number, L), and false otherwise.

Interface Rules

Subtypes of AbstractSciMLOperator should define has_concretization(L) as true only when their current state can be materialized without changing the operator action. Composite operators should return true only when every component needed for the materialization also has concretization.

This trait is intentionally separate from isconvertible(L): an operator can have a correct concrete representation but avoid cheap eager fusion in generic algebra paths.

Examples

using SciMLOperators

A = MatrixOperator([1.0 2.0; 3.0 4.0])
has_concretization(A) # true

F = FunctionOperator((y, x, u, p, t) -> copyto!(y, x), zeros(2), zeros(2);
    isinplace = true, T = Float64, islinear = true)
has_concretization(F) # false

This type indicates to preprocess_update_func to not to filter keyword arguments. Required in implementation of lazy Base.adjoint, Base.conj, Base.transpose.

has_tensor_outer_mul_fast(outer) -> Bool

Return whether outer provides the specialized tensor_outer_mul_fast! contract used by batched TensorProductOperator multiplication.

Developer API

This hook is for extension authors implementing an allocation-free fast path for an outer operator type. Return true only when the corresponding tensor_outer_mul_fast! methods are defined for both the unscaled and scaled call signatures. End users should rely on mul! or TensorProductOperator instead of calling or extending this hook directly.

tensor_outer_mul_fast!(w, outer, C, mi, mo, no, k[, α, β]) -> w

Write the batched outer multiplication used by TensorProductOperator into w without allocating intermediate arrays.

Arguments

• w: destination with mi * mo rows and k columns.
• outer: operator of size (mo, no).
• C: cached intermediate data with mi * no rows and k columns.
• mi, mo, no, k: dimensions derived from the tensor-product factors and the batch size.
• α, β: optional scaling coefficients; the scaled method must compute w = α * outer_product + β * w.

Developer API

Only implement this hook together with has_tensor_outer_mul_fast returning true for the same outer type. Implement both call signatures, preserve the stated destination shape, and return w. This contract exists for package extensions; ordinary callers should use mul! on the enclosing TensorProductOperator.

getcache(_)

Return the cache held by op, or nothing when it holds none.

Defining this method opts op's type into the scratch sharing AddedOperator does between its summands, which is sound because mul! applies them strictly serially and no two of their caches are ever live at once. A summand may then be given the cache of an earlier one whose cache has the identical type and identically sized slots, so only define this for operators where two such caches really are interchangeable — including which slots deliberately alias each other. The nothing default keeps a type out of it entirely.

Two further conditions come with defining it. The first is on the operator's own cache_self: it must treat the buffers it is handed as a template to similar, zero or reshape and never retain a reference to them, since the whole point is that those buffers may afterwards be replaced by an earlier summand's. Every cache_self in this package does. The second is on the caller: applying the summands of an AddedOperator concurrently, by hand, is no longer safe once they share scratch — a cached operator was never safe to apply concurrently, but now its summands are not independently safe either.

A wrapper that holds no scratch of its own — ScaledOperator, AdjointOperator, TransposedOperator — is transparent to all of this and forwards getcache, update_cache and adopt_cache to the operator it wraps. A new wrapper of that kind should do the same; without it the operator underneath drops out of sharing even though its cache would have qualified.

See also adopt_cache, which additionally spares the operator from allocating a cache it would only discard.

update_cache(op, new_cache)

Replace op's cache with new_cache. Only ever called with a new_cache of exactly the same type as the one it replaces, so op's type is unchanged.

adopt_cache(_, cache, v)

Return op set up for in-place use with v using cache — the cache of an operator already established as interchangeable with it — instead of allocating buffers of its own. Return nothing to decline, in which case op is cached the ordinary way and then offered the same cache after the fact, which costs an allocation that is immediately discarded. What comes back must report cache as its own — getcache on it is checked against what was handed over, and anything else is treated as having declined, since an implementation that quietly allocated instead would otherwise skip the ordinary path too.

Declining is the default, because there is no implementation that is correct for every operator: one that follows the cache_self/cache_internals split opts in with

adopt_cache(op::MyOperator, cache, v) = cache_internals(update_cache(op, cache), v)

whereas one that defines its own cache_operator must instead do whatever that does, minus the allocation — the line above would silently skip its internal caching entirely.

Defining this asks more of an operator than getcache alone. getcache only claims that two caches of the same type with the same slot sizes are interchangeable, and that is checked against buffers that already exist. This additionally claims that what cache_self builds is a function of the operator's type, its operands' sizes and v, since that is what the caller compares before any buffer exists. Leave it undefined when the cache depends on anything else — FunctionOperator does, because it sizes buffers from traits.sizes, which size does not expose.