Multiplicative Latent Force Model (MLFM)

• MLFM is a framework that integrates mechanistic ODEs with latent Gaussian processes to model curved trajectories under geometric constraints.
• It employs multiplicative interactions and Lie algebra structures to ensure state evolution on manifolds such as rotations and quaternions.
• Advanced inference strategies like adaptive gradient matching and mixture models enable robust estimation under both dense and sparse sampling regimes.

The Multiplicative Latent Force Model (MLFM) is an extension of the latent force model framework designed to embed intricate geometric constraints and manifold-valued dynamics within hybrid state-space models. MLFM couples mechanistic ordinary differential equations (ODEs) with multiplicative interactions between the state and latent Gaussian process (GP) dynamics, enabling the principled modeling of curved or group-structured trajectories—but introducing significant inference challenges absent in classical, additive LFM settings (Tait et al., 2018, Tait et al., 2018).

1. Model Definition and Formalism

The MLFM describes a $K$-dimensional latent state $x(t) \in \mathbb{R}^K$ evolving according to the non-autonomous linear ODE

$\frac{d x(t)}{dt} = A(t)\,x(t),$

where the time-varying coefficient matrix $A(t)$ is itself a linear function of latent Gaussian processes: $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$ with each $A_r \in \mathbb{R}^{K \times K}$. The latent forces $g_r(t)$ are independent GPs, typically with squared-exponential (RBF) kernels

$$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$

This model generalizes the additive LFM, where $A(t)$ would be replaced by $D$ (diagonal) and an additive $x(t) \in \mathbb{R}^K$0 term, resulting in solutions where $x(t) \in \mathbb{R}^K$1 remains jointly Gaussian with the GPs. In contrast, the multiplicative form forces the dependency between $x(t) \in \mathbb{R}^K$2 and $x(t) \in \mathbb{R}^K$3 to be nonlinear due to the interaction of $x(t) \in \mathbb{R}^K$4 and $x(t) \in \mathbb{R}^K$5 (Tait et al., 2018, Tait et al., 2018).

Often, the structure matrices $x(t) \in \mathbb{R}^K$6 are further expanded as

$x(t) \in \mathbb{R}^K$7

where the “connection” matrix $x(t) \in \mathbb{R}^K$8 combines a fixed matrix basis $x(t) \in \mathbb{R}^K$9, often drawn from a Lie algebra $\frac{d x(t)}{dt} = A(t)\,x(t),$0. This guarantees that $\frac{d x(t)}{dt} = A(t)\,x(t),$1 for all $\frac{d x(t)}{dt} = A(t)\,x(t),$2, so the solution map (the fundamental solution) $\frac{d x(t)}{dt} = A(t)\,x(t),$3 lies in the associated Lie group $\frac{d x(t)}{dt} = A(t)\,x(t),$4, constraining the trajectories to a specific manifold or group orbit (Tait et al., 2018).

2. Geometric Implications and Trajectory Structure

The intrinsic geometry of MLFM trajectories is determined by the algebraic structure of the $\frac{d x(t)}{dt} = A(t)\,x(t),$5 matrices. When these belong to a Lie algebra $\frac{d x(t)}{dt} = A(t)\,x(t),$6 (e.g., $\frac{d x(t)}{dt} = A(t)\,x(t),$7 for rotation), the solution

$\frac{d x(t)}{dt} = A(t)\,x(t),$8

ensures $\frac{d x(t)}{dt} = A(t)\,x(t),$9 evolves within the Lie group $A(t)$0. This property builds geometric constraints such as norm preservation or group actions directly into the model. The trajectories do not generally constitute a jointly Gaussian process over $A(t)$1, and their structure can encode rotations (Kubo oscillators), quaternions ($A(t)$2 for orientation), and other non-Euclidean evolutions. By contrast, in standard LFM, $A(t)$3 describes a linear-Gaussian process in Euclidean space (Tait et al., 2018, Tait et al., 2018).

3. Inference: Successive Approximations and Mixtures

Inference in MLFMs is complex due to the loss of closed-form solutions for the marginal likelihood of $A(t)$4 given observations and the nonlinear dependence on the latent GPs. Multiple tractable approximate inference strategies have been proposed:

a. Successive Approximations (Picard Iteration)

The solution of the ODE can be represented as the Picard integral: $A(t)$5 Discretizing this integral and applying $A(t)$6 iterations yields a block-matrix mapping $A(t)$7, which is then used to construct the $A(t)$8th-order approximation to the trajectory and its covariance via recursion: $A(t)$9 Likelihood evaluation and optimization over $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$0 (and optionally hyperparameters) proceed by recursively propagating this covariance and forming Laplace-style approximations for the joint posterior (Tait et al., 2018).

b. Mixtures of Successive Approximations (MLFM–MixSA)

To improve local accuracy for long or non-uniformly spaced trajectories, expansions are mixed over multiple anchor points (“local centers”). At each anchor $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$1, a local Picard expansion of mean $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$2 is built; the conditional density is then

$A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$3

and the full likelihood is a mixture

$A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$4

An EM algorithm alternates between estimating responsibilities and maximizing the likelihood (Tait et al., 2018).

4. Adaptive Gradient Matching

As an alternative, the adaptive gradient matching strategy (MLFM–AG) operates by enforcing consistency between an interpolating GP prior placed on the state path $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$5 and the ODE-implied derivatives:

• Each coordinate $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$6 receives a GP prior, resulting in a conditional Gaussian law for the derivative vector $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$7.
• A “regression expert” likelihood for the ODE model error is specified as a Gaussian about the modelled derivative $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$8.
• The product of the GP and regression expert yields a joint Gaussian approximating the likelihood of the observed path.
• In the linear MLFM case (where $A(t) = A_0 + \sum_{r=1}^R A_r\,g_r(t),$9 is affine in $A_r \in \mathbb{R}^{K \times K}$0 with fixed $A_r \in \mathbb{R}^{K \times K}$1 and $A_r \in \mathbb{R}^{K \times K}$2), the full conditionals for $A_r \in \mathbb{R}^{K \times K}$3, $A_r \in \mathbb{R}^{K \times K}$4, and $A_r \in \mathbb{R}^{K \times K}$5 given the others are Gaussian, permitting efficient Gibbs or mean-field updates. For each $A_r \in \mathbb{R}^{K \times K}$6, the joint posterior is explicitly computable; the same applies for $A_r \in \mathbb{R}^{K \times K}$7 (Tait et al., 2018).

5. Empirical Evaluation and Benchmarking

Quantitative assessment of MLFM inference approaches has been conducted using both simulated and real-world data:

• Kubo Oscillator ($A_r \in \mathbb{R}^{K \times K}$8): Simulation studies compare MLFM–AG and MLFM–MixSA by $A_r \in \mathbb{R}^{K \times K}$9-distance between the true $g_r(t)$0 and its MAP estimate. MLFM–AG achieves high accuracy at small sampling intervals, but deteriorates as data become sparse. MLFM–MixSA retains accuracy for sparse sampling, particularly with increased mixture centers and expansion order (Tait et al., 2018, Tait et al., 2018).
• Dynamics on $g_r(t)$1: MLFM–MixSA achieves lower reconstruction errors than MLFM–AG on sparsely observed rotational data, with improvements saturating at expansion orders $g_r(t)$2.
• Human Motion-Capture (MOCAP) Data: MLFM is applied to the modeling of joint trajectories represented as quaternions ($g_r(t)$3). Marginal MLFMs with $g_r(t)$4 latent forces suffice for individual joints. Coupled ("product MLFM") models with shared latent GPs ($g_r(t)$5) yield superior predictive accuracy, reflecting the pooling of information across coupled manifold-valued paths.

Benchmarking results for the Wasserstein-2 distance between the true and approximate GP posterior for the Kubo oscillator, as a function of truncation order $g_r(t)$6, interval $g_r(t)$7, and sample step $g_r(t)$8, show a monotonic drop in error with increasing $g_r(t)$9, and better performance for shorter time horizons (Tait et al., 2018).

6. Computational Properties and Practical Considerations

A comparison of inference schemes is summarized below:

Method	Accuracy	Computational Complexity	Scalability & Suitability
MLFM–AG	High (dense data)	$$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$0 for GP plus low-dimensional updates	Good for large $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$1, $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$2, well-sampled data
MLFM–MixSA	High (sparse data)	$$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$3, larger for more anchors, $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$4	Challenging for high $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$5 or long series

MLFM–AG is fast and effective when observation grids are dense enough that GP interpolation is accurate, but loses geometric fidelity when the state manifold is strongly curved (e.g., spheres, rotations) or sampling is sparse. MLFM–MixSA, using localized Taylor–Picard expansion mixtures, remains robust in these challenging regimes at the price of increased computation (Tait et al., 2018). In practice, MLFM–AG is recommended as a default, with MLFM–MixSA employed in high-curvature or sparsely sampled scenarios.

7. Summary of Implementation Approach

Implementation proceeds via:

1. Discretization of observation times $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$6, quadrature nodes and weights.
2. Initialization of prior distributions for $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$7, GP on $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$8.
3. Construction of block-matrix operators $$g_r(t) \sim \mathcal{GP}\bigl(0,\,k_r(t,t';\phi_r)\bigr), \quad\text{where}\quad k_r(t, t') = \sigma_r^2 \exp\left(-\frac{(t-t')^2}{2\ell_r^2}\right).$$9 for each $A(t)$0.
4. Recursion of covariance propagation for $A(t)$1 steps, yielding marginal likelihood via sequenced Gaussian updates or direct optimization.
5. Parameter and hyperparameter updating by (gradient) optimization or EM, as dictated by chosen approximation (Tait et al., 2018).

These procedure details are summarized in implementation pseudocode in (Tait et al., 2018), supporting Laplace-style MAP optimization.

---

The MLFM provides a flexible, geometrically expressive class of latent dynamical models for manifold-valued systems, at the cost of nontrivial but tractable approximate inference. Empirical results confirm suitability for rotational dynamics, coupled manifold states, and scenarios requiring the integration of differential geometry with modern machine learning techniques (Tait et al., 2018, Tait et al., 2018).