1
- # SciMLOperators.jl: Unified operator interface for ` SciML.ai ` and beyond
1
+ # SciMLOperators.jl: Unified operator interface for Julia and SciML
2
2
3
3
` SciMLOperators ` is a package for managing linear, nonlinear,
4
4
time-dependent, and parameter dependent operators acting on vectors,
@@ -25,69 +25,51 @@ using Pkg
25
25
Pkg. add (" SciMLOperators" )
26
26
```
27
27
28
- ## Examples
29
-
30
- Let ` M ` , ` D ` , ` F ` be matrix-based, diagonal-matrix-based, and function-based
31
- ` SciMLOperators ` respectively.
32
-
33
- ``` julia
34
- N = 4
35
- f = (v, u, p, t) -> u .* v
36
-
37
- M = MatrixOperator (rand (N, N))
38
- D = DiagonalOperator (rand (N))
39
- F = FunctionOperator (f, zeros (N), zeros (N))
40
- ```
41
-
42
- Then, the following codes just work.
43
-
44
- ``` julia
45
- L1 = 2 M + 3 F + LinearAlgebra. I + rand (N, N)
46
- L2 = D * F * M'
47
- L3 = kron (M, D, F)
48
- L4 = M \ D
49
- L5 = [M; D]' * [M F; F D] * [F; D]
50
- ```
51
-
52
- Each ` L# ` can be applied to ` AbstractVector ` s of appropriate sizes:
53
-
54
- ``` julia
55
- p = nothing # parameter struct
56
- t = 0.0 # time
57
-
58
- u = rand (N)
59
- v = rand (N)
60
- w = L1 (v, u, p, t) # == L1 * v
61
-
62
- v_kron = rand (N^ 3 )
63
- w_kron = L3 (v_kron, u, p, t) # == L3 * v_kron
64
- ```
65
-
66
- For mutating operator evaluations, call ` cache_operator ` to generate an
67
- in-place cache, so the operation is nonallocating.
68
-
69
- ``` julia
70
- α, β = rand (2 )
71
-
72
- # allocate cache
73
- L2 = cache_operator (L2, u)
74
- L4 = cache_operator (L4, u)
75
-
76
- # allocation-free evaluation
77
- L2 (w, v, u, p, t) # == mul!(w, L2, v)
78
- L4 (w, v, u, p, t, α, β) # == mul!(w, L4, v, α, β)
79
- ```
80
-
81
- The calling signature ` L(v, u, p, t) ` , for out-of-place evaluations, is
82
- equivalent to ` L * v ` , and the in-place evaluation ` L(w, v, u, p, t, args...) `
83
- is equivalent to ` LinearAlgebra.mul!(w, L, v, args...) ` , where the arguments
84
- ` u, p, t ` are passed to ` L ` to update its state. More details are provided
85
- in the operator update section below.
86
-
87
- The ` (v, u, p, t) ` calling signature is standardized over the ` SciML `
88
- ecosystem and is flexible enough to support use cases such as time-evolution
89
- in ODEs, as well as sensitivity computation with respect to the parameter
90
- object ` p ` .
28
+ ## Why ` SciMLOperators ` ?
29
+
30
+ Many functions, from linear solvers to differential equations, require
31
+ the use of matrix-free operators to achieve maximum performance in
32
+ many scenarios. ` SciMLOperators.jl ` defines the abstract interface for how
33
+ operators in the SciML ecosystem are supposed to be defined. It gives the
34
+ common set of functions and traits that solvers can rely on for properly
35
+ performing their tasks. Along with that, ` SciMLOperators.jl ` provides
36
+ definitions for the basic standard operators that are used as building
37
+ blocks for most tasks, simplifying the use of operators while also
38
+ demonstrating to users how such operators can be built and used in practice.
39
+
40
+ ` SciMLOperators.jl ` has the design that is required to be used in
41
+ all scenarios of equation solvers. For example, Magnus integrators for
42
+ differential equations require defining an operator `` u' = A(t) u `` , while
43
+ Munthe-Kaas methods require defining operators of the form `` u' = A(u) u `` .
44
+ Thus, the operators need some form of time and state dependence, which the
45
+ solvers can update and query when they are non-constant
46
+ (` update_coefficients! ` ). Additionally, the operators need the ability to
47
+ act like “normal” functions for equation solvers. For example, if ` A(v,u,p,t) `
48
+ has the same operation as ` update_coefficients(A, u, p, t); A * v ` , then ` A `
49
+ can be used in any place where a differential equation definition
50
+ ` (u,p,t) -> A(u, u, p, t) ` is used without requiring the user or solver to do any extra
51
+ work.
52
+
53
+ Another example is state-dependent mass matrices. ` M(u,p,t)*u' = f(u,p,t) ` .
54
+ When solving such an equation, the solver must understand how to "update M"
55
+ during operations, and thus the ability to update the state of ` M ` is a required
56
+ function in the interface. This is also required for the definition of Jacobians
57
+ ` J(u,p,t) ` in order to be properly used with Krylov methods inside of ODE solves
58
+ without reconstructing the matrix-free operator at each step.
59
+
60
+ Thus while previous good efforts for matrix-free operators have existed
61
+ in the Julia ecosystem, such as
62
+ [ LinearMaps.jl] ( https://github.com/JuliaLinearAlgebra/LinearMaps.jl ) , those
63
+ operator interfaces lack these aspects to actually be fully seamless
64
+ with downstream equation solvers. This necessitates the definition and use of
65
+ an extended operator interface with all of these properties, hence the
66
+ ` AbstractSciMLOperator ` interface.
67
+
68
+ !!! warn
69
+
70
+ This means that LinearMaps.jl is fundamentally lacking and is incompatible
71
+ with many of the tools in the SciML ecosystem, except for the specific cases
72
+ where the matrix-free operator is a constant!
91
73
92
74
## Features
93
75
0 commit comments