ODEPrior¶

class probnum.diffeq.ODEPrior(driftmat, dispmat, ordint, spatialdim, precond_step=1.0)[source]¶

Bases: probnum.filtsmooth.statespace.continuous.linearsdemodel.LTISDEModel

Prior dynamic model for ODE filtering and smoothing.

An ODE prior is an continuous LTI state space model with attributes:

order of integration \(q\)
spatial dimension of the underlying ODE
projection to \(X_i(t)\) (the \((i-1)\)-th derivative estimate)
A preconditioner \(P\) (see below)

Instead of the LTI SDE

\[d X(t) = [F X(t) + u] dt + L dB(t)\]

the prior for the ODE Dynamics is given by

\[dX(t) = P F P^{-1} X(t) dt + P L dB(t)\]

where \(P\) is a preconditioner matrix ensuring stability of the iterations. Note that ODE priors do not have a drift term. By default, we choose \(P\) to be the matrix that maps to filtering iteration to the Nordsieck vector,

\[P = \text{diag }(h^{-q}, h^{-q+1}, ..., 1).\]

Here, \(h\) is some expected average step size. Note that we ignored the factorials in this matrix. Our setting makes it easy to recover “no preconditioning” by choosing \(h=1\).

If no expected step size is available we choose \(h=1.0\). This recovers \(P=I_{d(q+1)}\), hence no preconditioning.
For fixed step size algorithms this quantity \(h\) is easy to choose
For adaptive steps it is a bit more involved.

Since it doesn’t have to be exact, any more or less appropriate choice will do well. The main effect of this preconditioning is that the predictive covariances inside each filter iteration are well-conditioned: for IBM(\(q\)) priors, the condition number of the predictive covariances only depends on order of integration \(q\) and not on the step size anymore. Nb: this only holds if all required derivatives of the RHS vector field of the ODE are specified: None for IBM(1), Jacobian of \(f\) for IBM(2), Hessian of \(f\) for IBM(3). If this is not the case the preconditioner still helps but is not as powerful anymore.

Without preconditioning they can be numerically singular for small steps and higher order methods which especially makes smoothing algorithms unstable.

Another advantage of this preconditioning is that the smallest value that appears inside the algorithm is \(h^{q}\) (with preconditioning) instead of \(h^{2q+1}\) (without preconditioning).

The matrices \(F, u, L\) are the usual matrices for IBM(\(q\)), IOUP(\(q\)) or Matern(\(q+1/2\)) processes. As always, \(B(t)\) is s-dimensional Brownian motion with unit diffusion matrix \(Q\).

Attributes Summary

`diffusionmatrix`	Evaluates Q.
`dispersionmatrix`
`driftmatrix`
`force`
`inverse_preconditioner`	Convenience property to return the readily-computed inverse preconditioner without having to remember abbreviations.
`ndim`	Spatial dimension (utility attribute).
`preconditioner`	Convenience property to return the readily-computed preconditioner without having to remember abbreviations.

Methods Summary

`chapmankolmogorov`(start, stop, step, …)	Solves Chapman-Kolmogorov equation from start to stop via step.
`dispersion`(time, state, **kwargs)	Evaluates l(t, x(t)) = L(t).
`drift`(time, state, **kwargs)	Evaluates f(t, x(t)) = F(t) x(t) + u(t).
`jacobian`(time, state, **kwargs)	maps t -> F(t)
`precond2nordsieck`(step)	Computes preconditioner inspired by Nordsieck.
`proj2coord`(coord)	Projection matrix to \(i\)-th coordinates.
`sample`(start, stop, step, initstate, **kwargs)	Samples from `initstate` at `start` to `stop` with stepsize `step`.

Attributes Documentation

diffusionmatrix¶: Evaluates Q.

dispersionmatrix¶

driftmatrix¶

force¶

inverse_preconditioner¶

Convenience property to return the readily-computed inverse preconditioner without having to remember abbreviations.

Returns:	Inverse preconditioner matrix \(P^{-1}\)
Return type:	np.ndarray, shape=(d(q+1), d(q+1))

ndim¶: Spatial dimension (utility attribute).

preconditioner¶

Convenience property to return the readily-computed preconditioner without having to remember abbreviations.

Returns:	Preconditioner matrix \(P\)
Return type:	np.ndarray, shape=(d(q+1), d(q+1))

Methods Documentation

chapmankolmogorov(start, stop, step, randvar, **kwargs)¶

Solves Chapman-Kolmogorov equation from start to stop via step.

For LTISDEs, there is a closed form solutions to the ODE for mean and kernels (see super().chapmankolmogorov(…)). We exploit this for [(stop - start)/step] steps.

References

Eq. (8) in http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.390.380&rep=rep1&type=pdf and Eq. 6.41 and Eq. 6.42 in Applied SDEs.

dispersion(time, state, **kwargs)¶: Evaluates l(t, x(t)) = L(t).

drift(time, state, **kwargs)¶: Evaluates f(t, x(t)) = F(t) x(t) + u(t).

jacobian(time, state, **kwargs)¶: maps t -> F(t)

precond2nordsieck(step)[source]¶

Computes preconditioner inspired by Nordsieck.

Computes the matrix \(P\) given by

\[P = I_d \otimes diag (1, h, h^2, ..., h^q)\]

as well as its inverse \(P^{-1}\).

Parameters:	step (float) – Step size \(h\) used for preconditioning. If \(h\) is so small that \(h^q! < 10^{-15}\), it is being set to \(h = (\cdot 10^{-15})^{1/q}\).
Returns:	precond (np.ndarray, shape=(d(q+1), d(q+1))) – Preconditioner matrix \(P\). invprecond (np.ndarray, shape=(d(q+1), d(q+1))) – Inverse preconditioner matrix \(P^{-1}\).

proj2coord(coord)[source]¶

Projection matrix to \(i\)-th coordinates.

Computes the matrix

\[H_i = \left[ I_d \otimes e_i \right] P^{-1},\]

where \(e_i\) is the \(i\)-th unit vector, that projects to the \(i\)-th coordinate of a vector. If the ODE is multidimensional, it projects to each of the \(i\)-th coordinates of each ODE dimension.

Parameters:	coord (int) – Coordinate index \(i\) which to project to. Expected to be in range \(0 \leq i \leq q + 1\).
Returns:	Projection matrix \(H_i\).
Return type:	np.ndarray, shape=(d, d*(q+1))

sample(start, stop, step, initstate, **kwargs)¶

Samples from initstate at start to stop with stepsize step.

Start, stop and step lead to a np.arange-like interface. Returns a single element at the end of the time, not the entire array!