From 165965e0d094b17dd800013a7734d4bd25c8ae6b Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Sat, 28 Sep 2024 12:43:54 +0000 Subject: [PATCH] build based on 6d4ca4a --- dev/.documenter-siteinfo.json | 2 +- dev/benchmarks/hodgkinhuxley/index.html | 2 +- dev/benchmarks/lotkavolterra/index.html | 2 +- .../multi-language-wrappers/index.html | 2 +- dev/benchmarks/orego/index.html | 2 +- dev/benchmarks/pleiades/index.html | 2 +- dev/benchmarks/rober/index.html | 2 +- dev/benchmarks/vanderpol/index.html | 2 +- dev/diffusions/index.html | 2 +- dev/filtering/index.html | 14 +- dev/implementation/index.html | 2 +- dev/index.html | 2 +- dev/initialization/index.html | 2 +- dev/likelihoods/index.html | 4 +- dev/objects.inv | Bin 3055 -> 3055 bytes dev/priors/49473a66.svg | 174 ----------------- dev/priors/53661c5c.svg | 174 +++++++++++++++++ dev/priors/69865de8.svg | 174 ----------------- dev/priors/78c44fc4.svg | 106 ++++++++++ dev/priors/8d30ced6.svg | 174 ----------------- dev/priors/a292ceb1.svg | 104 ---------- dev/priors/a350c826.svg | 174 +++++++++++++++++ dev/priors/d41d5e0e.svg | 174 +++++++++++++++++ dev/priors/d84323e8.svg | 178 ----------------- dev/priors/e224c6c6.svg | 182 ++++++++++++++++++ dev/priors/index.html | 20 +- dev/references/index.html | 2 +- dev/solvers/index.html | 8 +- dev/tutorials/dae/3542b1a5.svg | 90 +++++++++ dev/tutorials/dae/b933a08f.svg | 59 ++++++ dev/tutorials/dae/bbfc7e14.svg | 59 ------ dev/tutorials/dae/eec28457.svg | 90 --------- dev/tutorials/dae/index.html | 86 ++++----- .../{5a45fb09.svg => 2353f20a.svg} | 2 +- dev/tutorials/dynamical_odes/239c19ed.svg | 52 ----- dev/tutorials/dynamical_odes/5884f0f2.svg | 50 +++++ .../{08393b82.svg => 64959796.svg} | 2 +- dev/tutorials/dynamical_odes/afc0b3df.svg | 50 ----- dev/tutorials/dynamical_odes/ecd44a8c.svg | 52 +++++ dev/tutorials/dynamical_odes/index.html | 8 +- .../exponential_integrators/index.html | 2 +- dev/tutorials/getting_started/7dc5a135.svg | 50 ----- dev/tutorials/getting_started/b2752536.svg | 50 +++++ dev/tutorials/getting_started/b98fcdaf.svg | 50 +++++ dev/tutorials/getting_started/bd9f291c.svg | 54 ------ dev/tutorials/getting_started/c8d5dd68.svg | 54 ++++++ dev/tutorials/getting_started/d81735df.svg | 50 ----- dev/tutorials/getting_started/index.html | 118 ++++++------ .../ode_parameter_inference/48a1b67f.svg | 87 --------- .../ode_parameter_inference/689abaa5.svg | 92 --------- .../ode_parameter_inference/6f790629.svg | 87 +++++++++ .../ode_parameter_inference/74285845.svg | 92 +++++++++ .../ode_parameter_inference/9a3cb630.svg | 92 +++++++++ .../ode_parameter_inference/ea82dc47.svg | 92 --------- .../ode_parameter_inference/index.html | 14 +- 55 files changed, 1638 insertions(+), 1632 deletions(-) delete mode 100644 dev/priors/49473a66.svg create mode 100644 dev/priors/53661c5c.svg delete mode 100644 dev/priors/69865de8.svg create mode 100644 dev/priors/78c44fc4.svg delete mode 100644 dev/priors/8d30ced6.svg delete mode 100644 dev/priors/a292ceb1.svg create mode 100644 dev/priors/a350c826.svg create mode 100644 dev/priors/d41d5e0e.svg delete mode 100644 dev/priors/d84323e8.svg create mode 100644 dev/priors/e224c6c6.svg create mode 100644 dev/tutorials/dae/3542b1a5.svg create mode 100644 dev/tutorials/dae/b933a08f.svg delete mode 100644 dev/tutorials/dae/bbfc7e14.svg delete mode 100644 dev/tutorials/dae/eec28457.svg rename dev/tutorials/dynamical_odes/{5a45fb09.svg => 2353f20a.svg} (99%) delete mode 100644 dev/tutorials/dynamical_odes/239c19ed.svg create mode 100644 dev/tutorials/dynamical_odes/5884f0f2.svg rename dev/tutorials/dynamical_odes/{08393b82.svg => 64959796.svg} (86%) delete mode 100644 dev/tutorials/dynamical_odes/afc0b3df.svg create mode 100644 dev/tutorials/dynamical_odes/ecd44a8c.svg delete mode 100644 dev/tutorials/getting_started/7dc5a135.svg create mode 100644 dev/tutorials/getting_started/b2752536.svg create mode 100644 dev/tutorials/getting_started/b98fcdaf.svg delete mode 100644 dev/tutorials/getting_started/bd9f291c.svg create mode 100644 dev/tutorials/getting_started/c8d5dd68.svg delete mode 100644 dev/tutorials/getting_started/d81735df.svg delete mode 100644 dev/tutorials/ode_parameter_inference/48a1b67f.svg delete mode 100644 dev/tutorials/ode_parameter_inference/689abaa5.svg create mode 100644 dev/tutorials/ode_parameter_inference/6f790629.svg create mode 100644 dev/tutorials/ode_parameter_inference/74285845.svg create mode 100644 dev/tutorials/ode_parameter_inference/9a3cb630.svg delete mode 100644 dev/tutorials/ode_parameter_inference/ea82dc47.svg diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index fae51f8eb..27c029bf9 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-09-24T15:07:24","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-09-28T12:43:44","documenter_version":"1.7.0"}} \ No newline at end of file diff --git a/dev/benchmarks/hodgkinhuxley/index.html b/dev/benchmarks/hodgkinhuxley/index.html index ab64b0260..a410bc06b 100644 --- a/dev/benchmarks/hodgkinhuxley/index.html +++ b/dev/benchmarks/hodgkinhuxley/index.html @@ -653,4 +653,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/benchmarks/lotkavolterra/index.html b/dev/benchmarks/lotkavolterra/index.html index 7304bf27c..2c7686d83 100644 --- a/dev/benchmarks/lotkavolterra/index.html +++ b/dev/benchmarks/lotkavolterra/index.html @@ -748,4 +748,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/benchmarks/multi-language-wrappers/index.html b/dev/benchmarks/multi-language-wrappers/index.html index a09f608b7..7c8b0c87d 100644 --- a/dev/benchmarks/multi-language-wrappers/index.html +++ b/dev/benchmarks/multi-language-wrappers/index.html @@ -692,4 +692,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/benchmarks/orego/index.html b/dev/benchmarks/orego/index.html index c2bc9dda0..1a515d4e5 100644 --- a/dev/benchmarks/orego/index.html +++ b/dev/benchmarks/orego/index.html @@ -521,4 +521,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/benchmarks/pleiades/index.html b/dev/benchmarks/pleiades/index.html index 7809023fb..8879c00c4 100644 --- a/dev/benchmarks/pleiades/index.html +++ b/dev/benchmarks/pleiades/index.html @@ -564,4 +564,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/benchmarks/rober/index.html b/dev/benchmarks/rober/index.html index bb5cde8c7..a2a5a6361 100644 --- a/dev/benchmarks/rober/index.html +++ b/dev/benchmarks/rober/index.html @@ -517,4 +517,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/benchmarks/vanderpol/index.html b/dev/benchmarks/vanderpol/index.html index 783ec22a6..19aef22a0 100644 --- a/dev/benchmarks/vanderpol/index.html +++ b/dev/benchmarks/vanderpol/index.html @@ -658,4 +658,4 @@ [8e850b90] libblastrampoline_jll v5.8.0+1 [8e850ede] nghttp2_jll v1.52.0+1 [3f19e933] p7zip_jll v17.4.0+2 -Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` +Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m` diff --git a/dev/diffusions/index.html b/dev/diffusions/index.html index b5446781c..5ebfe28df 100644 --- a/dev/diffusions/index.html +++ b/dev/diffusions/index.html @@ -5,4 +5,4 @@ \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1, \\ \text{d} Y^{(q)}(t) &= \textcolor{#389826}{A} Y(t) \ \text{d}t + \textcolor{#4063D8}{\Gamma} \ \text{d}W(t), \\ Y(0) &\sim \textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }. -\end{aligned}\]

Then $Y^{(i)}(t)$ models the $i$-th derivative of $y(t)$. In this section, we consider choices relating to the "diffusion" $\textcolor{#4063D8}{\Gamma}$. If you're more interested in the drift matrix $\textcolor{#389826}{A}$ check out the Priors section, and for info on the initial distribution $\textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }$ check out the Initialization section.

Diffusion and calibration

We call $\textcolor{#4063D8}{\Gamma}$ the "diffusion" parameter. Since it is typically not known we need to estimate it; this is called "calibration".

ProbNumDiffEq.jl provides a few different choices for how to model and estimate $\textcolor{#4063D8}{\Gamma}$:

Or more compactly:

Isotropic:Diagonal (only for the EK0)
Time-varyingDynamicDiffusionDynamicMVDiffusion
Time-fixedFixedDiffusionFixedMVDiffusion

For more details on diffusions and calibration, check out this paper [9].

API

ProbNumDiffEq.DynamicDiffusion — Type
DynamicDiffusion()

Time-varying, isotropic diffusion, which is quasi-maximum-likelihood-estimated at each step.

This is the recommended diffusion when using adaptive step-size selection, and in particular also when solving stiff systems.

source
ProbNumDiffEq.FixedDiffusion — Type
FixedDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, isotropic diffusion, which is (optionally) quasi-maximum-likelihood-estimated.

This is the recommended diffusion when using fixed steps.

By default with calibrate=true, all covariances are re-scaled at the end of the solve with the MLE diffusion. Set calibrate=false to skip this step, e.g. when setting the initial_diffusion and then estimating the diffusion outside of the solver (e.g. with Fenrir.jl).

source
ProbNumDiffEq.DynamicMVDiffusion — Type
DynamicMVDiffusion()

Time-varying, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of DynamicDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [9] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source
ProbNumDiffEq.FixedMVDiffusion — Type
FixedMVDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of FixedDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [9] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source

References

[9]
N. Bosch, P. Hennig and F. Tronarp. Calibrated Adaptive Probabilistic ODE Solvers. In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, Vol. 130 of Proceedings of Machine Learning Research, edited by A. Banerjee and K. Fukumizu (PMLR, 13–15 Apr 2021); pp. 3466–3474.
+\end{aligned}\]

Then $Y^{(i)}(t)$ models the $i$-th derivative of $y(t)$. In this section, we consider choices relating to the "diffusion" $\textcolor{#4063D8}{\Gamma}$. If you're more interested in the drift matrix $\textcolor{#389826}{A}$ check out the Priors section, and for info on the initial distribution $\textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }$ check out the Initialization section.

Diffusion and calibration

We call $\textcolor{#4063D8}{\Gamma}$ the "diffusion" parameter. Since it is typically not known we need to estimate it; this is called "calibration".

ProbNumDiffEq.jl provides a few different choices for how to model and estimate $\textcolor{#4063D8}{\Gamma}$:

Or more compactly:

Isotropic:Diagonal (only for the EK0)
Time-varyingDynamicDiffusionDynamicMVDiffusion
Time-fixedFixedDiffusionFixedMVDiffusion

For more details on diffusions and calibration, check out this paper [9].

API

ProbNumDiffEq.DynamicDiffusion — Type
DynamicDiffusion()

Time-varying, isotropic diffusion, which is quasi-maximum-likelihood-estimated at each step.

This is the recommended diffusion when using adaptive step-size selection, and in particular also when solving stiff systems.

source
ProbNumDiffEq.FixedDiffusion — Type
FixedDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, isotropic diffusion, which is (optionally) quasi-maximum-likelihood-estimated.

This is the recommended diffusion when using fixed steps.

By default with calibrate=true, all covariances are re-scaled at the end of the solve with the MLE diffusion. Set calibrate=false to skip this step, e.g. when setting the initial_diffusion and then estimating the diffusion outside of the solver (e.g. with Fenrir.jl).

source
ProbNumDiffEq.DynamicMVDiffusion — Type
DynamicMVDiffusion()

Time-varying, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of DynamicDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [9] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source
ProbNumDiffEq.FixedMVDiffusion — Type
FixedMVDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of FixedDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [9] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source

References

[9]
N. Bosch, P. Hennig and F. Tronarp. Calibrated Adaptive Probabilistic ODE Solvers. In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, Vol. 130 of Proceedings of Machine Learning Research, edited by A. Banerjee and K. Fukumizu (PMLR, 13–15 Apr 2021); pp. 3466–3474.
diff --git a/dev/filtering/index.html b/dev/filtering/index.html index 420d3c25a..e0ec77a6d 100644 --- a/dev/filtering/index.html +++ b/dev/filtering/index.html @@ -1,26 +1,26 @@ -Filtering and Smoothing · ProbNumDiffEq.jl
GitHub

Gaussian Filtering and Smoothing

Predict

ProbNumDiffEq.predict — Function
predict(x::Gaussian, A::AbstractMatrix, Q::AbstractMatrix)

Prediction step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$, compute and return $\mathcal{N}(A μ, A Σ A^T + Q)$.

See also the non-allocating square-root version predict!.

source
ProbNumDiffEq.predict! — Function
predict!(x_out, x_curr, Ah, Qh, cachemat)

In-place and square-root implementation of predict which saves the result into x_out.

Only works with PSDMatrices.PSDMatrix types as Ah, Qh, and in the covariances of x_curr and x_out (both of type Gaussian). To prevent allocations, a cache matrix cachemat of size $D \times 2D$ (where $D \times D$ is the size of Ah and Qh) needs to be passed.

See also: predict.

source

Update

ProbNumDiffEq.update — Function
update(x, measurement, H)

Update step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$ and a measurement $z = \mathcal{N}(\hat{z}, S)$, with $S = H Σ H^T$, compute

\[\begin{aligned} +Filtering and Smoothing · ProbNumDiffEq.jl

GitHub

Gaussian Filtering and Smoothing

Predict

ProbNumDiffEq.predict — Function
predict(x::Gaussian, A::AbstractMatrix, Q::AbstractMatrix)

Prediction step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$, compute and return $\mathcal{N}(A μ, A Σ A^T + Q)$.

See also the non-allocating square-root version predict!.

source
ProbNumDiffEq.predict! — Function
predict!(x_out, x_curr, Ah, Qh, cachemat)

In-place and square-root implementation of predict which saves the result into x_out.

Only works with PSDMatrices.PSDMatrix types as Ah, Qh, and in the covariances of x_curr and x_out (both of type Gaussian). To prevent allocations, a cache matrix cachemat of size $D \times 2D$ (where $D \times D$ is the size of Ah and Qh) needs to be passed.

See also: predict.

source

Update

ProbNumDiffEq.update — Function
update(x, measurement, H)

Update step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$ and a measurement $z = \mathcal{N}(\hat{z}, S)$, with $S = H Σ H^T$, compute

\[\begin{aligned} K &= Σ^P H^T S^{-1}, \\ μ^F &= μ + K (0 - \hat{z}), \\ Σ^F &= Σ - K S K^T, -\end{aligned}\]

and return an updated state \mathcal{N}(μ^F, Σ^F). Note that this assumes zero-measurements. When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

For better performance, we recommend to use the non-allocating update!.

source
ProbNumDiffEq.update! — Function
update!(x_out, x_pred, measurement, H, K_cache, M_cache, S_cache)

In-place and square-root implementation of update which saves the result into x_out.

Implemented in Joseph Form to retain the PSDMatrix covariances:

\[\begin{aligned} +\end{aligned}\]

and return an updated state \mathcal{N}(μ^F, Σ^F). Note that this assumes zero-measurements. When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

For better performance, we recommend to use the non-allocating update!.

source
ProbNumDiffEq.update! — Function
update!(x_out, x_pred, measurement, H, K_cache, M_cache, S_cache)

In-place and square-root implementation of update which saves the result into x_out.

Implemented in Joseph Form to retain the PSDMatrix covariances:

\[\begin{aligned} K &= Σ^P H^T S^{-1}, \\ μ^F &= μ + K (0 - \hat{z}), \\ \sqrt{Σ}^F &= (I - KH) \sqrt(Σ), -\end{aligned}\]

where $\sqrt{M}$ denotes the left square-root of a matrix M, i.e. $M = \sqrt{M} \sqrt{M}^T$.

To prevent allocations, write into caches K_cache and M_cache, both of size D × D, and S_cache of same type as measurement.Σ.

See also: update.

source

Smooth

ProbNumDiffEq.smooth — Function
smooth(x_curr, x_next_smoothed, A, Q)

Update step of the Kalman smoother, aka. Rauch-Tung-Striebel smoother, for linear dynamics models.

Given Gaussians $x_n = \mathcal{N}(μ_{n}, Σ_{n})$ and $x_{n+1} = \mathcal{N}(μ_{n+1}^S, Σ_{n+1}^S)$, compute

\[\begin{aligned} +\end{aligned}\]

where $\sqrt{M}$ denotes the left square-root of a matrix M, i.e. $M = \sqrt{M} \sqrt{M}^T$.

To prevent allocations, write into caches K_cache and M_cache, both of size D × D, and S_cache of same type as measurement.Σ.

See also: update.

source

Smooth

ProbNumDiffEq.smooth — Function
smooth(x_curr, x_next_smoothed, A, Q)

Update step of the Kalman smoother, aka. Rauch-Tung-Striebel smoother, for linear dynamics models.

Given Gaussians $x_n = \mathcal{N}(μ_{n}, Σ_{n})$ and $x_{n+1} = \mathcal{N}(μ_{n+1}^S, Σ_{n+1}^S)$, compute

\[\begin{aligned} μ_{n+1}^P &= A μ_n^F, \\ P_{n+1}^P &= A Σ_n^F A + Q, \\ G &= Σ_n^S A^T (Σ_{n+1}^P)^{-1}, \\ μ_n^S &= μ_n^F + G (μ_{n+1}^S - μ_{n+1}^P), \\ Σ_n^S &= (I - G A) Σ_n^F (I - G A)^T + G Q G^T + G Σ_{n+1}^S G^T, -\end{aligned}\]

and return a smoothed state \mathcal{N}(μ_n^S, Σ_n^S). When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

source

Markov Kernels

ProbNumDiffEq.AffineNormalKernel — Type
AffineNormalKernel(A[, b], C)

Structure to represent affine Normal Markov kernels, i.e. conditional distributions of the form

\[\begin{aligned} +\end{aligned}\]

and return a smoothed state \mathcal{N}(μ_n^S, Σ_n^S). When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

source

Markov Kernels

ProbNumDiffEq.AffineNormalKernel — Type
AffineNormalKernel(A[, b], C)

Structure to represent affine Normal Markov kernels, i.e. conditional distributions of the form

\[\begin{aligned} y \mid x \sim \mathcal{N} \left( y; A x + b, C \right). -\end{aligned}\]

At the point of writing, AffineNormalKernels are only used to precompute and store the backward representation of the posterior (via compute_backward_kernel!) and for smoothing (via marginalize!).

source
ProbNumDiffEq.marginalize! — Function
marginalize!(
     xout::Gaussian{Vector{T},PSDMatrix{T,S}}
     x::Gaussian{Vector{T},PSDMatrix{T,S}},
     K::AffineNormalKernel{<:AbstractMatrix,Union{<:Number,<:AbstractVector,Missing},<:PSDMatrix};
     C_DxD, C_3DxD
-)

Basically the same as predict!), but in kernel language and with support for affine transitions. At the time of writing, this is only used to smooth the posterior using it's backward representation, where the kernels are precomputed with compute_backward_kernel!.

Note that this function assumes certain shapes:

  • size(x.μ) == (D, D)
  • size(x.Σ) == (D, D)
  • size(K.A) == (D, D)
  • size(K.b) == (D,), or missing
  • size(K.C) == (D, D), _but with a tall square-root size(K.C.R) == (3D, D)

xout is assumes to have the same shapes as x.

source
ProbNumDiffEq.compute_backward_kernel! — Function
compute_backward_kernel!(Kout, xpred, x, K; C_DxD[, diffusion=1])

Compute the backward representation of the posterior, i.e. the conditional distribution of the current state given the next state and the transition kernel.

More precisely, given a distribution (x)

\[\begin{aligned} +)

Basically the same as predict!), but in kernel language and with support for affine transitions. At the time of writing, this is only used to smooth the posterior using it's backward representation, where the kernels are precomputed with compute_backward_kernel!.

Note that this function assumes certain shapes:

  • size(x.μ) == (D, D)
  • size(x.Σ) == (D, D)
  • size(K.A) == (D, D)
  • size(K.b) == (D,), or missing
  • size(K.C) == (D, D), _but with a tall square-root size(K.C.R) == (3D, D)

xout is assumes to have the same shapes as x.

source
ProbNumDiffEq.compute_backward_kernel! — Function
compute_backward_kernel!(Kout, xpred, x, K; C_DxD[, diffusion=1])

Compute the backward representation of the posterior, i.e. the conditional distribution of the current state given the next state and the transition kernel.

More precisely, given a distribution (x)

\[\begin{aligned} x \sim \mathcal{N} \left( x; μ, Σ \right), \end{aligned}\]

a kernel (K)

\[\begin{aligned} y \mid x \sim \mathcal{N} \left( y; A x + b, C \right), @@ -34,4 +34,4 @@ G &= Σ A^\top (Σ^P)^{-1}, \\ d &= μ - G μ^P, \\ Λ &= Σ - G Σ^P G^\top. -\end{aligned}\]

Everything is computed in square-root form and with minimal allocations (thus the cache C_DxD), so the actual formulas implemented here differ a bit.

The resulting backward kernels are used to smooth the posterior, via marginalize!.

source
+\end{aligned}\]

Everything is computed in square-root form and with minimal allocations (thus the cache C_DxD), so the actual formulas implemented here differ a bit.

The resulting backward kernels are used to smooth the posterior, via marginalize!.

source
diff --git a/dev/implementation/index.html b/dev/implementation/index.html index c96fb30b3..8e20da6e5 100644 --- a/dev/implementation/index.html +++ b/dev/implementation/index.html @@ -1,2 +1,2 @@ -Implementation via OrdinaryDiffEq.jl · ProbNumDiffEq.jl
GitHub

Solver Implementation via OrdinaryDiffEq.jl

ProbNumDiffEq.jl builds directly on OrdinaryDiffEq.jl to benefit from its iterator interface, flexible step-size control, and efficient Jacobian calculations. But, this requires extending non-public APIs. This page is meant to provide an overview on which parts exactly ProbNumDiffEq.jl builds on.

For more discussion on the pros and cons of building on OrdinaryDiffEq.jl, see this thread on discourse.

Building on OrdinaryDiffEq.jl

ProbNumDiffEq.jl shares most of OrdinaryDiffEq.jl's implementation. In particular:

  1. OrdinaryDiffEq.__init builds the cache and the integrator, and calls OrdinaryDiffEq.initialize!
  2. OrdinaryDiffEq.solve! implements the actual iterator structure, with
    • OrdinaryDiffEq.loopheader!
    • OrdinaryDiffEq.perform_step!
    • OrdinaryDiffEq.loopfooter!
    • OrdinaryDiffEq.postamble!

ProbNumDiffEq.jl builds around this structure and overloads some of the parts:

  • Algorithms: EK0/EK1 <: AbstractEK <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm
    • ./src/algorithms.jl provides the algorithms themselves
    • ./src/alg_utils.jl implements many traits (relating to automatic differentiation, implicitness, step-size control, etc)
  • Cache: EKCache <: AbstractODEFilterCache <: OrdinaryDiffEq.OrdinaryDiffEqCache
    • ./src/caches.jl implements the cache and its main constructor: OrdinaryDiffEq.alg_cache
  • Initialization and perform_step!: via OrdinaryDiffEq.initialize! and OrdinaryDiffEq.perform_step!. Implemented in ./src/perform_step.jl.
  • Custom postamble by overloading OrdinaryDiffEq.postamble! (which should always call OrdinaryDiffEq._postamble!). This is where we do the "smoothing" of the solution. Implemented in ./src/integrator_utils.jl.
  • Custom saving by overloading OrdinaryDiffEq.savevalues! (which should always call OrdinaryDiffEq._savevalues!). Implemented in ./src/integrator_utils.jl.

Building on DiffEqBase.jl

  • DiffEqBase.__init is currently overloaded to transform OOP problems into IIP problems (in ./src/solve.jl).
  • The solution object: ProbODESolution <: AbstractProbODESolution <: DiffEqBase.AbstractODESolution
    • ./src/solution.jl implements the main parts. Note that the main constructor DiffEqBase.build_solution is called by OrdinaryDiffEq.__init, so OrdinaryDiffEq.jl has control over its inputs.
    • MeanProbODESolution <: DiffEqBase.AbstractODESolution is a wrapper that allows handling the mean of a probabilistic ODE solution the same way one would handle any "standard" ODE solution, by just ignoring the covariances.
    • AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation handles the interpolation.
    • Plot recipe in ./ext/RecipesBaseExt.jl
    • Sampling in ./src/solution_sampling.jl
  • DiffEqBase.prepare_alg(::EK1{0}); closely follows a similar function implemented in OrdinaryDiffEq.jl ./src/alg_utils.jl
    • this also required DiffEqBase.remake(::EK1)

Other packages

  • DiffEqDevTools.appxtrue: We extend this function to work with ProbODESolution. This also enables DiffEqDevTools.WorkPrecision to work out of the box.
+Implementation via OrdinaryDiffEq.jl · ProbNumDiffEq.jl
GitHub

Solver Implementation via OrdinaryDiffEq.jl

ProbNumDiffEq.jl builds directly on OrdinaryDiffEq.jl to benefit from its iterator interface, flexible step-size control, and efficient Jacobian calculations. But, this requires extending non-public APIs. This page is meant to provide an overview on which parts exactly ProbNumDiffEq.jl builds on.

For more discussion on the pros and cons of building on OrdinaryDiffEq.jl, see this thread on discourse.

Building on OrdinaryDiffEq.jl

ProbNumDiffEq.jl shares most of OrdinaryDiffEq.jl's implementation. In particular:

  1. OrdinaryDiffEq.__init builds the cache and the integrator, and calls OrdinaryDiffEq.initialize!
  2. OrdinaryDiffEq.solve! implements the actual iterator structure, with
    • OrdinaryDiffEq.loopheader!
    • OrdinaryDiffEq.perform_step!
    • OrdinaryDiffEq.loopfooter!
    • OrdinaryDiffEq.postamble!

ProbNumDiffEq.jl builds around this structure and overloads some of the parts:

  • Algorithms: EK0/EK1 <: AbstractEK <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm
    • ./src/algorithms.jl provides the algorithms themselves
    • ./src/alg_utils.jl implements many traits (relating to automatic differentiation, implicitness, step-size control, etc)
  • Cache: EKCache <: AbstractODEFilterCache <: OrdinaryDiffEq.OrdinaryDiffEqCache
    • ./src/caches.jl implements the cache and its main constructor: OrdinaryDiffEq.alg_cache
  • Initialization and perform_step!: via OrdinaryDiffEq.initialize! and OrdinaryDiffEq.perform_step!. Implemented in ./src/perform_step.jl.
  • Custom postamble by overloading OrdinaryDiffEq.postamble! (which should always call OrdinaryDiffEq._postamble!). This is where we do the "smoothing" of the solution. Implemented in ./src/integrator_utils.jl.
  • Custom saving by overloading OrdinaryDiffEq.savevalues! (which should always call OrdinaryDiffEq._savevalues!). Implemented in ./src/integrator_utils.jl.

Building on DiffEqBase.jl

  • DiffEqBase.__init is currently overloaded to transform OOP problems into IIP problems (in ./src/solve.jl).
  • The solution object: ProbODESolution <: AbstractProbODESolution <: DiffEqBase.AbstractODESolution
    • ./src/solution.jl implements the main parts. Note that the main constructor DiffEqBase.build_solution is called by OrdinaryDiffEq.__init, so OrdinaryDiffEq.jl has control over its inputs.
    • MeanProbODESolution <: DiffEqBase.AbstractODESolution is a wrapper that allows handling the mean of a probabilistic ODE solution the same way one would handle any "standard" ODE solution, by just ignoring the covariances.
    • AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation handles the interpolation.
    • Plot recipe in ./ext/RecipesBaseExt.jl
    • Sampling in ./src/solution_sampling.jl
  • DiffEqBase.prepare_alg(::EK1{0}); closely follows a similar function implemented in OrdinaryDiffEq.jl ./src/alg_utils.jl
    • this also required DiffEqBase.remake(::EK1)

Other packages

  • DiffEqDevTools.appxtrue: We extend this function to work with ProbODESolution. This also enables DiffEqDevTools.WorkPrecision to work out of the box.
diff --git a/dev/index.html b/dev/index.html index b0949a5bb..546f542ad 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ Home · ProbNumDiffEq.jl
GitHub

Probabilistic Numerical Differential Equation Solvers

Banner

ProbNumDiffEq.jl provides probabilistic numerical solvers to the DifferentialEquations.jl ecosystem. The implemented ODE filters solve differential equations via Bayesian filtering and smoothing and compute not just a single point estimate of the true solution, but a posterior distribution that contains an estimate of its numerical approximation error.

For a short intro video, check out our poster presentation at JuliaCon2021.

Installation

Run Julia, enter ] to bring up Julia's package manager, and add the ProbNumDiffEq.jl package:

julia> ]
-(v1.10) pkg> add ProbNumDiffEq

Getting Started

For a quick introduction check out the "Solving ODEs with Probabilistic Numerics" tutorial.

Features

  • ProbDiffEq is similar in scope to ProbNumDiffEq.jl and it provides fast and feature-rich probabilistic ODE solvers but is implemented in Python and built on JAX.
  • ProbNum implements a wide range of probabilistic numerical methods, not only for ODEs but also for linear algebra, quadrature, and filtering/smoothing. It is implemented in Python and NumPy, and it focuses more on breadth and didactic purposes than on performance.
+(v1.10) pkg> add ProbNumDiffEq

Getting Started

For a quick introduction check out the "Solving ODEs with Probabilistic Numerics" tutorial.

Features

diff --git a/dev/initialization/index.html b/dev/initialization/index.html index de80814df..547a1bed3 100644 --- a/dev/initialization/index.html +++ b/dev/initialization/index.html @@ -11,4 +11,4 @@ \end{aligned}\]

It is clear that this contains quite some information for $Y(0)$: The initial value $y_0$ and the vector field $f$ imply

\[\begin{aligned} Y^{(0)}(0) &= y_0, \\ Y^{(1)}(0) &= f(y_0, 0). -\end{aligned}\]

It turns out that we can also compute higher-order derivatives of $y$ with the chain rule, and then use these to better initialize $Y^{(i)}(0)$. This, done efficiently with Taylor-mode automatic differentiation by using TaylorIntegration.jl, is what TaylorModeInit does. See also [1].

In the vast majority of cases, just stick to the exact Taylor-mode initialization TaylorModeInit!

API

ProbNumDiffEq.TaylorModeInit — Type
TaylorModeInit(order)

Exact initialization via Taylor-mode automatic differentiation up to order order.

This is the recommended initialization method!

It uses TaylorIntegration.jl to efficiently compute the higher-order derivatives of the solution at the initial value, via Taylor-mode automatic differentiation.

In some special cases it can happen that TaylorIntegration.jl is incompatible with the given problem (typically because the problem definition does not allow for elements of type Taylor). If this happens, try one of SimpleInit, ForwardDiffInit (for low enough orders), ClassicSolverInit.

References

  • [7] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
source
ProbNumDiffEq.ForwardDiffInit — Type
ForwardDiffInit(order)

Exact initialization via ForwardDiff.jl up to order order.

Warning: This does not scale well to high orders! For orders > 3, TaylorModeInit most likely performs better.

source
ProbNumDiffEq.SimpleInit — Type
SimpleInit()

Simple initialization, only with the given initial value and derivative.

The remaining derivatives are set to zero with unit covariance (unless specified otherwise by setting a custom FixedDiffusion).

source
ProbNumDiffEq.ClassicSolverInit — Type
ClassicSolverInit(; alg=OrdinaryDiffEqCore.Tsit5(), init_on_ddu=false)

Initialization via regression on a few steps of a classic ODE solver.

In a nutshell, instead of specifying $\mu_0$ exactly and setting $\Sigma_0=0$ (which is what TaylorModeInit does), use a classic ODE solver to compute a few steps of the solution, and then regress on the computed values (by running a smoother) to compute $\mu_0$ and $\Sigma_0$ as the mean and covariance of the smoothing posterior at time 0. See also [2].

The initial value and derivative are set directly from the given initial value problem; optionally the second derivative can also be set via automatic differentiation by setting init_on_ddu=true.

Arguments

  • alg: The solver to be used. Can be any solver from OrdinaryDiffEqCore.jl.
  • init_on_ddu: If true, the second derivative is also initialized exactly via automatic differentiation with ForwardDiff.jl.

References

  • [7] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
  • [8] Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)
source

References

[7]
N. Krämer and P. Hennig. Stable Implementation of Probabilistic ODE Solvers. CoRR (2020), arXiv:2012.10106 [stat.ML].
[8]
M. Schober, S. Särkkä and P. Hennig. A probabilistic model for the numerical solution of initial value problems. Statistics and Computing 29, 99–122 (2019).
+\end{aligned}\]

It turns out that we can also compute higher-order derivatives of $y$ with the chain rule, and then use these to better initialize $Y^{(i)}(0)$. This, done efficiently with Taylor-mode automatic differentiation by using TaylorIntegration.jl, is what TaylorModeInit does. See also [1].

In the vast majority of cases, just stick to the exact Taylor-mode initialization TaylorModeInit!

API

ProbNumDiffEq.TaylorModeInit — Type
TaylorModeInit(order)

Exact initialization via Taylor-mode automatic differentiation up to order order.

This is the recommended initialization method!

It uses TaylorIntegration.jl to efficiently compute the higher-order derivatives of the solution at the initial value, via Taylor-mode automatic differentiation.

In some special cases it can happen that TaylorIntegration.jl is incompatible with the given problem (typically because the problem definition does not allow for elements of type Taylor). If this happens, try one of SimpleInit, ForwardDiffInit (for low enough orders), ClassicSolverInit.

References

  • [7] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
source
ProbNumDiffEq.ForwardDiffInit — Type
ForwardDiffInit(order)

Exact initialization via ForwardDiff.jl up to order order.

Warning: This does not scale well to high orders! For orders > 3, TaylorModeInit most likely performs better.

source
ProbNumDiffEq.SimpleInit — Type
SimpleInit()

Simple initialization, only with the given initial value and derivative.

The remaining derivatives are set to zero with unit covariance (unless specified otherwise by setting a custom FixedDiffusion).

source
ProbNumDiffEq.ClassicSolverInit — Type
ClassicSolverInit(; alg=OrdinaryDiffEqCore.Tsit5(), init_on_ddu=false)

Initialization via regression on a few steps of a classic ODE solver.

In a nutshell, instead of specifying $\mu_0$ exactly and setting $\Sigma_0=0$ (which is what TaylorModeInit does), use a classic ODE solver to compute a few steps of the solution, and then regress on the computed values (by running a smoother) to compute $\mu_0$ and $\Sigma_0$ as the mean and covariance of the smoothing posterior at time 0. See also [2].

The initial value and derivative are set directly from the given initial value problem; optionally the second derivative can also be set via automatic differentiation by setting init_on_ddu=true.

Arguments

  • alg: The solver to be used. Can be any solver from OrdinaryDiffEqCore.jl.
  • init_on_ddu: If true, the second derivative is also initialized exactly via automatic differentiation with ForwardDiff.jl.

References

  • [7] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
  • [8] Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)
source

References

[7]
N. Krämer and P. Hennig. Stable Implementation of Probabilistic ODE Solvers. CoRR (2020), arXiv:2012.10106 [stat.ML].
[8]
M. Schober, S. Särkkä and P. Hennig. A probabilistic model for the numerical solution of initial value problems. Statistics and Computing 29, 99–122 (2019).
diff --git a/dev/likelihoods/index.html b/dev/likelihoods/index.html index 145a89d28..65974a947 100644 --- a/dev/likelihoods/index.html +++ b/dev/likelihoods/index.html @@ -8,7 +8,7 @@ data, kwargs... ) -

Compute the Fenrir [3] approximate negative log-likelihood (NLL) of the data.

This is a convenience function that

  1. Solves the ODE with a ProbNumDiffEq.EK1 of the specified order and with a diffusion as provided by the diffusion_var argument, and
  2. Fits the ODE posterior to the data via Kalman filtering and thereby computes the log-likelihood of the data on the way.

You can control the step-size behaviour of the solver as you would for a standard ODE solve, but additionally the solver always steps through the data.t locations by adding them to tstops.

You can also choose steps adaptively by setting adaptive=true, but this is not well-tested so use at your own risk!

Arguments

Reference

source
ProbNumDiffEq.dalton_data_loglik — Function
dalton_data_loglik(
+

Compute the Fenrir [3] approximate negative log-likelihood (NLL) of the data.

This is a convenience function that

  1. Solves the ODE with a ProbNumDiffEq.EK1 of the specified order and with a diffusion as provided by the diffusion_var argument, and
  2. Fits the ODE posterior to the data via Kalman filtering and thereby computes the log-likelihood of the data on the way.

You can control the step-size behaviour of the solver as you would for a standard ODE solve, but additionally the solver always steps through the data.t locations by adding them to tstops.

You can also choose steps adaptively by setting adaptive=true, but this is not well-tested so use at your own risk!

Arguments

  • prob::SciMLBase.AbstractODEProblem: the initial value problem of interest
  • alg::AbstractEK: the probabilistic ODE solver to be used; use EK1 for best results.
  • data::NamedTuple{(:t, :u)}: the data to be fitted
  • observation_matrix::Union{AbstractMatrix,UniformScaling}: the matrix which maps the ODE state to the measurements; typically a projection matrix
  • observation_noise_cov::Union{Number,AbstractMatrix}: the scalar observation noise variance

Reference

  • [3] Tronarp et al, "Fenrir: Physics-Enhanced Regression for Initial Value Problems", ICML (2022)
source
ProbNumDiffEq.dalton_data_loglik — Function
dalton_data_loglik(
     prob::SciMLBase.AbstractODEProblem,
     alg::ProbNumDiffEq.AbstractEK,
     args...;
@@ -17,4 +17,4 @@
     data,
     kwargs...
 )
-

Compute the DALTON [10] approximate negative log-likelihood (NLL) of the data.

You can control the step-size behaviour of the solver as you would for a standard ODE solve, but additionally the solver always steps through the data.t locations by adding them to tstops. You can also choose steps adaptively by setting adaptive=true, but this is not well-tested so use at your own risk!

Arguments

  • prob::SciMLBase.AbstractODEProblem: the initial value problem of interest
  • alg::AbstractEK: the probabilistic ODE solver to be used; use EK1 for best results.
  • data::NamedTuple{(:t, :u)}: the data to be fitted
  • observation_matrix::Union{AbstractMatrix,UniformScaling}: the matrix which maps the ODE state to the measurements; typically a projection matrix
  • observation_noise_cov::Union{Number,AbstractMatrix}: the scalar observation noise variance

Reference

  • [10] Wu et al, "Data-Adaptive Probabilistic Likelihood Approximation for Ordinary Differential Equations", arXiv (2023)
source
+

Compute the DALTON [10] approximate negative log-likelihood (NLL) of the data.

You can control the step-size behaviour of the solver as you would for a standard ODE solve, but additionally the solver always steps through the data.t locations by adding them to tstops. You can also choose steps adaptively by setting adaptive=true, but this is not well-tested so use at your own risk!

Arguments

Reference

source diff --git a/dev/objects.inv b/dev/objects.inv index 208f5127c9964f21c5268d669c3ec585cfae3701..ba3d96e3ac979df1adce5a8e7d71e81e4e6d62ac 100644 GIT binary patch delta 12 TcmaDa{$6~7C!^6uuch1oBG?4P delta 12 TcmaDa{$6~7C!^s;uch1oBGUxJ diff --git a/dev/priors/49473a66.svg b/dev/priors/49473a66.svg deleted file mode 100644 index 1b24893aa..000000000 --- a/dev/priors/49473a66.svg +++ /dev/null @@ -1,174 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/priors/53661c5c.svg b/dev/priors/53661c5c.svg new file mode 100644 index 000000000..a12b562b7 --- /dev/null +++ b/dev/priors/53661c5c.svg @@ -0,0 +1,174 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/priors/69865de8.svg b/dev/priors/69865de8.svg deleted file mode 100644 index 8f0066f8a..000000000 --- a/dev/priors/69865de8.svg +++ /dev/null @@ -1,174 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/priors/78c44fc4.svg b/dev/priors/78c44fc4.svg new file mode 100644 index 000000000..082adc66b --- /dev/null +++ b/dev/priors/78c44fc4.svg @@ -0,0 +1,106 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/priors/8d30ced6.svg b/dev/priors/8d30ced6.svg deleted file mode 100644 index a015a27bc..000000000 --- a/dev/priors/8d30ced6.svg +++ /dev/null @@ -1,174 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/priors/a292ceb1.svg b/dev/priors/a292ceb1.svg deleted file mode 100644 index d2f5ad13c..000000000 --- a/dev/priors/a292ceb1.svg +++ /dev/null @@ -1,104 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/priors/a350c826.svg b/dev/priors/a350c826.svg new file mode 100644 index 000000000..39cf95ca6 --- /dev/null +++ b/dev/priors/a350c826.svg @@ -0,0 +1,174 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/priors/d41d5e0e.svg b/dev/priors/d41d5e0e.svg new file mode 100644 index 000000000..063b53a14 --- /dev/null +++ b/dev/priors/d41d5e0e.svg @@ -0,0 +1,174 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/priors/d84323e8.svg b/dev/priors/d84323e8.svg deleted file mode 100644 index 789e6e0a0..000000000 --- a/dev/priors/d84323e8.svg +++ /dev/null @@ -1,178 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/priors/e224c6c6.svg b/dev/priors/e224c6c6.svg new file mode 100644 index 000000000..e6447d192 --- /dev/null +++ b/dev/priors/e224c6c6.svg @@ -0,0 +1,182 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/priors/index.html b/dev/priors/index.html index 2fff8522e..9e73f9f1d 100644 --- a/dev/priors/index.html +++ b/dev/priors/index.html @@ -10,7 +10,7 @@ \end{aligned}\]

defined as the solution of the stochastic differential equation

\[\begin{aligned} \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1 \\ \text{d} Y^{(q)}(t) &= \Gamma \ \text{d}W(t). -\end{aligned}\]

Then, $Y^{(0)}(t)$ is the $q$-times integrated Wiener process (IWP) and $Y^{(i)}(t)$ is the $i$-th derivative of the IWP, for $i = 1, \dots, q$.

Examples

julia> solve(prob, EK1(prior=IWP(2)))
source

Here is how the IWP looks for varying smoothness parameters $q$:

using ProbNumDiffEq, Plots
+\end{aligned}\]
Then, $Y^{(0)}(t)$ is the $q$-times integrated Wiener process (IWP) and $Y^{(i)}(t)$ is the $i$-th derivative of the IWP, for $i = 1, \dots, q$.
Examples
julia> solve(prob, EK1(prior=IWP(2)))
source
Here is how the IWP looks for varying smoothness parameters $q$:
using ProbNumDiffEq, Plots
 plotrange = range(0, 10, length=250)
 plot(
     plot(IWP(1), plotrange; title="q=1"),
@@ -18,22 +18,22 @@
     plot(IWP(3), plotrange; title="q=3"),
     plot(IWP(4), plotrange; title="q=4");
     ylims=(-20,20),
-)
Example block output
In the context of ODE solvers, the smoothness parameter $q$ influences the convergence rate of the solver, and so it is typically chose similarly to the order of a Runge–Kutta method: lower order for low accuracy, higher order for high accuracy.
Integrated Ornstein–Uhlenbeck process (IOUP)
The $q$-times integrated Ornstein–Uhlenbeck process prior IOUP is a generalization of the IWP prior, where the drift matrix $\textcolor{#389826}{A}$ is not zero:
ProbNumDiffEq.IOUP — Type
IOUP([dim::Integer=1,]
+)
Example block output
In the context of ODE solvers, the smoothness parameter $q$ influences the convergence rate of the solver, and so it is typically chose similarly to the order of a Runge–Kutta method: lower order for low accuracy, higher order for high accuracy.
Integrated Ornstein–Uhlenbeck process (IOUP)
The $q$-times integrated Ornstein–Uhlenbeck process prior IOUP is a generalization of the IWP prior, where the drift matrix $\textcolor{#389826}{A}$ is not zero:
ProbNumDiffEq.IOUP — Type
IOUP([dim::Integer=1,]
      num_derivatives::Integer,
      rate_parameter::Union{Number,AbstractVector,AbstractMatrix})
Integrated Ornstein–Uhlenbeck process.
This prior is mostly used in the context of Probabilistic Exponential Integrators to include the linear part of a semi-linear ODE in the prior, so it is used in the ExpEK and the RosenbrockExpEK.
In math
The IOUP is a Gauss–Markov process, which we model with a state representation
\[\begin{aligned}
 Y(t) = \left[ Y^{(0)}(t), Y^{(1)}(t), \dots Y^{(q)}(t) \right],
 \end{aligned}\]
defined as the solution of the stochastic differential equation
\[\begin{aligned}
 \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1 \\
 \text{d} Y^{(q)}(t) &= L Y^{(q)}(t) \ \text{d}t + \Gamma \ \text{d}W(t),
-\end{aligned}\]
where $L$ is called the rate_parameter. Then, $Y^{(0)}(t)$ is the $q$-times integrated Ornstein–Uhlenbeck process (IOUP) and $Y^{(i)}(t)$ is the $i$-th derivative of the IOUP, for $i = 1, \dots, q$.
Examples
julia> solve(prob, EK1(prior=IOUP(2, -1)))
source
Here is how the IOUP looks for varying rate parameters:
using ProbNumDiffEq, Plots
+\end{aligned}\]
where $L$ is called the rate_parameter. Then, $Y^{(0)}(t)$ is the $q$-times integrated Ornstein–Uhlenbeck process (IOUP) and $Y^{(i)}(t)$ is the $i$-th derivative of the IOUP, for $i = 1, \dots, q$.
Examples
julia> solve(prob, EK1(prior=IOUP(2, -1)))
source
Here is how the IOUP looks for varying rate parameters:
using ProbNumDiffEq, Plots
 plotrange = range(0, 10, length=250)
 plot(
     plot(IOUP(1, -1), plotrange; title="q=1,L=-1", ylims=(-20,20)),
     plot(IOUP(1, 1), plotrange; title="q=1,L=1", ylims=(-20,20)),
     plot(IOUP(4, -1), plotrange; title="q=4,L=-1", ylims=(-50,50)),
     plot(IOUP(4, 1), plotrange; title="q=4,L=1", ylims=(-50,50));
-)
Example block output
In the context of Probabilistic Exponential Integrators, the rate parameter is often chosen according to the given ODE. Here is an example for a damped oscillator:
plot(IOUP(dim=2, num_derivatives=1, rate_parameter=[-0.2 -2π; 2π -0.2]),
-     plotrange; plot_title="damped oscillator prior")
Example block output
Matérn process (Matern)
ProbNumDiffEq.Matern — Type
Matern([dim::Integer=1,]
+)
Example block output
In the context of Probabilistic Exponential Integrators, the rate parameter is often chosen according to the given ODE. Here is an example for a damped oscillator:
plot(IOUP(dim=2, num_derivatives=1, rate_parameter=[-0.2 -2π; 2π -0.2]),
+     plotrange; plot_title="damped oscillator prior")
Example block output
Matérn process (Matern)
ProbNumDiffEq.Matern — Type
Matern([dim::Integer=1,]
        num_derivatives::Integer,
        lengthscale::Number)
Matern process.
The class of Matern processes is well-known in the Gaussian process literature, and they also have a corresponding SDE representation similarly to the IWP and the IOUP. See also [6] for more details.
In math
A Matern process is a Gauss–Markov process, which we model with a state representation
\[\begin{aligned}
 Y(t) = \left[ Y^{(0)}(t), Y^{(1)}(t), \dots Y^{(q)}(t) \right],
@@ -43,22 +43,22 @@
   \begin{pmatrix} q+1 \\ j \end{pmatrix}
   \left( \frac{\sqrt{2q - 1}}{l} \right)^{q-j}
   Y^{(j)}(t) \right) \ \text{d}t + \Gamma \ \text{d}W(t).
-\end{aligned}\]
where $l$ is called the lengthscale parameter. Then, $Y^{(0)}(t)$ is a Matern process and $Y^{(i)}(t)$ is the $i$-th derivative of this process, for $i = 1, \dots, q$.
Examples
julia> solve(prob, EK1(prior=Matern(2, 1)))
source
Here is how the Matern looks for varying smoothness parameters $q$:
using ProbNumDiffEq, Plots
+\end{aligned}\]
where $l$ is called the lengthscale parameter. Then, $Y^{(0)}(t)$ is a Matern process and $Y^{(i)}(t)$ is the $i$-th derivative of this process, for $i = 1, \dots, q$.
Examples
julia> solve(prob, EK1(prior=Matern(2, 1)))
source
Here is how the Matern looks for varying smoothness parameters $q$:
using ProbNumDiffEq, Plots
 plotrange = range(0, 10, length=250)
 plot(
     plot(Matern(1, 1), plotrange; title="q=1"),
     plot(Matern(2, 1), plotrange; title="q=2"),
     plot(Matern(3, 1), plotrange; title="q=3"),
     plot(Matern(4, 1), plotrange; title="q=4");
-)
Example block output
and for varying length scales $\ell$:
plot(
+)
Example block output
and for varying length scales $\ell$:
plot(
     plot(Matern(2, 5), plotrange; title="l=5"),
     plot(Matern(2, 2), plotrange; title="l=2"),
     plot(Matern(2, 1), plotrange; title="l=1"),
     plot(Matern(2, 0.5), plotrange; title="l=0.5"),
-)
Example block output
API
ProbNumDiffEq.AbstractGaussMarkovProcess — Type
AbstractGaussMarkovProcess{elType}
Abstract type for Gauss-Markov processes.
Gauss-Markov processes are solutions to linear time-invariant stochastic differential equations (SDEs). Here we assume SDEs of the form
\[\begin{aligned}
+)
Example block output

API

ProbNumDiffEq.AbstractGaussMarkovProcess — Type
AbstractGaussMarkovProcess{elType}

Abstract type for Gauss-Markov processes.

Gauss-Markov processes are solutions to linear time-invariant stochastic differential equations (SDEs). Here we assume SDEs of the form

\[\begin{aligned} dX_t &= F X_t dt + L dW_t \\ X_0 &= \mathcal{N} \left( X_0; \mu_0, \Sigma_0 \right) -\end{aligned}\]

where $X_t$ is the state, $W_t$ is a Wiener process, and $F$ and $L$ are matrices.

Currently, ProbNumDiffEq.jl makes many assumptions about the structure of the SDEs that it can solve. In particular, it assumes that the state vector $X_t$ contains a range of dervatives, and that the Wiener process only enters the highest one. It also assumes a certain ordering of dimensions and derivatives. This is not a limitation of the underlying mathematics, but rather a limitation of the current implementation. In the future, we hope to remove these limitations.

Warning

We currently strongly recommended to not implement your own Gauss-Markov process by subtyping this type! The interface is not yet stable, and the implementation is not yet sufficiently documented. Proceed at your own risk.

source
ProbNumDiffEq.LTISDE — Type
LTISDE(F::AbstractMatrix, L::AbstractMatrix)

Linear time-invariant stochastic differential equation.

A LTI-SDE is a stochastic differential equation of the form

\[dX_t = F X_t dt + L dW_t\]

where $X_t$ is the state, $W_t$ is a Wiener process, and $F$ and $L$ are matrices. This LTISDE object holds the matrices $F$ and $L$. It also provides some functionality to discretize the SDE via a matrix-fraction decomposition. See: discretize(::LTISDE, ::Real).

source
ProbNumDiffEq.dim — Function
dim(p::AbstractGaussMarkovProcess)

Return the dimension of the process.

This is not the dimension of the "state" that is used to efficiently model the prior process as a state-space model, but it is the dimension of the process itself that we aim to model.

See AbstractGaussMarkovProcess for more details on Gauss-Markov processes in ProbNumDiffEq.

source
ProbNumDiffEq.to_sde — Function
to_sde(p::AbstractGaussMarkovProcess)

Convert the prior to the corresponding SDE.

Gauss-Markov processes are solutions to linear time-invariant stochastic differential equations (SDEs) of the form

\[\begin{aligned} +\end{aligned}\]

where $X_t$ is the state, $W_t$ is a Wiener process, and $F$ and $L$ are matrices.

Currently, ProbNumDiffEq.jl makes many assumptions about the structure of the SDEs that it can solve. In particular, it assumes that the state vector $X_t$ contains a range of dervatives, and that the Wiener process only enters the highest one. It also assumes a certain ordering of dimensions and derivatives. This is not a limitation of the underlying mathematics, but rather a limitation of the current implementation. In the future, we hope to remove these limitations.

Warning

We currently strongly recommended to not implement your own Gauss-Markov process by subtyping this type! The interface is not yet stable, and the implementation is not yet sufficiently documented. Proceed at your own risk.

source
ProbNumDiffEq.LTISDE — Type
LTISDE(F::AbstractMatrix, L::AbstractMatrix)

Linear time-invariant stochastic differential equation.

A LTI-SDE is a stochastic differential equation of the form

\[dX_t = F X_t dt + L dW_t\]

where $X_t$ is the state, $W_t$ is a Wiener process, and $F$ and $L$ are matrices. This LTISDE object holds the matrices $F$ and $L$. It also provides some functionality to discretize the SDE via a matrix-fraction decomposition. See: discretize(::LTISDE, ::Real).

source
ProbNumDiffEq.dim — Function
dim(p::AbstractGaussMarkovProcess)

Return the dimension of the process.

This is not the dimension of the "state" that is used to efficiently model the prior process as a state-space model, but it is the dimension of the process itself that we aim to model.

See AbstractGaussMarkovProcess for more details on Gauss-Markov processes in ProbNumDiffEq.

source
ProbNumDiffEq.to_sde — Function
to_sde(p::AbstractGaussMarkovProcess)

Convert the prior to the corresponding SDE.

Gauss-Markov processes are solutions to linear time-invariant stochastic differential equations (SDEs) of the form

\[\begin{aligned} dX_t &= F X_t dt + L dW_t \\ X_0 &= \mathcal{N} \left( X_0; \mu_0, \Sigma_0 \right) -\end{aligned}\]

where $X_t$ is the state, $W_t$ is a Wiener process, and $F$ and $L$ are matrices. This function returns the corresponding SDE, i.e. the matrices $F$ and $L$, as a LTISDE.

source
ProbNumDiffEq.discretize — Function
discretize(p::AbstractGaussMarkovProcess, step_size::Real)

Compute the transition matrices of the process for a given step size.

source
discretize(p::LTISDE, step_size::Real)

Compute the transition matrices of the SDE solution for a given step size.

source
ProbNumDiffEq.initial_distribution — Function
initial_distribution(p::AbstractGaussMarkovProcess)

Return the initial distribution of the process.

Currently this is always a Gaussian distribution with zero mean and unit variance, unless explicitly overwitten (e.g. for Matern processes to have the stationary distribution). This implementation is likely to change in the future to allow for more flexibility.

source

Convenience functions to analyze and visualize priors

ProbNumDiffEq.marginalize — Function
marginalize(process::AbstractGaussMarkovProcess, times)

Compute the marginal distributions of the process at the given time points.

This function computes the marginal distributions of the process at the given times. It does so by discretizing the process with the given step sizes (using ProbNumDiffEq.discretize), and then computing the marginal distributions of the resulting Gaussian distributions.

See also: sample.

source
+\end{aligned}\]

where $X_t$ is the state, $W_t$ is a Wiener process, and $F$ and $L$ are matrices. This function returns the corresponding SDE, i.e. the matrices $F$ and $L$, as a LTISDE.

source
ProbNumDiffEq.discretize — Function
discretize(p::AbstractGaussMarkovProcess, step_size::Real)

Compute the transition matrices of the process for a given step size.

source
discretize(p::LTISDE, step_size::Real)

Compute the transition matrices of the SDE solution for a given step size.

source
ProbNumDiffEq.initial_distribution — Function
initial_distribution(p::AbstractGaussMarkovProcess)

Return the initial distribution of the process.

Currently this is always a Gaussian distribution with zero mean and unit variance, unless explicitly overwitten (e.g. for Matern processes to have the stationary distribution). This implementation is likely to change in the future to allow for more flexibility.

source

Convenience functions to analyze and visualize priors

ProbNumDiffEq.marginalize — Function
marginalize(process::AbstractGaussMarkovProcess, times)

Compute the marginal distributions of the process at the given time points.

This function computes the marginal distributions of the process at the given times. It does so by discretizing the process with the given step sizes (using ProbNumDiffEq.discretize), and then computing the marginal distributions of the resulting Gaussian distributions.

See also: sample.

source
diff --git a/dev/references/index.html b/dev/references/index.html index 3b456b993..f73d4f887 100644 --- a/dev/references/index.html +++ b/dev/references/index.html @@ -1,2 +1,2 @@ -References · ProbNumDiffEq.jl
GitHub

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Vol. 151 of Proceedings of Machine Learning Research, edited by G. Camps-Valls, F. J. Ruiz and I. Valera (PMLR, 28–30 Mar 2022); pp. 10015–10027.
[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).
[3]
F. Tronarp, N. Bosch and P. Hennig. Fenrir: Physics-Enhanced Regression for Initial Value Problems. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 21776–21794.
[4]
[5]
N. Krämer, N. Bosch, J. Schmidt and P. Hennig. Probabilistic ODE Solutions in Millions of Dimensions. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 11634–11649.
[6]
S. Särkkä and A. Solin. Applied Stochastic Differential Equations. Institute of Mathematical Statistics Textbooks (Cambridge University Press, 2019).
[7]
[8]
[9]
N. Bosch, P. Hennig and F. Tronarp. Calibrated Adaptive Probabilistic ODE Solvers. In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, Vol. 130 of Proceedings of Machine Learning Research, edited by A. Banerjee and K. Fukumizu (PMLR, 13–15 Apr 2021); pp. 3466–3474.
[10]
[11]
P. Hennig, M. A. Osborne and H. P. Kersting. Probabilistic Numerics: Computation as Machine Learning (Cambridge University Press, 2022).
[12]
H. Kersting and P. Hennig. Active Uncertainty Calibration in Bayesian ODE Solvers. In: Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence, UAI'16 (AUAI Press, 2016); pp. 309–318.
[13]
H. Kersting, T. J. Sullivan and P. Hennig. Convergence rates of Gaussian ODE filters. Statistics and Computing 30, 1791–1816 (2020).
[14]
+References · ProbNumDiffEq.jl
GitHub

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Vol. 151 of Proceedings of Machine Learning Research, edited by G. Camps-Valls, F. J. Ruiz and I. Valera (PMLR, 28–30 Mar 2022); pp. 10015–10027.
[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).
[3]
F. Tronarp, N. Bosch and P. Hennig. Fenrir: Physics-Enhanced Regression for Initial Value Problems. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 21776–21794.
[4]
[5]
N. Krämer, N. Bosch, J. Schmidt and P. Hennig. Probabilistic ODE Solutions in Millions of Dimensions. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 11634–11649.
[6]
S. Särkkä and A. Solin. Applied Stochastic Differential Equations. Institute of Mathematical Statistics Textbooks (Cambridge University Press, 2019).
[7]
[8]
[9]
N. Bosch, P. Hennig and F. Tronarp. Calibrated Adaptive Probabilistic ODE Solvers. In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, Vol. 130 of Proceedings of Machine Learning Research, edited by A. Banerjee and K. Fukumizu (PMLR, 13–15 Apr 2021); pp. 3466–3474.
[10]
[11]
P. Hennig, M. A. Osborne and H. P. Kersting. Probabilistic Numerics: Computation as Machine Learning (Cambridge University Press, 2022).
[12]
H. Kersting and P. Hennig. Active Uncertainty Calibration in Bayesian ODE Solvers. In: Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence, UAI'16 (AUAI Press, 2016); pp. 309–318.
[13]
H. Kersting, T. J. Sullivan and P. Hennig. Convergence rates of Gaussian ODE filters. Statistics and Computing 30, 1791–1816 (2020).
[14]
diff --git a/dev/solvers/index.html b/dev/solvers/index.html index 22cf8254f..1fb324e45 100644 --- a/dev/solvers/index.html +++ b/dev/solvers/index.html @@ -4,10 +4,10 @@ prior=IWP(order), diffusionmodel=DynamicDiffusion(), initialization=TaylorModeInit(num_derivatives(prior)), - kwargs...)

Gaussian ODE filter with first-order vector field linearization.

This is a semi-implicit, L-stable ODE solver so it can handle stiffness quite well [4], and it generally produces more expressive posterior covariances than the EK0. However, as typical implicit ODE solvers it scales cubically with the ODE dimension [5], so if you're solving a high-dimensional non-stiff problem you might want to give the EK0 a try.

Arguments

Some additional kwargs relating to implicit solvers are supported; check out DifferentialEquations.jl's Extra Options page. Right now, we support autodiff, chunk_size, and diff_type. In particular, autodiff=false can come in handy to use finite differences instead of ForwardDiff.jl to compute Jacobians.

Examples

julia> solve(prob, EK1())

References

source
ProbNumDiffEq.EK0 — Type
EK0(; order=3,
+      kwargs...)

Gaussian ODE filter with first-order vector field linearization.

This is a semi-implicit, L-stable ODE solver so it can handle stiffness quite well [4], and it generally produces more expressive posterior covariances than the EK0. However, as typical implicit ODE solvers it scales cubically with the ODE dimension [5], so if you're solving a high-dimensional non-stiff problem you might want to give the EK0 a try.

Arguments

  • order::Integer: Order of the integrated Wiener process (IWP) prior.
  • smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
  • prior::AbstractGaussMarkovProcess: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
  • diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
  • initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Some additional kwargs relating to implicit solvers are supported; check out DifferentialEquations.jl's Extra Options page. Right now, we support autodiff, chunk_size, and diff_type. In particular, autodiff=false can come in handy to use finite differences instead of ForwardDiff.jl to compute Jacobians.

Examples

julia> solve(prob, EK1())

References

source
ProbNumDiffEq.EK0 — Type
EK0(; order=3,
       smooth=true,
       prior=IWP(order),
       diffusionmodel=DynamicDiffusion(),
-      initialization=TaylorModeInit(num_derivatives(prior)))

Gaussian ODE filter with zeroth-order vector field linearization.

This is an explicit ODE solver. It is fast and scales well to high-dimensional problems [5], but it is not L-stable [4]. So for stiff problems, use the EK1.

Whenever possible this solver will use a Kronecker-factored implementation to achieve its linear scaling and to get the best runtimes. This can currently be done only with an IWP prior (default), with a scalar diffusion model (either DynamicDiffusion or FixedDiffusion). For other configurations the solver falls back to a dense implementation which scales cubically with the problem size.

Arguments

  • order::Integer: Order of the integrated Wiener process (IWP) prior.
  • smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
  • prior::AbstractGaussMarkovProcess: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
  • diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
  • initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Examples

julia> solve(prob, EK0())

References

source

Probabilistic Exponential Integrators

ProbNumDiffEq.ExpEK — Function
ExpEK(; L, order=3, kwargs...)

Probabilistic exponential integrator

Probabilistic exponential integrators are a class of integrators for semi-linear stiff ODEs that provide improved stability by essentially solving the linear part of the ODE exactly. In probabilistic numerics, this amounts to including the linear part into the prior model of the solver.

ExpEK is therefore just a short-hand for EK0 with IOUP prior:

ExpEK(; order=3, L, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)

See also RosenbrockExpEK, EK0, EK1.

Arguments

See EK0 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
-julia> solve(prob, ExpEK(L=-1))

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
source
ProbNumDiffEq.RosenbrockExpEK — Function
RosenbrockExpEK(; order=3, kwargs...)

Probabilistic Rosenbrock-type exponential integrator

A probabilistic exponential integrator similar to ExpEK, but with automatic linearization along the mean numerical solution. This brings the advantage that the linearity does not need to be specified manually, and the more accurate local linearization can sometimes also improve stability; but since the "prior" is adjusted at each step the probabilistic interpretation becomes more complicated.

RosenbrockExpEK is just a short-hand for EK1 with locally-updated IOUP prior:

RosenbrockExpEK(; order=3, kwargs...) = EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)

See also ExpEK, EK0, EK1.

Arguments

See EK1 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
-julia> solve(prob, RosenbrockExpEK())

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", NeurIPS (2023)
source

References

[4]
F. Tronarp, H. Kersting, S. Särkkä and P. Hennig. Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filtering: a new perspective. Statistics and Computing 29, 1297–1315 (2019).
[5]
N. Krämer, N. Bosch, J. Schmidt and P. Hennig. Probabilistic ODE Solutions in Millions of Dimensions. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 11634–11649.
[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).
+ initialization=TaylorModeInit(num_derivatives(prior)))

Gaussian ODE filter with zeroth-order vector field linearization.

This is an explicit ODE solver. It is fast and scales well to high-dimensional problems [5], but it is not L-stable [4]. So for stiff problems, use the EK1.

Whenever possible this solver will use a Kronecker-factored implementation to achieve its linear scaling and to get the best runtimes. This can currently be done only with an IWP prior (default), with a scalar diffusion model (either DynamicDiffusion or FixedDiffusion). For other configurations the solver falls back to a dense implementation which scales cubically with the problem size.

Arguments

Examples

julia> solve(prob, EK0())

References

source

Probabilistic Exponential Integrators

ProbNumDiffEq.ExpEK — Function
ExpEK(; L, order=3, kwargs...)

Probabilistic exponential integrator

Probabilistic exponential integrators are a class of integrators for semi-linear stiff ODEs that provide improved stability by essentially solving the linear part of the ODE exactly. In probabilistic numerics, this amounts to including the linear part into the prior model of the solver.

ExpEK is therefore just a short-hand for EK0 with IOUP prior:

ExpEK(; order=3, L, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)

See also RosenbrockExpEK, EK0, EK1.

Arguments

See EK0 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
+julia> solve(prob, ExpEK(L=-1))

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
source
ProbNumDiffEq.RosenbrockExpEK — Function
RosenbrockExpEK(; order=3, kwargs...)

Probabilistic Rosenbrock-type exponential integrator

A probabilistic exponential integrator similar to ExpEK, but with automatic linearization along the mean numerical solution. This brings the advantage that the linearity does not need to be specified manually, and the more accurate local linearization can sometimes also improve stability; but since the "prior" is adjusted at each step the probabilistic interpretation becomes more complicated.

RosenbrockExpEK is just a short-hand for EK1 with locally-updated IOUP prior:

RosenbrockExpEK(; order=3, kwargs...) = EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)

See also ExpEK, EK0, EK1.

Arguments

See EK1 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
+julia> solve(prob, RosenbrockExpEK())

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", NeurIPS (2023)
source

References

[4]
F. Tronarp, H. Kersting, S. Särkkä and P. Hennig. Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filtering: a new perspective. Statistics and Computing 29, 1297–1315 (2019).
[5]
N. Krämer, N. Bosch, J. Schmidt and P. Hennig. Probabilistic ODE Solutions in Millions of Dimensions. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 11634–11649.
[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).
diff --git a/dev/tutorials/dae/3542b1a5.svg b/dev/tutorials/dae/3542b1a5.svg new file mode 100644 index 000000000..3b1415436 --- /dev/null +++ b/dev/tutorials/dae/3542b1a5.svg @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dae/b933a08f.svg b/dev/tutorials/dae/b933a08f.svg new file mode 100644 index 000000000..2dd3ef5d3 --- /dev/null +++ b/dev/tutorials/dae/b933a08f.svg @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dae/bbfc7e14.svg b/dev/tutorials/dae/bbfc7e14.svg deleted file mode 100644 index 45bea5c96..000000000 --- a/dev/tutorials/dae/bbfc7e14.svg +++ /dev/null @@ -1,59 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/dae/eec28457.svg b/dev/tutorials/dae/eec28457.svg deleted file mode 100644 index e0402ef96..000000000 --- a/dev/tutorials/dae/eec28457.svg +++ /dev/null @@ -1,90 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/dae/index.html b/dev/tutorials/dae/index.html index ca334e813..d4a998d8a 100644 --- a/dev/tutorials/dae/index.html +++ b/dev/tutorials/dae/index.html @@ -28,7 +28,7 @@ ylabel=["u₁(t)" "u₂(t)" "u₃(t)"], xlabel=["" "" "t"], denseplot=false, -)Example block output

Looks good!

Solving an Index-3 DAE directly

The following is based on the "Automatic Index Reduction of DAEs" tutorial by ModelingToolkit.jl, which demonstrates how the classic Rodas4 solver fails to solve a DAE due to the fact that it is of index 3; which is why ModelingToolkit's automatic index reduction is so useful.

It turns out that our probabilistic numerical solvers can directly solve the index-3 DAE!

First, define the pendulum problem as in the tutorial:

function pendulum!(du, u, p, t)
+)
Example block output

Looks good!

Solving an Index-3 DAE directly

The following is based on the "Automatic Index Reduction of DAEs" tutorial by ModelingToolkit.jl, which demonstrates how the classic Rodas4 solver fails to solve a DAE due to the fact that it is of index 3; which is why ModelingToolkit's automatic index reduction is so useful.

It turns out that our probabilistic numerical solvers can directly solve the index-3 DAE!

First, define the pendulum problem as in the tutorial:

function pendulum!(du, u, p, t)
     x, dx, y, dy, T = u
     g, L = p
     du[1] = dx
@@ -54,57 +54,57 @@
  0.0
  1.0e-6
  1.1e-5
- 3.806859510273279e-5
- 4.367846401990628e-5
+ 3.806859541877657e-5
+ 4.5496933032029966e-5
 u: 5-element Vector{Vector{Float64}}:
  [1.0, 0.0, 0.0, 0.0, 0.0]
  [1.0, 0.0, -4.899999999999997e-12, -9.799999999999994e-6, 0.0]
  [1.0, -5.282199999999992e-16, -5.928999999999995e-10, -0.00010779999999999988, -9.04785641150889e-11]
- [1.0, -2.2119452635682024e-13, -7.1011678721694636e-9, -0.000373072232006781, -2.6515338160403213e-8]
- [1.0, -3.0396474329579032e-12, -9.348260273777398e-9, -0.00042804894739508116, -3.896506115963402e-6]

It does not work! This is because of the index of the DAE; see for example this explanation from the tutorial.

Does this also hold for the EK1 solver? Let's find out:

sol = solve(pendulum_prob, EK1())
retcode: Success
+ [1.0, -2.2119452819316723e-13, -7.101167990076622e-9, -0.00037307223510401005, -2.651533791865963e-8]
+ [1.0, -3.1661973720383655e-12, -1.0142857485072981e-8, -0.00044586994371389327, -2.9596714810803877e-6]

It does not work! This is because of the index of the DAE; see for example this explanation from the tutorial.

Does this also hold for the EK1 solver? Let's find out:

sol = solve(pendulum_prob, EK1())
retcode: Success
 Interpolation: ODE Filter Posterior
-t: 650-element Vector{Float64}:
+t: 637-element Vector{Float64}:
  0.0
  1.0e-6
  1.1e-5
- 6.293457400357088e-5
- 0.0001717966984100629
- 0.00041700371867394907
- 0.0010432618632924343
- 0.0023836933769503994
- 0.00376397894518298
- 0.00519132985680733
+ 6.293457247564385e-5
+ 0.0001717966902698356
+ 0.00041700284549204755
+ 0.001043224623113885
+ 0.0023835910941641687
+ 0.003763856534047043
+ 0.005191196217624695
  ⋮
- 4.907142111362168
- 4.92444553498383
- 4.9411116199936265
- 4.956803287770239
- 4.97098674548783
- 4.9822968246165
- 4.99143218593388
- 4.998834745818967
+ 4.905981142344071
+ 4.923072046379877
+ 4.939712741084385
+ 4.955319566223264
+ 4.9694706451654245
+ 4.982014085482488
+ 4.991556791257769
+ 4.999028623538522
  5.0
-u: 650-element Vector{Vector{Float64}}:
+u: 637-element Vector{Vector{Float64}}:
  [1.0, 0.0, 0.0, 0.0, 0.0]
- [1.0, -2.4314864246071144e-18, -4.899999999937332e-12, -9.800000000000003e-6, -3.870090693653406e-16]
- [1.0, 1.778914640672942e-14, -5.929000004468074e-10, -0.00010780000000000091, -3.017617371928309e-12]
- [0.9999999999999998, -3.704393921873193e-12, -1.940772692127836e-8, -0.0006167588252348035, -1.3035790293759336e-9]
- [0.9999999999999896, 1.2314212502343228e-11, -1.4461911744623022e-7, -0.0016836076444193238, -3.157756035951345e-8]
- [0.9999999999996368, -9.146450513208143e-11, -8.520712963690976e-7, -0.004086636443002532, -4.845250915236145e-7]
- [0.9999999999857788, -1.8913629380116342e-9, -5.333137047544597e-6, -0.01022396626026856, -7.830195678929709e-6]
- [0.9999999996124179, -5.727074650413097e-8, -2.7841771139450103e-5, -0.023360195092905718, -9.453725069261283e-5]
- [0.9999999975903671, -3.4756550064331693e-7, -6.942093370364985e-5, -0.03688699364700587, -0.00037351474250189485]
- [0.9999999912808006, -1.274813756871893e-6, -0.0001320545364564354, -0.05087503248488846, -0.0009815791650569213]
+ [1.0, -2.4320529399291366e-18, -4.899999999937319e-12, -9.800000000000003e-6, -3.8704386586000607e-16]
+ [1.0, 1.779329406237433e-14, -5.929000004468982e-10, -0.00010780000000000091, -3.0178910935929276e-12]
+ [0.9999999999999998, -3.705258905142665e-12, -1.9407725978899373e-8, -0.0006167588102611186, -1.3036995276752237e-9]
+ [0.9999999999999896, 1.2317344504921522e-11, -1.4461910374134008e-7, -0.001683607564645096, -3.158050959193969e-8]
+ [0.9999999999996371, -9.14785818709879e-11, -8.520677279956676e-7, -0.0040866278858198955, -4.845676419218797e-7]
+ [0.9999999999857808, -1.89123361432229e-9, -5.332756312011664e-6, -0.010223601306518774, -7.83009150843562e-6]
+ [0.9999999996124844, -5.726676667507785e-8, -2.7839381844866408e-5, -0.023359192721600835, -9.453399402972135e-5]
+ [0.9999999975906804, -3.475524673632411e-7, -6.941641839846304e-5, -0.0368857940178753, -0.00037351360952485107]
+ [0.9999999912816985, -1.2748041584224705e-6, -0.0001320477376462588, -0.050873722820905354, -0.0009815962497527718]
  ⋮
- [0.9902628022671714, -0.230012218970361, -0.13921057057397462, -1.6366422719422729, -4.0945422713871436]
- [0.9856333315526842, -0.3073677975051621, -0.16889919052551622, -1.7941346919319325, -4.968908564858052]
- [0.9797897702820494, -0.39617679785706794, -0.20003001727924735, -1.9407717381750471, -5.888064423522993]
- [0.972827745793681, -0.49345240804801094, -0.23153008794499258, -2.073030808278413, -6.821514288820706]
- [0.9651366055329357, -0.593329573335952, -0.2617467158743893, -2.186589633729141, -7.72111506539165]
- [0.957940002319977, -0.6813737682784675, -0.28696857832357914, -2.2722969869869667, -8.475874423432991]
- [0.9513764315907425, -0.7580420941403236, -0.3080306846666298, -2.337925795788078, -9.10937542703355]
- [0.9455309355620407, -0.8238465818165532, -0.325532285560256, -2.388489111563977, -9.638156836814543]
- [0.944566419026587, -0.8345069039015303, -0.3283203924231612, -2.3962208016104234, -9.722642388563708]

Nope! The EK1 is able to solve the index-3 DAE directly. Pretty cool!

plot(sol)
Example block output

Is index-reduction still worth it?

The point of the "Automatic Index Reduction of DAEs" tutorial is to demonstrate ModelingToolkit's utility for automatic index reduction, which enables the classic implicit Runge-Kutta solvers such as Rodas5 to solve this DAE. Let's see if that still helps in this context here.

First, modelingtoolkitize the problem:

traced_sys = modelingtoolkitize(pendulum_prob)

\[ \begin{align} + [0.9905247799261719, -0.22537551686665527, -0.13733412142246734, -1.6259873619138987, -4.039397159711885] + [0.9860478635340651, -0.30073881653778056, -0.1664620427217003, -1.781877069088252, -4.897256533167963] + [0.9803344626636025, -0.3882225173693333, -0.197343211472203, -1.9287504836515301, -5.808941760007055] + [0.9735476385359553, -0.48374584567209383, -0.22848413502658996, -2.0608598494557318, -6.7314773610507235] + [0.9660215638346084, -0.5821903599880927, -0.25846150309914806, -2.174808501343366, -7.623501407344294] + [0.9581252462517632, -0.6791831555269633, -0.2863494943871072, -2.2702581674292546, -8.457578518679306] + [0.9512746785824608, -0.7592332348196944, -0.30834478058129816, -2.338838284797042, -9.118979833194777] + [0.945362978941864, -0.8257295231332944, -0.32601971896089704, -2.3898180239992906, -9.65281669795631] + [0.9445580045387526, -0.8346223374638927, -0.32834459959791545, -2.3962584682163732, -9.723238999078445]

Nope! The EK1 is able to solve the index-3 DAE directly. Pretty cool!

plot(sol)
Example block output

Is index-reduction still worth it?

The point of the "Automatic Index Reduction of DAEs" tutorial is to demonstrate ModelingToolkit's utility for automatic index reduction, which enables the classic implicit Runge-Kutta solvers such as Rodas5 to solve this DAE. Let's see if that still helps in this context here.

First, modelingtoolkitize the problem:

traced_sys = modelingtoolkitize(pendulum_prob)

\[ \begin{align} \frac{\mathrm{d} x_1\left( t \right)}{\mathrm{d}t} &= x_2\left( t \right) \\ \frac{\mathrm{d} x_2\left( t \right)}{\mathrm{d}t} &= x_1\left( t \right) x_5\left( t \right) \\ \frac{\mathrm{d} x_3\left( t \right)}{\mathrm{d}t} &= x_4\left( t \right) \\ @@ -131,7 +131,7 @@ sol3_final_error = norm(sol3.u[end] - truesol.u[end]) sol3_f_evals = sol3.stats.nf @info "Results" sol1_final_error sol1_f_evals sol3_final_error sol3_f_evals

┌ Info: Results
-│   sol1_final_error = 0.010930749098469738
+│   sol1_final_error = 0.010926370202213426
 │   sol1_f_evals = 1085
-│   sol3_final_error = 0.06122741667364596
-└   sol3_f_evals = 1409

The error for the index-1 DAE solve is much lower. Thus it seems that, even if the index-3 DAE could also be solved directly, index lowering might still be beneficial when solving DAEs with the EK1!

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Vol. 151 of Proceedings of Machine Learning Research, edited by G. Camps-Valls, F. J. Ruiz and I. Valera (PMLR, 28–30 Mar 2022); pp. 10015–10027.
+│ sol3_final_error = 0.06182444138049775 +└ sol3_f_evals = 1377

The error for the index-1 DAE solve is much lower. Thus it seems that, even if the index-3 DAE could also be solved directly, index lowering might still be beneficial when solving DAEs with the EK1!

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Vol. 151 of Proceedings of Machine Learning Research, edited by G. Camps-Valls, F. J. Ruiz and I. Valera (PMLR, 28–30 Mar 2022); pp. 10015–10027.
diff --git a/dev/tutorials/dynamical_odes/5a45fb09.svg b/dev/tutorials/dynamical_odes/2353f20a.svg similarity index 99% rename from dev/tutorials/dynamical_odes/5a45fb09.svg rename to dev/tutorials/dynamical_odes/2353f20a.svg index 808255410..b2b17ceae 100644 --- a/dev/tutorials/dynamical_odes/5a45fb09.svg +++ b/dev/tutorials/dynamical_odes/2353f20a.svg @@ -39,7 +39,7 @@ - + diff --git a/dev/tutorials/dynamical_odes/239c19ed.svg b/dev/tutorials/dynamical_odes/239c19ed.svg deleted file mode 100644 index 07055745e..000000000 --- a/dev/tutorials/dynamical_odes/239c19ed.svg +++ /dev/null @@ -1,52 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/dynamical_odes/5884f0f2.svg b/dev/tutorials/dynamical_odes/5884f0f2.svg new file mode 100644 index 000000000..ff72fdf48 --- /dev/null +++ b/dev/tutorials/dynamical_odes/5884f0f2.svg @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dynamical_odes/08393b82.svg b/dev/tutorials/dynamical_odes/64959796.svg similarity index 86% rename from dev/tutorials/dynamical_odes/08393b82.svg rename to dev/tutorials/dynamical_odes/64959796.svg index 8481923cc..2b094fcbf 100644 --- a/dev/tutorials/dynamical_odes/08393b82.svg +++ b/dev/tutorials/dynamical_odes/64959796.svg @@ -39,7 +39,7 @@ - + diff --git a/dev/tutorials/dynamical_odes/afc0b3df.svg b/dev/tutorials/dynamical_odes/afc0b3df.svg deleted file mode 100644 index 8e93c1500..000000000 --- a/dev/tutorials/dynamical_odes/afc0b3df.svg +++ /dev/null @@ -1,50 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/dynamical_odes/ecd44a8c.svg b/dev/tutorials/dynamical_odes/ecd44a8c.svg new file mode 100644 index 000000000..3a18e71cb --- /dev/null +++ b/dev/tutorials/dynamical_odes/ecd44a8c.svg @@ -0,0 +1,52 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dynamical_odes/index.html b/dev/tutorials/dynamical_odes/index.html index f61e029c7..c73044331 100644 --- a/dev/tutorials/dynamical_odes/index.html +++ b/dev/tutorials/dynamical_odes/index.html @@ -19,13 +19,13 @@ tspan = (0.0, 100.0) prob = ODEProblem(Hénon_Heiles, [du0; u0], tspan) sol = solve(prob, EK1()); -plot(sol, idxs=(3, 4)) # where `idxs=(3,4)` is used to plot x agains yExample block output

Solving the second-order ODE directly

Instead of first transforming the problem, we can also solve it directly as a second-order ODE, by defining it as a SecondOrderODEProblem.

Note

The SecondOrderODEProblem type is not defined in ProbNumDiffEq.jl but is provided by SciMLBase.jl. For more information, check out the DifferentialEquations.jl documentation on Dynamical, Hamiltonian and 2nd Order ODE Problems.

function Hénon_Heiles2(ddu, du, u, p, t)
+plot(sol, idxs=(3, 4)) # where `idxs=(3,4)` is used to plot x agains y
Example block output

Solving the second-order ODE directly

Instead of first transforming the problem, we can also solve it directly as a second-order ODE, by defining it as a SecondOrderODEProblem.

Note

The SecondOrderODEProblem type is not defined in ProbNumDiffEq.jl but is provided by SciMLBase.jl. For more information, check out the DifferentialEquations.jl documentation on Dynamical, Hamiltonian and 2nd Order ODE Problems.

function Hénon_Heiles2(ddu, du, u, p, t)
     ddu[1] = -u[1] - 2 * u[1] * u[2]
     ddu[2] = u[2]^2 - u[2] - u[1]^2
 end
 prob2 = SecondOrderODEProblem(Hénon_Heiles2, du0, u0, tspan)
 sol2 = solve(prob2, EK1());
-plot(sol2, idxs=(3, 4))
Example block output

Benchmark: Solving second order ODEs is faster

Solving second-order ODEs is not just a matter of convenience - in fact, SciMLBase's SecondOrderODEProblem is neatly designed in such a way that all the classic solvers from OrdinaryDiffEq.jl can handle it by solving the corresponding first-order ODE. But, transforming the ODE to first order increases the dimensionality of the problem, and comes therefore at increased computational cost; this also motivates classic specialized solvers for second-order ODEs.

The probabilistic numerical solvers from ProbNumDiffEq.jl have the same internal state representation for first and second order ODEs; all that changes is the measurement model [1]. As a result, we can use the EK1 both for first and second order ODEs, but it automatically specializes on the latter to provide a 2x performance boost:

julia> using BenchmarkTools
+plot(sol2, idxs=(3, 4))
Example block output

Benchmark: Solving second order ODEs is faster

Solving second-order ODEs is not just a matter of convenience - in fact, SciMLBase's SecondOrderODEProblem is neatly designed in such a way that all the classic solvers from OrdinaryDiffEq.jl can handle it by solving the corresponding first-order ODE. But, transforming the ODE to first order increases the dimensionality of the problem, and comes therefore at increased computational cost; this also motivates classic specialized solvers for second-order ODEs.

The probabilistic numerical solvers from ProbNumDiffEq.jl have the same internal state representation for first and second order ODEs; all that changes is the measurement model [1]. As a result, we can use the EK1 both for first and second order ODEs, but it automatically specializes on the latter to provide a 2x performance boost:

julia> using BenchmarkTools
 
 julia> @btime solve(prob, EK1(order=3), adaptive=false, dt=1e-2);
   317.336 ms (140561 allocations: 140.41 MiB)
@@ -39,7 +39,7 @@
 E(dx, dy, x, y) = PotentialEnergy(x, y) + KineticEnergy(dx, dy)
 E(u) = E(u...); # convenient shorthand
E (generic function with 2 methods)

Let's have a look at how the total energy changes over time when we numerically simulate the Hénon-Heiles model over a long period of time: Standard solve

longprob = remake(prob2, tspan=(0.0, 1e3))
 longsol = solve(longprob, EK1(smooth=false), dense=false)
-plot(longsol.t, E.(longsol.u))
Example block output

It visibly loses energy over time, from an initial 0.12967 to a final 0.12899. Let's fix this to get a physically more meaningful solution.

Energy preservation with the ManifoldUpdate callback

In the language of ODE filters, preserving energy over time amounts to just another measurement model [1]. The most convenient way of updating on this additional zero measurement with ProbNumDiffEq.jl is with the ManifoldUpdate callback.

Note

The ManifoldUpdate callback can be thought of a probabilistic counterpart to the ManifoldProjection callback provided by DiffEqCallbacks.jl.

To do so, first define a (vector-valued) residual function, here chosen to be the difference between the current energy and the initial energy, and build a ManifoldUpdate callback

residual(u) = [E(u) - E(du0..., u0...)]
+plot(longsol.t, E.(longsol.u))
Example block output

It visibly loses energy over time, from an initial 0.12967 to a final 0.12899. Let's fix this to get a physically more meaningful solution.

Energy preservation with the ManifoldUpdate callback

In the language of ODE filters, preserving energy over time amounts to just another measurement model [1]. The most convenient way of updating on this additional zero measurement with ProbNumDiffEq.jl is with the ManifoldUpdate callback.

Note

The ManifoldUpdate callback can be thought of a probabilistic counterpart to the ManifoldProjection callback provided by DiffEqCallbacks.jl.

To do so, first define a (vector-valued) residual function, here chosen to be the difference between the current energy and the initial energy, and build a ManifoldUpdate callback

residual(u) = [E(u) - E(du0..., u0...)]
 cb = ManifoldUpdate(residual)
DiscreteCallback{ProbNumDiffEq.var"#condition#140", ProbNumDiffEq.var"#affect!#141"{Int64, Float64, Float64, typeof(Main.residual)}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Nothing}(ProbNumDiffEq.var"#condition#140"(), ProbNumDiffEq.var"#affect!#141"{Int64, Float64, Float64, typeof(Main.residual)}(100, 1.0e-25, 1.0e-15, Main.residual), SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing)

Then, solve the ODE with this callback

longsol_preserving = solve(longprob, EK1(smooth=false), dense=false, callback=cb)
 plot(longsol.t, E.(longsol.u))
-plot!(longsol_preserving.t, E.(longsol_preserving.u))
Example block output

Voilà! With the ManifoldUpdate callback we could preserve the energy over time and obtain a more truthful probabilistic numerical long-term simulation of the Hénon-Heiles model.

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Vol. 151 of Proceedings of Machine Learning Research, edited by G. Camps-Valls, F. J. Ruiz and I. Valera (PMLR, 28–30 Mar 2022); pp. 10015–10027.
+plot!(longsol_preserving.t, E.(longsol_preserving.u))Example block output

Voilà! With the ManifoldUpdate callback we could preserve the energy over time and obtain a more truthful probabilistic numerical long-term simulation of the Hénon-Heiles model.

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Vol. 151 of Proceedings of Machine Learning Research, edited by G. Camps-Valls, F. J. Ruiz and I. Valera (PMLR, 28–30 Mar 2022); pp. 10015–10027.
diff --git a/dev/tutorials/exponential_integrators/index.html b/dev/tutorials/exponential_integrators/index.html index f55e4139b..080b5dd5e 100644 --- a/dev/tutorials/exponential_integrators/index.html +++ b/dev/tutorials/exponential_integrators/index.html @@ -32,4 +32,4 @@ plot(ylims=(0.3, 1.05)) plot!(ref, color=:black, linestyle=:dash, label="Reference") -plot!(sol_expek1, denseplot=false, marker=:o, markersize=2, label="EK1 + IOUP")Example block output

References

[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).
+plot!(sol_expek1, denseplot=false, marker=:o, markersize=2, label="EK1 + IOUP")Example block output

References

[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators. In: Thirty-seventh Conference on Neural Information Processing Systems (2023).
diff --git a/dev/tutorials/getting_started/7dc5a135.svg b/dev/tutorials/getting_started/7dc5a135.svg deleted file mode 100644 index 2e46238cf..000000000 --- a/dev/tutorials/getting_started/7dc5a135.svg +++ /dev/null @@ -1,50 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/getting_started/b2752536.svg b/dev/tutorials/getting_started/b2752536.svg new file mode 100644 index 000000000..eb01c28fc --- /dev/null +++ b/dev/tutorials/getting_started/b2752536.svg @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/getting_started/b98fcdaf.svg b/dev/tutorials/getting_started/b98fcdaf.svg new file mode 100644 index 000000000..a74cb8378 --- /dev/null +++ b/dev/tutorials/getting_started/b98fcdaf.svg @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/getting_started/bd9f291c.svg b/dev/tutorials/getting_started/bd9f291c.svg deleted file mode 100644 index e0c4932e9..000000000 --- a/dev/tutorials/getting_started/bd9f291c.svg +++ /dev/null @@ -1,54 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/getting_started/c8d5dd68.svg b/dev/tutorials/getting_started/c8d5dd68.svg new file mode 100644 index 000000000..5a6d1a639 --- /dev/null +++ b/dev/tutorials/getting_started/c8d5dd68.svg @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/getting_started/d81735df.svg b/dev/tutorials/getting_started/d81735df.svg deleted file mode 100644 index 48094a3e3..000000000 --- a/dev/tutorials/getting_started/d81735df.svg +++ /dev/null @@ -1,50 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/getting_started/index.html b/dev/tutorials/getting_started/index.html index b6cca961a..37ea1d8f5 100644 --- a/dev/tutorials/getting_started/index.html +++ b/dev/tutorials/getting_started/index.html @@ -15,7 +15,7 @@ prob = ODEProblem(fitz, u0, tspan, p) sol = solve(prob, EK1()) -plot(sol)Example block output

Step 1: Define the problem

First, import ProbNumDiffEq.jl

using ProbNumDiffEq

Then, set up the ODEProblem exactly as you would in DifferentialEquations.jl. Define the vector field

function fitz(du, u, p, t)
+plot(sol)
Example block output

Step 1: Define the problem

First, import ProbNumDiffEq.jl

using ProbNumDiffEq

Then, set up the ODEProblem exactly as you would in DifferentialEquations.jl. Define the vector field

function fitz(du, u, p, t)
     a, b, c = p
     du[1] = c * (u[1] - u[1]^3 / 3 + u[2])
     du[2] = -(1 / c) * (u[1] - a - b * u[2])
@@ -24,77 +24,77 @@
 p = (0.2, 0.2, 3.0)
 prob = ODEProblem(fitz, u0, tspan, p)

Step 2: Solve the problem

To solve the ODE we just use DifferentialEquations.jl's solve interface, together with one of the algorithms implemented in this package. For now, let's use EK1:

sol = solve(prob, EK1())
retcode: Success
 Interpolation: ODE Filter Posterior
-t: 267-element Vector{Float64}:
+t: 266-element Vector{Float64}:
   0.0
   0.021276864853851562
-  0.05530062522770621
-  0.09069833974374658
-  0.13926827781702378
-  0.18486670656185625
-  0.24179366397911012
-  0.29051816272523345
-  0.3490769287679599
-  0.39571551758025386
+  0.055300624266158616
+  0.09069833746250563
+  0.13926827172051978
+  0.18486670143008446
+  0.24179365923629548
+  0.29051816084594195
+  0.349076930314825
+  0.3957155180134495
   ⋮
- 19.52620471671648
- 19.57118993725592
- 19.618145281565827
- 19.669683561096146
- 19.72660662373306
- 19.789272758467693
- 19.860208126620186
- 19.94055004101715
+ 19.50477232090273
+ 19.56373664631041
+ 19.601537671872666
+ 19.65357250448659
+ 19.709230758519457
+ 19.767685129234636
+ 19.83797606872198
+ 19.913564325034223
  20.0
-u: 267-element Vector{Vector{Float64}}:
+u: 266-element Vector{Vector{Float64}}:
  [-1.0, 1.0]
- [-0.9783978986607919, 1.0098599972789337]
- [-0.9424079091803095, 1.0253304289278917]
- [-0.9028542905882088, 1.0410170752594412]
- [-0.8445349615757097, 1.0618117231749666]
- [-0.7847703259504635, 1.0804970332988892]
- [-0.7018976697192769, 1.1025562993298472]
- [-0.6220689201271815, 1.1201790771323805]
- [-0.5125846235780994, 1.139594841674336]
- [-0.412302426024602, 1.1534759756505175]
+ [-0.9783978986604488, 1.0098599972786104]
+ [-0.942407910238602, 1.0253304285109315]
+ [-0.9028542931496027, 1.0410170741923244]
+ [-0.8445349697989268, 1.061811721267787]
+ [-0.7847703324141683, 1.0804970303389172]
+ [-0.7018976822120261, 1.1025563047040523]
+ [-0.6220689234312479, 1.1201790754906849]
+ [-0.5125846363437299, 1.1395948721268157]
+ [-0.4123024317571745, 1.1534759911322578]
  ⋮
- [2.0826740294562938, 0.9061101524201137]
- [2.078744560756923, 0.8805858429639286]
- [2.0732021275130457, 0.853936833241539]
- [2.0660462633712045, 0.8247002716083807]
- [2.057329129463748, 0.7924433159689349]
- [2.0471229128144666, 0.7569883434597815]
- [2.0350861534147566, 0.7169394922557]
- [2.0210436519484087, 0.6717007174964914]
- [2.0104405118668827, 0.6383145073764079]

That's it! we just computed a probabilistic numerical ODE solution!

Step 3: Analyze the solution

Let's plot the result with Plots.jl.

using Plots
-plot(sol)
Example block output

Looks good! Looks like the EK1 managed to solve the Fitzhugh-Nagumo problem quite well.

Tip

To learn more about plotting ODE solutions, check out the plotting tutorial for DifferentialEquations.jl + Plots.jl provided here. Most of that works exactly as expected with ProbNumDiffEq.jl.

Plot the probabilistic error estimates

The plot above looks like a standard ODE solution – but it's not! The numerical errors are just so small that we can't see them in the plot, and the probabilistic error estimates are too. We can visualize them by plotting the errors and error estimates directly:

using OrdinaryDiffEq, Statistics
+ [2.083872897192169, 0.918268165750397]
+ [2.0795125000065657, 0.8848226145152666]
+ [2.075297316215094, 0.8633685451854366]
+ [2.068379627055637, 0.833843863276972]
+ [2.0600615084210836, 0.8022918066068052]
+ [2.050696723446599, 0.7692012836461666]
+ [2.0389038967500825, 0.729487596596834]
+ [2.0258028691294294, 0.686887353466895]
+ [2.0104415505713544, 0.6383210413653303]

That's it! we just computed a probabilistic numerical ODE solution!

Step 3: Analyze the solution

Let's plot the result with Plots.jl.

using Plots
+plot(sol)
Example block output

Looks good! Looks like the EK1 managed to solve the Fitzhugh-Nagumo problem quite well.

Tip

To learn more about plotting ODE solutions, check out the plotting tutorial for DifferentialEquations.jl + Plots.jl provided here. Most of that works exactly as expected with ProbNumDiffEq.jl.

Plot the probabilistic error estimates

The plot above looks like a standard ODE solution – but it's not! The numerical errors are just so small that we can't see them in the plot, and the probabilistic error estimates are too. We can visualize them by plotting the errors and error estimates directly:

using OrdinaryDiffEq, Statistics
 reference = solve(prob, Vern9(), abstol=1e-9, reltol=1e-9, saveat=sol.t)
 errors = reduce(hcat, mean.(sol.pu) .- reference.u)'
 error_estimates = reduce(hcat, std.(sol.pu))'
 plot(sol.t, errors, label="error", color=[1 2], xlabel="t", ylabel="err")
 plot!(sol.t, zero(errors), ribbon=3error_estimates, label="error estimate",
-      color=[1 2], alpha=0.2)
Example block output

More about the ProbabilisticODESolution

The solution object returned by ProbNumDiffEq.jl mostly behaves just like any other ODESolution in DifferentialEquations.jl – with some added uncertainties and related functionality on top. The ProbabilisticODESolution can be indexed with

julia> sol.u[1]2-element Vector{Float64}:
+      color=[1 2], alpha=0.2)
Example block output

More about the ProbabilisticODESolution

The solution object returned by ProbNumDiffEq.jl mostly behaves just like any other ODESolution in DifferentialEquations.jl – with some added uncertainties and related functionality on top. The ProbabilisticODESolution can be indexed with

julia> sol.u[1]2-element Vector{Float64}:
  -1.0
   1.0
julia> sol.u[end]2-element Vector{Float64}:
- 2.0104405118668827
- 0.6383145073764079
julia> sol.t[end]20.0

But since sol is a probabilistic numerical ODE solution, it contains a Gaussian distributions over solution values. The marginals of this posterior are stored in sol.pu:

julia> sol.pu[end]Gaussian{Vector{Float64},PSDMatrix{Float64, Matrix{Float64}}}(
-  μ=[2.0104405118668827, 0.6383145073764079],
-  Σ=2x2 PSDMatrix{Float64, Matrix{Float64}}; R=[4.819328866957359e-5 0.00015008628165071627; 3.994225962938706e-5 0.00012478474684357176; -1.3359035909208665e-5 -4.029546230185406e-5; -2.361011408675205e-7 -8.04110127435485e-6; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
+ 2.0104415505713544
+ 0.6383210413653303
julia> sol.t[end]20.0

But since sol is a probabilistic numerical ODE solution, it contains a Gaussian distributions over solution values. The marginals of this posterior are stored in sol.pu:

julia> sol.pu[end]Gaussian{Vector{Float64},PSDMatrix{Float64, Matrix{Float64}}}(
+  μ=[2.0104415505713544, 0.6383210413653303],
+  Σ=2x2 PSDMatrix{Float64, Matrix{Float64}}; R=[2.2911755902881306e-5 7.030881226868274e-5; 4.7881113296024815e-5 0.000149914216586101; -3.848896497691175e-5 -0.00011288553520532613; -2.214055515418541e-6 -3.3121471754153486e-5; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
 )

You can compute means, covariances, and standard deviations via Statistics.jl:

julia> using Statistics
julia> mean(sol.pu[5])2-element Vector{Float64}:
- -0.8445349615757097
-  1.0618117231749666
julia> cov(sol.pu[5])2x2 PSDMatrix{Float64, Matrix{Float64}}
+ -0.8445349697989268
+  1.061811721267787
julia> cov(sol.pu[5])2x2 PSDMatrix{Float64, Matrix{Float64}}
  Right square root: R=8×2 Matrix{Float64}:
- -2.8014e-6  -1.24815e-7
-  0.0        -2.74772e-6
-  0.0         0.0
-  0.0         0.0
-  0.0         0.0
-  0.0         0.0
-  0.0         0.0
-  0.0         0.0
julia> std(sol.pu[5])2-element Vector{Float64}:
- 2.8014003609822313e-6
- 2.7505533924668052e-6

Dense output

Probabilistic numerical ODE solvers approximate the posterior distribution

\[p \Big( y(t) ~\big|~ y(0) = y_0, \{ \dot{y}(t_i) = f_\theta(y(t_i), t_i) \} \Big),\]

which describes a posterior not just for the discrete steps but for any $t$ in the continuous space $t \in [0, T]$; in classic ODE solvers, this is also known as "interpolation" or "dense output". The probabilistic solutions returned by our solvers can be interpolated as usual by treating them as functions, but they return Gaussian distributions

julia> sol(0.45)Gaussian{Vector{Float64},PSDMatrix{Float64, Matrix{Float64}}}(
-  μ=[-0.2773821283086794, 1.1675659430627088],
-  Σ=2x2 PSDMatrix{Float64, Matrix{Float64}}; R=[-3.2083544077891136e-5 -4.7899698453730076e-6; 0.0 -2.4206415622162586e-5; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
+ 2.8014e-6  1.24815e-7
+ 0.0        2.74772e-6
+ 0.0        0.0
+ 0.0        0.0
+ 0.0        0.0
+ 0.0        0.0
+ 0.0        0.0
+ 0.0        0.0
julia> std(sol.pu[5])2-element Vector{Float64}:
+ 2.8013996013753136e-6
+ 2.750552649946573e-6

Dense output

Probabilistic numerical ODE solvers approximate the posterior distribution

\[p \Big( y(t) ~\big|~ y(0) = y_0, \{ \dot{y}(t_i) = f_\theta(y(t_i), t_i) \} \Big),\]

which describes a posterior not just for the discrete steps but for any $t$ in the continuous space $t \in [0, T]$; in classic ODE solvers, this is also known as "interpolation" or "dense output". The probabilistic solutions returned by our solvers can be interpolated as usual by treating them as functions, but they return Gaussian distributions

julia> sol(0.45)Gaussian{Vector{Float64},PSDMatrix{Float64, Matrix{Float64}}}(
+  μ=[-0.2773821468415606, 1.1675659958232383],
+  Σ=2x2 PSDMatrix{Float64, Matrix{Float64}}; R=[3.2083546826003285e-5 4.789976253577815e-6; 0.0 2.4206406783751565e-5; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
 )
julia> mean(sol(0.45))2-element Vector{Float64}:
- -0.2773821283086794
-  1.1675659430627088

Next steps

Check out one of the other tutorials:

+ -0.2773821468415606 + 1.1675659958232383

Next steps

Check out one of the other tutorials:

diff --git a/dev/tutorials/ode_parameter_inference/48a1b67f.svg b/dev/tutorials/ode_parameter_inference/48a1b67f.svg deleted file mode 100644 index e18bdb653..000000000 --- a/dev/tutorials/ode_parameter_inference/48a1b67f.svg +++ /dev/null @@ -1,87 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/ode_parameter_inference/689abaa5.svg b/dev/tutorials/ode_parameter_inference/689abaa5.svg deleted file mode 100644 index 870ed3c8d..000000000 --- a/dev/tutorials/ode_parameter_inference/689abaa5.svg +++ /dev/null @@ -1,92 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/ode_parameter_inference/6f790629.svg b/dev/tutorials/ode_parameter_inference/6f790629.svg new file mode 100644 index 000000000..710a374c1 --- /dev/null +++ b/dev/tutorials/ode_parameter_inference/6f790629.svg @@ -0,0 +1,87 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/ode_parameter_inference/74285845.svg b/dev/tutorials/ode_parameter_inference/74285845.svg new file mode 100644 index 000000000..ca460b2aa --- /dev/null +++ b/dev/tutorials/ode_parameter_inference/74285845.svg @@ -0,0 +1,92 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/ode_parameter_inference/9a3cb630.svg b/dev/tutorials/ode_parameter_inference/9a3cb630.svg new file mode 100644 index 000000000..d63652d82 --- /dev/null +++ b/dev/tutorials/ode_parameter_inference/9a3cb630.svg @@ -0,0 +1,92 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/ode_parameter_inference/ea82dc47.svg b/dev/tutorials/ode_parameter_inference/ea82dc47.svg deleted file mode 100644 index a39790d88..000000000 --- a/dev/tutorials/ode_parameter_inference/ea82dc47.svg +++ /dev/null @@ -1,92 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/ode_parameter_inference/index.html b/dev/tutorials/ode_parameter_inference/index.html index 2310bee54..5cac37be9 100644 --- a/dev/tutorials/ode_parameter_inference/index.html +++ b/dev/tutorials/ode_parameter_inference/index.html @@ -26,18 +26,18 @@ odedata = [H*true_sol(t) .+ σ * randn() for t in times] plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""]) -scatter!(times, stack(odedata)', color=1, label=["Noisy Data" ""])Example block output

Our goal is then to recover the true parameter p (and thus also the true trajectory plotted above) the noisy data.

Computing the negative log-likelihood

To do parameter inference - be it maximum-likelihod, maximum a posteriori, or full Bayesian inference with MCMC - we need to evaluate the likelihood of given a parameter estimate $\theta_\text{est}$, which corresponds to the probability of the data under the trajectory returned by the ODE solver

θ_est = (0.1, 0.1, 2.0)
+scatter!(times, stack(odedata)',  color=1, label=["Noisy Data" ""])
Example block output

Our goal is then to recover the true parameter p (and thus also the true trajectory plotted above) the noisy data.

Computing the negative log-likelihood

To do parameter inference - be it maximum-likelihod, maximum a posteriori, or full Bayesian inference with MCMC - we need to evaluate the likelihood of given a parameter estimate $\theta_\text{est}$, which corresponds to the probability of the data under the trajectory returned by the ODE solver

θ_est = (0.1, 0.1, 2.0)
 prob = remake(true_prob, p=θ_est)
 plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""])
 scatter!(times, stack(odedata)', color=1, label=["Noisy Data" ""])
 sol = solve(prob, EK1(), adaptive=false, dt=1e-1)
-plot!(sol, color=2, label=["Numerical solution for θ_est" ""])
Example block output

This quantity can be computed in multiple ways; see Data Likelihoods. Here we use ProbNumDiffEq.DataLikelihoods.fenrir_data_loglik:

using ProbNumDiffEq.DataLikelihoods
+plot!(sol, color=2, label=["Numerical solution for θ_est" ""])
Example block output

This quantity can be computed in multiple ways; see Data Likelihoods. Here we use ProbNumDiffEq.DataLikelihoods.fenrir_data_loglik:

using ProbNumDiffEq.DataLikelihoods
 
 data = (t=times, u=odedata)
 nll = -fenrir_data_loglik(
     prob, EK1(smooth=true);
     data, observation_noise_cov=σ^2, observation_matrix=H,
-    adaptive=false, dt=1e-1)
6576.165525142743

This is the negative marginal log-likelihood of the parameter θ_est. You can use it as any other NLL: Optimize it to compute maximum-likelihood estimates or MAPs, or plug it into MCMC to sample from the posterior. In our paper [3] we compute MLEs by pairing Fenrir with Optimization.jl and ForwardDiff.jl. Let's quickly explore how to do this next.

Maximum-likelihood parameter inference

To compute a maximum-likelihood estimate (MLE), we just need to maximize $\theta \to p(\mathcal{D} \mid \theta)$ - that is, minimize the nll from above. We use Optimization.jl for this. First, define a loss function and create an OptimizationProblem

using Optimization, OptimizationOptimJL
+    adaptive=false, dt=1e-1)
6561.828515917825

This is the negative marginal log-likelihood of the parameter θ_est. You can use it as any other NLL: Optimize it to compute maximum-likelihood estimates or MAPs, or plug it into MCMC to sample from the posterior. In our paper [3] we compute MLEs by pairing Fenrir with Optimization.jl and ForwardDiff.jl. Let's quickly explore how to do this next.

Maximum-likelihood parameter inference

To compute a maximum-likelihood estimate (MLE), we just need to maximize $\theta \to p(\mathcal{D} \mid \theta)$ - that is, minimize the nll from above. We use Optimization.jl for this. First, define a loss function and create an OptimizationProblem

using Optimization, OptimizationOptimJL
 
 function loss(x, _)
     ode_params = x[begin:end-1]
@@ -61,9 +61,9 @@
  2.0
  1.0

Then, just solve it! Here we use LBFGS:

optsol = solve(optprob, LBFGS())
 p_mle = optsol.u[1:3]
3-element Vector{Float64}:
- 0.20068944932293867
- 0.2011221541489642
- 3.000075251147067

Success! The computed MLE is quite close to the true parameter which we used to generate the data. As a final step, let's plot the true solution, the data, and the result of the MLE:

plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""])
+ 0.20573039697989073
+ 0.2980971657217813
+ 2.9110574908388878

Success! The computed MLE is quite close to the true parameter which we used to generate the data. As a final step, let's plot the true solution, the data, and the result of the MLE:

plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""])
 scatter!(times, stack(odedata)', color=1, label=["Noisy Data" ""])
 mle_sol = solve(remake(true_prob, p=p_mle), EK1())
-plot!(mle_sol, color=3, label=["MLE-parameter Solution" ""])
Example block output

Looks good!

API Documentation

For more details, see the API documentation of ProbNumDiffEq.DataLikelihoods at Data Likelihoods.

References

[3]
F. Tronarp, N. Bosch and P. Hennig. Fenrir: Physics-Enhanced Regression for Initial Value Problems. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 21776–21794.
[10]
+plot!(mle_sol, color=3, label=["MLE-parameter Solution" ""])Example block output

Looks good!

API Documentation

For more details, see the API documentation of ProbNumDiffEq.DataLikelihoods at Data Likelihoods.

References

[3]
F. Tronarp, N. Bosch and P. Hennig. Fenrir: Physics-Enhanced Regression for Initial Value Problems. In: Proceedings of the 39th International Conference on Machine Learning, Vol. 162 of Proceedings of Machine Learning Research, edited by K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu and S. Sabato (PMLR, 17–23 Jul 2022); pp. 21776–21794.
[10]