ODEPrior¶

class probnum.diffeq.ODEPrior(driftmat: numpy.ndarray, dispmat: numpy.ndarray, ordint: int, spatialdim: int, precond_step: float = 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.
`dimension`	Spatial dimension (utility attribute).
`dispersionmatrix`
`driftmatrix`
`force`
`inverse_preconditioner`	Convenience property to return the readily-computed inverse preconditioner without having to remember abbreviations.
`preconditioner`	Convenience property to return the readily-computed preconditioner without having to remember abbreviations.

Methods Summary

`__call__`(arr_or_rv, RandomVariable], start, …)	Transition a random variable or a realization of one.
`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.
`transition_realization`(real, start, stop, …)	Transition a realization of a random variable from time \(t\) to time \(t+\Delta t\).
`transition_rv`(rv, start, stop, **kwargs)	Transition a random variable from time \(t\) to time \(t+\Delta t\).

Attributes Documentation

diffusionmatrix¶: Evaluates Q.

dimension¶: Spatial dimension (utility attribute).

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))

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

__call__(arr_or_rv: Union[numpy.ndarray, RandomVariable], start: float = None, stop: float = None, **kwargs) -> ('RandomVariable', typing.Dict)¶

Transition a random variable or a realization of one.

The input is either interpreted as a random variable or as a realization. Accordingly, the respective methods are called: transition_realization() or transition_rv().

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: float) -> (<class 'numpy.ndarray'>, <class 'numpy.ndarray'>)[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: int) → numpy.ndarray[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))

transition_realization(real, start, stop, **kwargs)¶

Transition a realization of a random variable from time \(t\) to time \(t+\Delta t\).

For random variable \(x_t\), it returns the random variable defined by

\[x_{t + \Delta t} \sim p(x_{t + \Delta t} | x_t = r) .\]

This is different to transition_rv() which computes the parametrization of \(x_{t + \Delta t}\) based on the parametrization of \(x_t\).

Nb: Think of transition as a verb, i.e. this method “transitions” a realization of a random variable.

Parameters:

real – Realization of the random variable.
start – Starting point \(t\).
stop – End point \(t + \Delta t\).

Returns:

RandomVariable – Random variable, describing the state at time \(t + \Delta t\) based on realization at time \(t\).
dict – Additional information in form of a dictionary, for instance the cross-covariance in the prediction step, access to which is useful in smoothing.