@@ -37,7 +37,7 @@ julia> Pkg.add("SciMLOperators")
3737
3838## Examples
3939
40- Let ` M ` , ` D ` , ` F ` be matrix, diagonal matrix, and function-based
40+ Let ` M ` , ` D ` , ` F ` be matrix-based , diagonal- matrix-based , and function-based
4141` SciMLOperators ` respectively.
4242
4343``` julia
@@ -87,220 +87,6 @@ L2(v, u, p, t) # == mul!(v, L2, u)
8787L4 (v, u, p, t, α, β) # == mul!(v, L4, u, α, β)
8888```
8989
90- The calling signature ` L(u, p, t) ` , for out-of-place evaluations is
91- equivalent to ` L * u ` , and the in-place evaluation ` L(v, u, p, t, args...) `
92- is equivalent to ` LinearAlgebra.mul!(v, L, u, args...) ` , where the arguments
93- ` p, t ` are passed to ` L ` to update its state. More details are provided
94- in the operator update section below. While overloads to ` Base.* `
95- and ` LinearAlgebra.mul! ` are available, where a ` SciMLOperator ` behaves
96- like an ` AbstractMatrix ` , we recommend sticking with the
97- ` L(u, p, t) ` , ` L(v, u, p, t) ` , ` L(v, u, p, t, α, β) ` calling signatures
98- as the latter internally update the operator state.
99-
100- The ` (u, p, t) ` calling signature is standardized over the ` SciML `
101- ecosystem and is flexible enough to support use cases such as time-evolution
102- in ODEs, as well as sensitivity computation with respect to the parameter
103- object ` p ` .
104-
105- Thanks to overloads defined for evaluation methods and traits in
106- ` Base ` , ` LinearAlgebra ` , the behaviour of a ` SciMLOperator ` is
107- indistinguishable from an ` AbstractMatrix ` . These operators can be
108- passed to linear solver packages, and even to ordinary differential
109- equation solvers. The list of overloads to the ` AbstractMatrix `
110- interface include, but are not limited, the following:
111-
112- - ` Base: size, zero, one, +, -, *, /, \, ∘, inv, adjoint, transpose, convert `
113- - ` LinearAlgebra: mul!, ldiv!, lmul!, rmul!, factorize, issymmetric, ishermitian, isposdef `
114- - ` SparseArrays: sparse, issparse `
115-
116- ## Multidimension arrays and batching
117-
118- SciMLOperator can also be applied to ` AbstractMatrix ` subtypes where
119- operator-evaluation is done column-wise.
120-
121- ``` julia
122- K = 10
123- u_mat = rand (N, K)
124-
125- v_mat = F (u_mat, p, t) # == mul!(v_mat, F, u_mat)
126- size (v_mat) == (N, K) # true
127- ```
128-
129- ` L# ` can also be applied to ` AbstractArray ` s that are not
130- ` AbstractVecOrMat ` s so long as their size in the first dimension is appropriate
131- for matrix-multiplication. Internally, ` SciMLOperator ` s reshapes an
132- ` N ` -dimensional array to an ` AbstractMatrix ` , and applies the operator via
133- matrix-multiplication.
134-
135- ## Operator update
136-
137- This package can also be used to write time-dependent, and
138- parameter-dependent operators, whose state can be updated per
139- a user-defined function.
140- The updates can be done in-place, i.e. by mutating the object,
141- or out-of-place, i.e. in a non-mutating, ` Zygote ` -compatible way.
142-
143- For example,
144-
145- ``` julia
146- u = rand (N)
147- p = rand (N)
148- t = rand ()
149-
150- # out-of-place update
151- mat_update_func = (A, u, p, t) -> t * (p * p' )
152- sca_update_func = (a, u, p, t) -> t * sum (p)
153-
154- M = MatrixOperator (zero (N, N); update_func = mat_update_func)
155- α = ScalarOperator (zero (Float64); update_func = sca_update_func)
156-
157- L = α * M
158- L = cache_operator (L, u)
159-
160- # L is initialized with zero state
161- L * u == zeros (N) # true
162-
163- # update operator state with `(u, p, t)`
164- L = update_coefficients (L, u, p, t)
165- # and multiply
166- L * u != zeros (N) # true
167-
168- # updates state and evaluates L at (u, p, t)
169- L (u, p, t) != zeros (N) # true
170- ```
171-
172- The out-of-place evaluation function ` L(u, p, t) ` calls
173- ` update_coefficients ` under the hood, which recursively calls
174- the ` update_func ` for each component ` SciMLOperator ` .
175- Therefore the out-of-place evaluation function is equivalent to
176- calling ` update_coefficients ` followed by ` Base.* ` . Notice that
177- the out-of-place evaluation does not return the updated operator.
178-
179- On the other hand,, the in-place evaluation function, ` L(v, u, p, t) ` ,
180- mutates ` L ` , and is equivalent to calling ` update_coefficients! `
181- followed by ` mul! ` . The in-place update behaviour works the same way
182- with a few ` <!> ` s appended here and there. For example,
183-
184- ``` julia
185- v = rand (N)
186- u = rand (N)
187- p = rand (N)
188- t = rand ()
189-
190- # in-place update
191- _A = rand (N, N)
192- _d = rand (N)
193- mat_update_func! = (A, u, p, t) -> (copy! (A, _A); lmul! (t, A); nothing )
194- diag_update_func! = (diag, u, p, t) -> copy! (diag, N)
195-
196- M = MatrixOperator (zero (N, N); update_func! = mat_update_func!)
197- D = DiagonalOperator (zero (N); update_func! = diag_update_func!)
198-
199- L = D * M
200- L = cache_operator (L, u)
201-
202- # L is initialized with zero state
203- L * u == zeros (N) # true
204-
205- # update L in-place
206- update_coefficients! (L, u, p, t)
207- # and multiply
208- mul! (v, u, p, t) != zero (N) # true
209-
210- # updates L in-place, and evaluates at (u, p, t)
211- L (v, u, p, t) != zero (N) # true
212- ```
213-
214- The update behaviour makes this package flexible enough to be used
215- in ` OrdianryDiffEq ` . As the parameter object ` p ` is often reserved
216- for sensitivy computation via automatic-differentiation, a user may
217- prefer to pass in state information via other arguments. For that
218- reason, we allow for update functions with arbitrary keyword arguments.
219-
220- ``` julia
221- mat_update_func = (A, u, p, t; scale = 0.0 ) -> scale * (p * p' )
222-
223- M = MatrixOperator (zero (N, N); update_func = mat_update_func,
224- accepted_kwargs = (:state ,))
225-
226- M (u, p, t) == zeros (N) # true
227- M (u, p, t; scale = 1.0 ) != zero (N)
228- ```
229-
230- ## Features
231-
232- * Matrix-free operators with ` FunctionOperator `
233- * Fast tensor product evaluation
234- * Mutating, nonmutating update behaviour (Zygote compatible)
235- * ` InvertibleOperator ` - pair fwd, bwd operators
236- * Lazy algebra: addition, subtraction, multiplication, inverse, adjoint
237- * Pre-caching methods for in-place evaluations
238-
239- ## Why ` SciMLOperators ` ?
240-
241- Many functions, from linear solvers to differential equations, require
242- the use of matrix-free operators in order to achieve maximum performance in
243- many scenarios. ` SciMLOperators.jl ` defines the abstract interface for how
244- operators in the SciML ecosystem are supposed to be defined. It gives the
245- common set of functions and traits which solvers can rely on for properly
246- performing their tasks. Along with that, ` SciMLOperators.jl ` provides
247- definitions for the basic standard operators which are used in building
248- blocks for most tasks, both simplifying the use of operators while also
249- demonstrating to users how such operators can be built and used in practice.
250-
251- ` SciMLOperators.jl ` has the design that is required in order to be used in
252- all scenarios of equation solvers. For example, Magnus integrators for
253- differential equations require defining an operator `` u' = A(t) u `` , while
254- Munthe-Kaas methods require defining operators of the form `` u' = A(u) u `` .
255- Thus the operators need some form of time and state dependence which the
256- solvers can update and query when they are non-constant
257- (` update_coefficients! ` ). Additionally, the operators need the ability to
258- act like "normal" functions for equation solvers. For example, if ` A(u,p,t) `
259- has the same operation as ` update_coefficients(A, u, p, t); A * u ` , then ` A `
260- can be used in any place where a differential equation definition
261- ` f(u, p, t) ` is used without requring the user or solver to do any extra
262- work. Thus while previous good efforts for matrix-free operators have existed
263- in the Julia ecosystem, such as
264- [ LinearMaps.jl] ( https://github.com/JuliaLinearAlgebra/LinearMaps.jl ) , those
265- operator interfaces lack these aspects in order to actually be fully seamless
266- with downstream equation solvers. This necessitates the definition and use of
267- an extended operator interface with all of these properties, hence the
268- ` AbstractSciMLOperator ` interface.
269-
270- Some packages providing similar functionality are
271- * [ LinearMaps.jl] ( https://github.com/JuliaLinearAlgebra/LinearMaps.jl )
272- * [ ` DiffEqOperators.jl ` ] ( https://github.com/SciML/DiffEqOperators.jl/tree/master ) (deprecated)
273-
274- ## Interoperability and extended Julia ecosystem
275-
276- ` SciMLOperator.jl ` overloads the ` AbstractMatrix ` interface for
277- ` AbstractSciMLOperator ` s, allowing seamless compatibility with
278- linear solves, and nonlinear solvers. Further, due to the update functionality,
279- ` AbstractSciMLOperator ` s can represent an ` ODEFunction ` in ` OrdinaryDiffEq.jl ` ,
280- and downstream packages. See tutorials for example of usage with
281- ` OrdinaryDiffEq.jl ` , ` LinearSolve.jl ` , ` NonlinearSolve.jl ` .
282-
283- Further, the nonmutating update functionality allows gradient propogation
284- through ` AbstractSciMLOperator ` s, and is compatible with
285- automatic-differentiation libraries like
286- [ ` Zygote.jl ` ] ( https://github.com/SciML/DiffEqOperators.jl/tree/master ) .
287- An example of ` Zygote.jl ` usage with
288- [ ` Lux.jl ` ] ( https://github.com/LuxDL/Lux.jl ) is also provided in the tutorials.
289-
290- Please make an issue [ here] ( https://github.com/SciML/SciMLOperators.jl/issues )
291- if you come across an unexpected issue while using ` SciMLOperators ` .
292-
293- We provide below a list of packages that make use of ` SciMLOperators ` .
294- If you are using ` SciMLOperators ` in your work, feel free to create a PR
295- and add your package to this list.
296-
297- * [ ` SciML.ai ` ] ( https://sciml.ai/ ) ecosystem: ` SciMLOperators ` is compatible with, and utilized by every ` SciML ` package.
298- * [ ` CalculustJL ` ] ( https://github.com/CalculustJL ) packages use ` SciMLOperators ` to define matrix-free vector-calculus operators for solving partial differential equations.
299- * [ ` CalculustCore.jl ` ] ( https://github.com/CalculustJL/CalculustCore.jl )
300- * [ ` FourierSpaces.jl ` ] ( https://github.com/CalculustJL/FourierSpaces.jl )
301- * [ ` NodalPolynomialSpaces.jl ` ] ( https://github.com/CalculustJL/NodalPolynomialSpaces.jl )
302- * ` SparseDiffTools.jl `
303-
30490## Roadmap
30591
30692- [ ] [ Complete integration with ` SciML ` ecosystem] ( https://github.com/SciML/SciMLOperators.jl/issues/142 )
@@ -309,15 +95,3 @@ and add your package to this list.
30995- [ ] [ Fast tensor-sum (` kronsum ` ) evaluation] ( https://github.com/SciML/SciMLOperators.jl/issues/53 )
31096- [ ] [ Fully flesh out operator array algebra] ( https://github.com/SciML/SciMLOperators.jl/issues/62 )
31197- [ ] [ Operator fusion/matrix chain multiplication at constant ` (u, p, t) ` -slices] ( https://github.com/SciML/SciMLOperators.jl/issues/51 )
312-
313- ## Contributing
314-
315- - Please refer to the
316- [ SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages] ( https://github.com/SciML/ColPrac/blob/master/README.md )
317- for guidance on PRs, issues, and other matters relating to contributing to SciML.
318- - There are a few community forums:
319- - The #diffeq-bridged and #sciml-bridged channels in the
320- [ Julia Slack] ( https://julialang.org/slack/ )
321- - [ JuliaDiffEq] ( https://gitter.im/JuliaDiffEq/Lobby ) on Gitter
322- - On the Julia Discourse forums (look for the [ modelingtoolkit tag] ( https://discourse.julialang.org/tag/modelingtoolkit )
323- - See also [ SciML Community page] ( https://sciml.ai/community/ )
0 commit comments