Premade SciMLOperators

Direct Operator Definitions

SciMLOperators.ScalarOperatorType
ScalarOperator(val; update_func, accepted_kwargs)

Represents a linear scaling operator that may be applied to a Number, or an AbstractArray subtype. Its state is updated by the user-provided update_func during operator evaluation (L([w,] v, u, p, t)), or by calls to update_coefficients[!]. Both recursively call the update function, update_func which is assumed to have the signature:

update_func(oldval::Number, u, p, t; <accepted kwargs>) -> newval

The set of keyword-arguments accepted by update_func must be provided to ScalarOperator via the kwarg accepted_kwargs as a tuple of Symbols. kwargs cannot be passed down to update_func if accepted_kwargs are not provided.

Warning

The user-provided update_func[!] must not use u in its computation. Positional argument (u, p, t) to update_func[!] are passed down by update_coefficients[!](L, u, p, t), where u is the input-vector to the composite AbstractSciMLOperator. For that reason, the values of u, or even shape, may not correspond to the input expected by update_func[!]. If an operator's state depends on its input vector, then it is, by definition, a nonlinear operator. We recommend sticking such nonlinearities in FunctionOperator. This topic is further discussed in (this issue)[https://github.com/SciML/SciMLOperators.jl/issues/159].

Interface

Lazy scalar algebra is defined for AbstractSciMLScalarOperators. The interface supports lazy addition, subtraction, multiplication and division.

Example

v = rand(4)
u = rand(4)
w = zeros(4)
p = nothing
t = 0.0

val_update = (a, u, p, t; scale = 0.0) -> scale
α = ScalarOperator(0.0; update_func = val_update, accepted_kwargs = (:scale,))
β = 2 * α + 3 / α

# Update β and evaluate with the new interface
result = β(v, u, p, t; scale = 1.0)

# In-place application
β(w, v, u, p, t; scale = 1.0)

# In-place with scaling
w_orig = copy(w)
α_val = 2.0
β_val = 0.5
β(w, v, u, p, t, α_val, β_val; scale = 1.0) # w = α_val*(β*v) + β_val*w
source
SciMLOperators.MatrixOperatorType

Represents a linear operator given by an AbstractMatrix that may be applied to an AbstractVecOrMat. Its state is updated by the user-provided update_func during operator evaluation (L([w,], v, u, p, t)), or by calls to update_coefficients[!](L, u, p, t). Both recursively call the update_function, update_func which is assumed to have the signature

update_func(A::AbstractMatrix, u, p, t; <accepted kwargs>) -> newA

or

update_func!(A::AbstractMatrix, u, p, t; <accepted kwargs>) -> [modifies A]

The set of keyword-arguments accepted by update_func[!] must be provided to MatrixOperator via the kwarg accepted_kwargs as a tuple of Symbols. kwargs cannot be passed down to update_func[!] if accepted_kwargs are not provided.

Warning

The user-provided update_func[!] must not use u in its computation. Positional argument (u, p, t) to update_func[!] are passed down by update_coefficients[!](L, u, p, t), where u is the input-vector to the composite AbstractSciMLOperator. For that reason, the values of u, or even shape, may not correspond to the input expected by update_func[!]. If an operator's state depends on its input vector, then it is, by definition, a nonlinear operator. We recommend sticking such nonlinearities in FunctionOperator. This topic is further discussed in (this issue)[https://github.com/SciML/SciMLOperators.jl/issues/159].

Interface

Lazy matrix algebra is defined for AbstractSciMLOperators. The Interface supports lazy addition, subtraction, multiplication, inversion, adjoints, transposes.

Example

Out-of-place update and usage

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

mat_update = (A, u, p, t; scale = 0.0) -> t * p
M = MatrixOperator(0.0; update_func = mat_update, accepted_kwargs = (:scale,))

L = M * M + 3I
L = cache_operator(L, v)

# update and evaluate 
w = L(v, u, p, t; scale = 1.0)

# In-place evaluation
w = similar(v)
L(w, v, u, p, t; scale = 1.0)

# In-place with scaling
β = 0.5
L(w, v, u, p, t, 2.0, β; scale = 1.0) # w = 2.0*(L*v) + 0.5*w

In-place update and usage

w = zeros(4)
v = zeros(4)
u = rand(4)
p = rand(4) # Must be non-nothing
t = rand()

mat_update! = (A, u, p, t; scale = 0.0) -> (A .= t * p * u' * scale)
M = MatrixOperator(zeros(4, 4); update_func! = mat_update!, accepted_kwargs = (:scale,))
L = M * M + 3I
L = cache_operator(L, v) 

# update L in-place and evaluate
update_coefficients!(L, u, p, t; scale = 1.0)
mul!(w, L, v)

# Or use the new interface that separates update and application
L(w, v, u, p, t; scale = 1.0)
source
SciMLOperators.DiagonalOperatorFunction
DiagonalOperator(
    diag;
    update_func,
    update_func!,
    accepted_kwargs
)

Represents an elementwise scaling (diagonal-scaling) operation that may be applied to an AbstractVecOrMat. When diag is an AbstractVector of length N, L = DiagonalOperator(diag, ...) can be applied to AbstractArrays with size(u, 1) == N. Each column of the v will be scaled by diag, as in LinearAlgebra.Diagonal(diag) * v.

When diag is a multidimensional array, L = DiagonalOperator(diag, ...) forms an operator of size (N, N) where N = size(diag, 1) is the leading length of diag. L then is the elementwise-scaling operation on arrays of length(v) = length(diag) with leading length size(u, 1) = N.

Its state is updated by the user-provided update_func during operator evaluation (L([w,], v, u, p, t)), or by calls to update_coefficients[!](L, u, p, t). Both recursively call the update_function, update_func which is assumed to have the signature

update_func(diag::AbstractVecOrMat, u, p, t; <accepted kwargs>) -> new_diag

or

update_func!(diag::AbstractVecOrMat, u, p, t; <accepted kwargs>) -> [modifies diag]

The set of keyword-arguments accepted by update_func[!] must be provided to MatrixOperator via the kwarg accepted_kwargs as a tuple of Symbols. kwargs cannot be passed down to update_func[!] if accepted_kwargs are not provided.

Warning

The user-provided update_func[!] must not use u in its computation. Positional argument (u, p, t) to update_func[!] are passed down by update_coefficients[!](L, u, p, t), where u is the input-vector to the composite AbstractSciMLOperator. For that reason, the values of u, or even shape, may not correspond to the input expected by update_func[!]. If an operator's state depends on its input vector, then it is, by definition, a nonlinear operator. We recommend sticking such nonlinearities in FunctionOperator. This topic is further discussed in (this issue)[https://github.com/SciML/SciMLOperators.jl/issues/159].

Example

source
SciMLOperators.AffineOperatorType

Represents a generalized affine operation (w = A * v + B * b) that may be applied to an AbstractVecOrMat. The user-provided update functions, update_func[!] update the AbstractVecOrMat b, and are called during operator evaluation (L([w,], v, u, p, t)), or by calls to update_coefficients[!](L, u, p, t). The update functions are assumed to have the syntax

update_func(b::AbstractVecOrMat, u, p, t; <accepted kwargs>) -> new_b

or

update_func!(b::AbstractVecOrMat, u ,p , t; <accepted kwargs>) -> [modifies b]

and B, b are expected to have an appropriate size so that A * v + B * b makes sense. Specifically, size(A, 1) == size(B, 1), and size(v, 2) == size(b, 2).

The set of keyword-arguments accepted by update_func[!] must be provided to AffineOperator via the kwarg accepted_kwargs as a tuple of Symbols. kwargs cannot be passed down to update_func[!] if accepted_kwargs are not provided.

Example

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

A = MatrixOperator(rand(4, 4))
B = MatrixOperator(rand(4, 4))

vec_update_func = (b, u, p, t) -> p .* u * t
L = AffineOperator(A, B, zeros(4); update_func = vec_update_func)
L = cache_operator(M, v)

# update L and evaluate
w = L(v, u, p, t) # == A * v + B * (p .* u * t)
source
SciMLOperators.AddVectorFunction
AddVector(b; update_func, update_func!, accepted_kwargs)

Represents the affine operation w = I * v + I * b. The update functions, update_func[!] update the state of AbstractVecOrMat b. See documentation of AffineOperator for more details.

source
AddVector(B, b; update_func, update_func!, accepted_kwargs)

Represents the affine operation w = I * v + B * b. The update functions, update_func[!] update the state of AbstractVecOrMat b. See documentation of AffineOperator for more details.

source
SciMLOperators.FunctionOperatorType

Matrix free operator given by a function

  • op: Function with signature op(v, u, p, t) and (if isinplace) op(w, v, u, p, t)

  • op_adjoint: Adjoint operator

  • op_inverse: Inverse operator

  • op_adjoint_inverse: Adjoint inverse operator

  • traits: Traits

  • u: State

  • p: Parameters

  • t: Time

  • cache: Cache

source
SciMLOperators.TensorProductOperatorType

Computes the lazy pairwise Kronecker product, or tensor product, operator of AbstractMatrix, and AbstractSciMLOperator subtypes. Calling ⊗(ops...) is equivalent to Base.kron(ops...). Fast operator evaluation is performed without forming the full tensor product operator.

TensorProductOperator(A, B) = A ⊗ B
TensorProductOperator(A, B, C) = A ⊗ B ⊗ C

(A ⊗ B)(v) = vec(B * reshape(v, M, N) * transpose(A))

where M = size(B, 2), and N = size(A, 2)

Example

using SciMLOperators, LinearAlgebra

# Create basic operators
A = rand(3, 3)
B = rand(4, 4)
A_op = MatrixOperator(A)
B_op = MatrixOperator(B)

# Create tensor product operator
T = A_op ⊗ B_op

# Apply to a vector using the new interface
v = rand(3*4)    # Action vector
u = rand(3*4)    # Update vector
p = nothing
t = 0.0

# Out-of-place application
result = T(v, u, p, t)

# For in-place operations, need to cache the operator first
T_cached = cache_operator(T, v)

# In-place application
w = zeros(size(T, 1))
T_cached(w, v, u, p, t)

# In-place with scaling
w_orig = copy(w)
α = 2.0
β = 0.5
T_cached(w, v, u, p, t, α, β) # w = α*(T*v) + β*w_orig
source

Lazy Scalar Operator Combination

Lazy Operator Combination