kipe.roukf module

The reduced-order unscented Kalman filter (ROUKF): its linear algebra.

This module is the filter itself, on plain arrays: it knows no forward solver, no measurements and no files, so it can be tested on its own. Running it over a study (propagating the sigma points with the solver, computing the innovations, outer iterations, output) is done in kipe.estimation, which calls initial_state(), sample() and update().

The method of [MC11], in the formulation and notation of [NAB26]. The filter estimates the parameters \(\theta \in \mathbb{R}^p\); the model state \(\chi\) (x in the code; \(n\) entries on this MPI rank) is corrected along, but has no uncertainty of its own at the start (reduced order: \(L_\chi^0 = 0\)). The uncertainty is carried by the factors \(L_\chi\) (\(n \times p\)), \(L_\theta\) and \(U\) (both \(p \times p\)), with the parameter covariance

\[P_\theta = L_\theta U^{-1} L_\theta^T.\]

The sigma points are the estimates perturbed along a fixed SigmaPointStencil \(I_{(*)}\), whose \(r\) columns \(I_{(i)} \in \mathbb{R}^p\) all have the weight \(\alpha = 1/r\), with the weighted mean

\[E_\alpha(Y_{(*)}) = \alpha \sum_{i=1}^r Y_{(i)}\]

and \(Y_{(*)}\) the matrix collecting the vectors \(Y_{(i)}\) column-wise.

The estimate and its uncertainty factors form the FilterState. One assimilation step consists of four stages. Propagation and innovations involve the forward solver and the measurements and are done in kipe.estimation; sampling and correction are done here. [NAB26] counts the means of the propagated sigma points to the prediction; here they start the correction, since they need no solver.

  1. Sampling (sample()): with \(C\) the lower Cholesky factor of \(U^{-1}\), the sigma points are

    \[\hat{\chi}_{(i)} = \hat{\chi}_+ + L_\chi C I_{(i)}, \qquad \hat{\theta}_{(i)} = \hat{\theta}_+ + L_\theta C I_{(i)}.\]

    Some publications, e.g., [Ber12] and [NB22], write \(C^T\) here, with \(C\) a Cholesky factor of \(U^{-1}\): that holds for the upper factor (\(C^T C = U^{-1}\), e.g., MATLAB’s chol), not for the lower one used here (\(C C^T = U^{-1}\), numpy’s cholesky).

  2. Propagation: each sigma point is propagated with the forward model to the next measurement time (kipe.estimation), the parameters stay constant.

  3. Innovation: per sigma point, with the spatial sampler and measurement model (kipe.estimation),

    \[\Gamma_{(i)} = Z - \mathcal{H}(\hat{\chi}_{(i)}).\]

  4. Correction (update()): the a priori estimates (mean() of the propagated sigma points), their sensitivity(), the gain() and the corrected (a posteriori) estimates,

    \[ \begin{align}\begin{aligned}\hat{\chi}_- = E_\alpha(\hat{\chi}_{(*)}), \qquad \hat{\theta}_- = E_\alpha(\hat{\theta}_{(*)}),\\L_\chi = \alpha \hat{\chi}_{(*)} I_{(*)}^T, \qquad L_\theta = \alpha \hat{\theta}_{(*)} I_{(*)}^T, \qquad L_\Gamma = \alpha \Gamma_{(*)} I_{(*)}^T,\\U = P_\alpha + L_\Gamma^T W^{-1} L_\Gamma, \qquad M = U^{-1} L_\Gamma^T W^{-1} E_\alpha(\Gamma_{(*)}),\\\hat{\chi}_+ = \hat{\chi}_- - L_\chi M, \qquad \hat{\theta}_+ = \hat{\theta}_- - L_\theta M.\end{aligned}\end{align} \]

Here the innovations are passed weighted with the noise, \(W^{-1/2} \Gamma_{(i)}\), so \(W\) does not appear in the code. States are rank-local vectors: the state fields, concatenated. All functions are pure; only gain() communicates (two small sums over the ranks).

[MC11] (1,2,3,4)

Moireau, P., & Chapelle, D. (2011). Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems. ESAIM: Control, Optimisation and Calculus of Variations, 17(2), 380-405.

[Ber12] (1,2,3)

Bertoglio, C. (2012). Forward and Inverse Problems in Fluid-Structure Interaction. Application to Hemodynamics (PhD thesis). Université Pierre et Marie Curie - Paris VI.

[NB22] (1,2)

Nolte, D., & Bertoglio, C. (2022). Inverse problems in blood flow modeling: A review. International Journal for Numerical Methods in Biomedical Engineering, 38(8), e3613. https://doi.org/10.1002/cnm.3613

[NAB26] (1,2,3)

Nolte, D., Aróstica, R., & Bertoglio, C. Parameter estimation and identifiability in biventricular cardiac fluid-structure interaction from displacement and flow measurements. Submitted to Computer Methods in Applied Mechanics and Engineering.

class kipe.roukf.SigmaPointStencil(points: NDArray_f64, alpha: float)[source]

Bases: object

The stencil \(I_{(*)}\) along which the sigma points are sampled, and its weight.

points: NDArray_f64

The stencil points \(I_{(i)}\), one per column, shape (p, r).

alpha: float

The weight \(\alpha = 1/r\) of every stencil point.

p_alpha() → NDArray_f64[source]

Return \(P_\alpha = \alpha I_{(*)} I_{(*)}^T\).

Returns:

\(P_\alpha\), shape (p, p)

class kipe.roukf.FilterState(x: NDArray_f64, theta: NDArray_f64, L_x: NDArray_f64, L_theta: NDArray_f64, U_inv: NDArray_f64)[source]

Bases: object

The estimate and the factors that carry its uncertainty.

x: NDArray_f64

State estimate \(\hat{\chi}_+\), shape (n,).

theta: NDArray_f64

Parameter estimate \(\hat{\theta}_+\), shape (p,).

L_x: NDArray_f64

State sensitivity \(L_\chi\), shape (n, p).

L_theta: NDArray_f64

Parameter sensitivity \(L_\theta\), shape (p, p).

U_inv: NDArray_f64

\(U^{-1}\), shape (p, p); zero for unique.

covariance() → NDArray_f64[source]

Return the parameter covariance \(P_\theta = L_\theta U^{-1} L_\theta^T\).

Returns:

\(P_\theta\), shape (p, p)

kipe.roukf.sigma_point_stencil(sigma_points: SigmaPoints, p: int) → SigmaPointStencil[source]

Construct the sigma-point stencil of the given kind, after [MC11].

  • simplex: \(p + 1\) points on a regular simplex; \(P_\alpha = I\).

  • canonical: \(2p\) points \(\pm\sqrt{p}\,e_j\); \(P_\alpha = I\).

  • unique: one point at the origin: no spread, hence no correction. A sanity check that runs the whole pipeline with fixed parameters.

Note

[MC11] also list star, the canonical points and the origin. Its weights are not given; for \(P_\alpha = I\) and weights summing to one, the origin must have weight zero, so it equals canonical. Equal weights \(\alpha=1/(2p + 1)\) would give \(P_\alpha = \frac{2p}{2p + 1} I\): the sigma points would represent too small a covariance, at every step.

Parameters:

  • sigma_points – kind of stencil

  • p – number of parameters

Returns:

the stencil

kipe.roukf.initial_state(x0: NDArray_f64, theta0: NDArray_f64, stddev_theta: NDArray_f64, stencil: SigmaPointStencil) → FilterState[source]

Return the initial filter state.

As in [NAB26], the initial factors are

\[L_\chi^0 = 0, \qquad L_\theta^0 = \sqrt{P_\theta^0}, \qquad U^0 = P_\alpha,\]

with \(\sqrt{P_\theta^0}\) the Cholesky factor of the diagonal initial covariance, so that \(P_\theta^0 = L_\theta^0 (U^0)^{-1} (L_\theta^0)^T\) for stencils with \(P_\alpha = I\). This differs from [MC11] and [Ber12], which use the inconsistent \((U^0)^{-1} = (P_\theta^0)^{-1/2}\), \(L_\theta^0 = I\).

Parameters:

  • x0 – initial state, shape (n,); known, without uncertainty of its own

  • theta0 – initial parameter estimate, shape (p,)

  • stddev_theta – initial standard deviations of \(\theta\), shape (p,)

  • stencil – the sigma-point stencil

Returns:

the initial filter state

kipe.roukf.sample(state: FilterState, stencil: SigmaPointStencil) → tuple[NDArray_f64, NDArray_f64][source]

Sample the sigma points around the estimate.

\[\hat{\chi}_{(i)} = \hat{\chi}_+ + L_\chi C I_{(i)}, \qquad \hat{\theta}_{(i)} = \hat{\theta}_+ + L_\theta C I_{(i)},\]

with \(C\) the lower Cholesky factor of \(U^{-1}\).

Some publications, e.g., [Ber12] and [NB22], write \(C^T\) here, with \(C\) a Cholesky factor of \(U^{-1}\): that holds for the upper factor (\(C^T C = U^{-1}\), e.g., MATLAB’s chol), not for the lower one used here (\(C C^T = U^{-1}\), numpy’s cholesky).

Parameters:

  • state – the filter state, around whose estimate the sigma points are sampled

  • stencil – the sigma-point stencil

Returns:

  • sigma-point states, one per column, shape (n, r)

  • sigma-point parameters, one per column, shape (p, r)

kipe.roukf.mean(ensemble: NDArray_f64, stencil: SigmaPointStencil) → NDArray_f64[source]

Return the weighted mean \(E_\alpha(Y_{(*)}) = \alpha \sum_i Y_{(i)}\).

Parameters:

  • ensemble – one sigma point per column, each of size k (e.g., n for states, p for parameters), shape (k, r)

  • stencil – the sigma-point stencil

Returns:

the mean, shape (k,)

kipe.roukf.sensitivity(ensemble: NDArray_f64, stencil: SigmaPointStencil) → NDArray_f64[source]

Return the sensitivity \(\alpha Y_{(*)} I_{(*)}^T\) of the sigma points.

Parameters:

  • ensemble – one sigma point per column, each of size k (e.g., n for states, p for parameters), shape (k, r)

  • stencil – the sigma-point stencil

Returns:

the sensitivity, shape (k, p)

kipe.roukf.gain(Gamma: NDArray_f64, stencil: SigmaPointStencil, comm: mpi4py.MPI.Comm) → tuple[NDArray_f64, NDArray_f64][source]

Compute \(U^{-1}\) and the correction factor \(M\) from the innovations.

\[L_\Gamma = \alpha \Gamma_{(*)} I_{(*)}^T, \qquad U = P_\alpha + L_\Gamma^T L_\Gamma, \qquad M = U^{-1} L_\Gamma^T E_\alpha(\Gamma_{(*)}),\]

for innovations weighted with the noise (so \(W = I\)). The products \(L_\Gamma^T L_\Gamma\) and \(L_\Gamma^T E_\alpha(\Gamma_{(*)})\) are sums over the innovation entries, summed over the ranks; each entry must be held by exactly one rank (see kipe.sampling.SpatialSampler.global_indices).

Parameters:

  • Gamma – innovations, weighted with the noise, one sigma point per column; the m innovation entries held by this rank (possibly none), shape (m, r)

  • stencil – the sigma-point stencil

  • comm – communicator over which the innovation entries are distributed

Returns:

  • \(U^{-1}\), shape (p, p); zero for unique

  • \(M\), shape (p,)

kipe.roukf.update(x_sigma: NDArray_f64, theta_sigma: NDArray_f64, Gamma: NDArray_f64, stencil: SigmaPointStencil, comm: mpi4py.MPI.Comm) → FilterState[source]

Return the corrected filter state from the propagated sigma points and their innovations.

The correction stage of the assimilation step: the a priori estimates (mean()), the sensitivities (sensitivity()), the gain() and the corrected estimates.

Parameters:

  • x_sigma – the propagated sigma-point states, one per column, shape (n, r)

  • theta_sigma – the sigma-point parameters, one per column, shape (p, r)

  • Gamma – their innovations, weighted with the noise, see gain()

  • stencil – the sigma-point stencil

  • comm – communicator over which the state and innovation entries are distributed

Returns:

the corrected filter state