Many modern tasks consist of prior, mechanistic knowledge that we want to use and hold on to, but we know is incomplete.

Among other things, it motivates the framework of universal differential equations [Rackauckas et al., 2020].

A hybrid model is the statistical realization of this pattern, combining a known and unknown mechanism.

For example, in continuous time:

\(\def\V#1{\boldsymbol{\mathbf{#1}}} \mathrm{d}\V{x}_t = \left(f_{\mathrm{known}}(\V{x}_t) + f_{\mathrm{unknown}}^{\V\theta}(\V{x}_t)\right)\mathrm{d}t + \left(g_{\mathrm{known}}(\V{x}_t) + g_{\mathrm{unknown}}^{\V\theta}(\V{x}_t)\right)\mathrm{d}\V\beta_t.\)

Writing down a hybrid model is easy. Making it behave is not:

• How do I write \(f_{\mathrm{unknown}}\) so that it does not overpower my prior knowledge \(f_{\mathrm{known}}\)?
• How do I set priors on such a thing?
• How do I do proper uncertainty quantification?

The solution I present is not fancy — from a Bayesian perspective it is the straightforward "right" thing to do.

The point is to show how it may work, and that it deals naturally with things like UQ, partial observations, noisy observations, and arbitrary missingness patterns.

Goal 1. Show a simple Bayesian recipe for hybrid dynamical systems: that it can work, and that you may want to use it.

Goal 2. Show how dynestyx [DW et al., 2026] makes this modeling explicit, efficient, and easy to use.

Goal 3. Convince you that you should work on these :).

We break our recipe into five steps:

1. Write down the known mechanistic component \(f_{\mathrm{known}}, g_{\mathrm{known}}\), and any priors it carries.
2. Put a prior on the unknown components \(f_{\mathrm{unknown}}^{\theta}\) and, where needed, \(g_{\mathrm{unknown}}^{\theta}\).
3. Specify the observation likelihood \(\def\V#1{\boldsymbol{\mathbf{#1}}}p(\V{y}_t \mid \V{x}_t, \V\theta)\) and priors for unknown observation parameters.
4. Choose a method for marginalizing the dynamical system. We use an EnKF, which we find very effective; others are certainly possible.
5. Choose a method for inference on \(\def\V#1{\boldsymbol{\mathbf{#1}}}\V\theta\). We use black-box stochastic variational inference.

This is really one part of a larger Bayesian workflow; today we focus on model building.

We use two generic families throughout.

Sparse dictionary. A small set of low-degree polynomial terms with a Laplace prior to shrink coefficients — analogous to SINDy [Brunton et al., 2016]. When it works, results are interpretable and data-efficient... it does not always work.

HSGP closure. A Gaussian process prior on \(f_{\mathrm{unknown}}\); we'll use a squared-exponential kernel, treated fully Bayesian.

To make this tractable we use the Hilbert space GP approximation [Solin & Särkkä, 2020] — a truncated Karhunen–Loève-like expansion on a bounded domain with Dirichlet boundary conditions. What matters here is that we end up with

\(\def\V#1{\boldsymbol{\mathbf{#1}}}f(\V{x}_t) \approx \V\theta^\top \V\phi(\V{x}_t; \ell, \sigma).\)

The two-scale Lorenz–96 system has \(K = 9\) slow variables, with \(J = 8\) fast variables attached to each slow coordinate:

\[\begin{aligned} \dot X_k &= -X_{k-1}(X_{k-2} - X_{k+1}) - X_k + F + h_x \bar Z_k, \\ \dot Z_{k,j} &= \frac{-Z_{k,j+1}(Z_{k,j+2} - Z_{k,j-1}) - Z_{k,j} + h_y X_k}{\varepsilon}, \qquad \bar Z_k = \frac{1}{J}\sum_{j=1}^J Z_{k,j}, \end{aligned}\]

with circular indexing and

\(\varepsilon = 2^{-5}, \qquad h_x = -0.8, \qquad h_y = 1, \qquad F_\star = 10.\)

The model is chosen so that a Markovian closure is possible [Levine & Stuart, 2022].

To specify a reduced model, we simply omit all observation and knowledge of the fast states. Naïvely, this leaves

\(\dot X_k = -X_{k-1}(X_{k-2} - X_{k+1}) - X_k + F + c(X_k).\)

We simulate from the full system with \(4.0\) time units of burn-in, \(10.0\) of training data, and \(4.0\) of test rollout. Observations arrive every \(0.1\) time units.

The baseline noisily observes every resolved coordinate \(X_1, \ldots, X_9\) with independent Gaussian noise of standard deviation \(\sigma_{y,\star} = 0.15\). All \(72\) fast coordinates remain hidden.

We leave the forcing \(F\) and the observation scale \(\sigma_y\) unknown, with priors

\(F \sim \mathcal{N}(10.0, 2.0^2), \qquad \sigma_y \sim \operatorname{HalfNormal}(0.25).\)

Recovering these interpretable quantities is a primary scientific target.

We also want to learn the closure \(c\), which is what lets us do state estimation. So we measure three things:

• parameter recovery,
• filtered-state accuracy, and
• held-out forecast accuracy.

Naïve Bayes. The fully-reduced model, \(c(x) = 0\).

SDE Relaxation. Model the missing components generically, as a stochastic term:

\(\mathrm{d}X_k = \left[-X_{k-1}(X_{k-2} - X_{k+1}) - X_k + F\right]\mathrm{d}t + q\,\mathrm{d}W_t.\)

Sparse closure. A dictionary in the rescaled state (\(x_k\) divided by \(10\)):

\(c_{\mathrm{sparse}}(x_k) = \beta_1 x_k + \beta_2 x_k^2 + \beta_3 x_k^3, \qquad \beta_j \sim \operatorname{Laplace}(0, 0.5).\)

HSGP closure. Ten squared-exponential basis modes on the scaled domain \([-1.5, 1.5]\), standard-normal basis weights, a half-normal prior on kernel variance, and a log-normal prior on length scale.

An unrestricted intercept in \(c\) would confound \(F\). We rule it out by construction: the dictionary has \(c(0) = 0\), and the HSGP subtracts off \(f(0)\).

We discretize the ODE and use an ensemble Kalman filter to marginalize the states, then SVI with clipped Adam and one Monte Carlo particle per ELBO estimate [Tran et al., 2017].

Both are interpretations wrapped around the same model code — and missing observations need no model change at all. We mark them with NaN.

One benefit of a closure is that we retain good inference and UQ on the interpretable parameters.

In the ODE case, \(F\) is statistically learnable even without a closure: if \(\bar Z\) is zero-mean, it simply averages out. The SDE relaxation, on the other hand, corrupts it.

Both hybrid models recover the sign and shape of the hidden forcing \(h_x \bar Z_k\), with some amount of uncertainty.

So what is the drawback of the mis-specified model?

The SDE relaxation and naive ODE method are also both poor at forecasting.

Alternating dimensions. We always observe \(X_1, X_3, X_5, X_7, X_9\) (and never \(X_2, X_4, X_6, X_8\)). i.e., the states and closure must be recoverable from a subset of components.

Salt-and-pepper. We randomly drop ~55% of observations, across time and component.

Neither our methods nor dynestyx care much either way; each model is refit by exactly the same recipe (and code!).

The hybrid models (and SDE relaxation) are able to perform this task well; the naive model doesn't.

Closure models outperform the baselines, even under missingness, while maintaining meaningfull parameters.

Consider the reaction–diffusion equation on \(x \in (0,1)\) with homogeneous Dirichlet conditions \(u(t,0) = u(t,1) = 0\):

\[\frac{\partial u}{\partial t}(t,x) = \nu\frac{\partial^2 u}{\partial x^2}(t,x) + c_\star(u(t,x)), \qquad c_\star(u) = r_\star u(1-u).\]

We assume that the diffusion term is known, but the logistic growth is hidden. We set \(\nu = 0.01\) and \(r_\star = 1.1\).

Let \(x_i = ih\), exclude the two fixed boundary values from the state, and write \(\def\V#1{\boldsymbol{\mathbf{#1}}}\V{u}(t) = (u(t,x_1), \ldots, u(t,x_m))^\top\).

The centered finite-difference approximation

\(\frac{\partial^2 u}{\partial x^2}(t,x_i) \approx \frac{u(t,x_{i-1}) - 2u(t,x_i) + u(t,x_{i+1})}{h^2}\)

gives the method-of-lines ODE

\(\def\V#1{\boldsymbol{\mathbf{#1}}}\dot{\V{u}}(t) = \nu L_h \V{u}(t) + c(\V{u}(t)), \qquad L_h = \frac{1}{h^2}\operatorname{tridiag}(1, -2, 1).\)

This is now a differential equation we can use in dynestyx!

We use 40 grid points (including boundaries), and observe every interior node is observed every \(0.1\) time units with Gaussian noise of standard deviation \(0.05\).

Eleven fields from \(t = 0\) to \(t = 1\) are used for training; we then forecast from \(t = 1.1\) to \(t = 2.5\).

The candidates mirror the ODE example:

\[\begin{aligned} c_{\mathrm{diffusion}}(u) &= 0, \\ c_{\mathrm{sparse}}(u) &= \beta_1 u + \beta_2 u^2 + \beta_3 u^3, \qquad \beta_j \sim \operatorname{Laplace}(0,1), \\ c_{\mathrm{HSGP}}(u) &= \mathrm{HSGP}(u). \end{aligned}\]

This is quite generous to the "sparse" solution, as the law really is quadratic!

We keep \(\nu\) known, and use the EnKF with SVI once again.

We recover the reactions well; on values during training, the reaction RMSE is \(0.2103\) for diffusion only, \(0.0190\) for the sparse dictionary, and \(0.0305\) for the HSGP.

The filtered states show clearly that, despite informative results, the mis-specified model loses quickly. Hybrid models learn well.

In forecasting using the learned reaction law, our hybrid models work quite well.

1. We presented a simple Bayesian recipe to deal with hybrid models.
2. Two simple patterns: a sparse dictionary closure and GP closure, both work well in our toy examples.
3. With dynestyx, we showed how these models are simple to implement, and deal with things like missingness and differing inference methods.
4. We don't get perfect UQ here! Though we get much better than if we didn't try anything.

1. Our examples where quite simple here; let's get this to work on bigger and harder problems!
2. We also had "easy" fixes to identifiability; how do we deal with this more generally?
3. How far can we take the PDE formulation?
4. What can we do for better calibration?
5. How can we use generative methods (e.g., diffusion models) as part of this framework?