1 Introduction↩︎

Conservation laws are a subclass of partial differential equations (PDEs) with applications ranging from weather prediction [1] and aerodynamics [2] to flood modeling [3] and subsurface contaminant transport [4]. They take the general form \[\frac{\partial u}{\partial t} + \nabla \cdot F(u) = s, \label{eq:conservation-law-general}\tag{1}\] where \(u\) denotes the conserved quantity, \(F(u)\) the (in general nonlinear) flux operator, and \(s\) the source term.

Real-world applications incur uncertainty. In the classical forward problem, \(u\) is the object of interest: Under initial and boundary conditions, it is inferred from the known flux and source terms. Practical applications however often operate in settings that lack information: Sensors measuring relevant quantities are noisy and often only available at sparse locations for cost reasons. Worse, in many applications the source itself is unknown—a contamination event must be localized, or an emission field reconstructed, from sparse measurements of \(u\). These inverse problems admit multiple solutions consistent with the data. Orthogonal to these issues, any solver used to approach these problems may only produce an approximate solution. Together, these factors demand methods that quantify uncertainty – whether epistemic, computational, or otherwise – to help practitioners make informed decisions.

In the machine learning community, neural approaches such as physics-informed neural networks [5] and neural operators [6] have gained popularity for PDE solving. While powerful, they are point estimators by design. Ensembles or Bayesian variants may quantify uncertainty, but at significant computational cost and often with poor calibration [7].

Probabilistic numerical methods frame numerical methods as Bayesian inference, thereby treating uncertainty naturally while matching the accuracy of classical methods in the mean [8]. Recent work has developed such methods for the solution of PDEs by imposing a Gaussian process (GP) prior over the object(s) of interest and conditioning on physical information [9], [10]. While this is principled, GP inference naively scales cubically in the discretization resolution. These approaches are mostly based on collocation, which corresponds to a pointwise discretization of the PDE.

Numerical conservation. Conservation laws derive their name from a defining property: quantities such as mass or energy are neither created nor destroyed. Numerical methods that fail to preserve this can produce unphysical solutions that drift over time. The finite volume method (FVM) [11] enforces conservation by construction, discretizing the integral balance rather than the PDE pointwise — making it the standard in computational fluid dynamics. Expressing FVM as a scalable probabilistic numerical scheme has remained an open problem.

Recent work has begun addressing scalability for GP-based PDE solvers. [12] apply a Vecchia approximation to solve nonlinear PDEs, but via collocation. [13] develop a probabilistic FVM, but their method is restricted to linear PDEs and their low-rank approximation does not yield reliable UQ at large scales.

Our contributions. We present a scalable method for Bayesian inference on nonlinear conservation laws. We formulate the finite volume method as inference under a Gaussian process prior — extending prior work on linear PDEs to nonlinear flux operators through iterative conditioning. To achieve scalability, we extend Vecchia approximations to integral observations and introduce a sparsity-preserving time-stepping scheme. Our framework can solve forward and inverse problems alike, while propagating uncertainties and explicitly embedding prior knowledge.

On a source identification inverse problem (1), our method recovers nonparametric source fields with quantified uncertainty in seconds, outperforming a physics-informed neural network baseline that requires minutes for a less accurate point estimate. 2 develops the high-level method, 3 addresses scalability, and 4 presents experiments.

2 Probabilistic Conservation Schemes↩︎

Our key insight is that the gold standard discretization for conservation laws, namely the finite volume method (FVM), may be expressed as Bayesian inference under a GP prior. We first review how GPs interact with linear functionals (2.1), then recall the classical FVM discretization (2.2), and finally show how to unify them into a coherent Bayesian framework (2.3).

2.1 Gaussian Processes and Linear Functionals↩︎

We model our prior knowledge about the solution function \(u\) in 1 through a Gaussian process (GP) prior with mean function \(\mu\) and kernel function \(k\): \[u \sim \mathcal{GP}(\mu, k). \label{eq:gp-prior}\tag{2}\] In the following, we discuss how to express the simulation of conservation laws as Bayesian inference under this prior.

We leverage that linear transformations of Gaussian processes remain Gaussian: If \(\ell\) is a bounded linear functional acting on \(u\), then \(\ell(u)\) is a Gaussian random variable with \[\ell(u) \sim \mathcal{N}\bigl(\ell \mu, \, \ell k \ell^\prime \bigr), \label{eq:single-functional}\tag{3}\] where \(\ell \mu := \ell[\mu(\cdot)]\) denotes the functional applied to the mean, and \(\ell_1 k \ell_2^\prime := \ell_1[x \mapsto \ell_2[k(x, \cdot)]]\) denotes the functionals applied to each argument of the kernel. More generally, given a collection of linear functionals \(\ell_1, \ldots, \ell_n\), the vector \([\ell_1(u), \ldots, \ell_n(u)]^\top\) is jointly Gaussian: \[\begin{bmatrix} \ell_1(u) \\ \vdots \\ \ell_n(u) \end{bmatrix} \sim \mathcal{N}\left( \begin{bmatrix} \ell_1 \mu \\ \vdots \\ \ell_n \mu \end{bmatrix}, \begin{bmatrix} \ell_1 k \ell_1^\prime & \cdots & \ell_1 k \ell_n^\prime \\ \vdots & \ddots & \vdots \\ \ell_n k \ell_1^\prime & \cdots & \ell_n k \ell_n^\prime \end{bmatrix} \right). \label{eq:joint-functional}\tag{4}\] As a concrete example, consider a GP prior over functions \(u: [0, 1] \to \mathbb{R}\). We can define an evaluation functional \(\ell_1[u] := u(0.5)\) that extracts the function value at the midpoint, and a second-derivative functional \(\ell_2[u] := u''(0.5)\) that extracts the curvature. Applying 4 , the pair \(\boldsymbol{s}_{\text{example}} := (u(0.5), u''(0.5))^\top\) is jointly Gaussian.

Since Gaussians are closed under linear conditioning, we can incorporate observations in closed form. Let \(\boldsymbol{s} \in \mathbb{R}^N\) be Gaussian with mean \(\boldsymbol{\mu}\in \mathbb{R}^N\) and precision matrix \(\boldsymbol{Q} \in \mathbb{R}^{N \times N}\), i.e.\(\boldsymbol{s} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{Q}^{-1})\). Here \(N\) denotes the dimension of the joint state \(\boldsymbol{s}\). Given a matrix \(\boldsymbol{A} \in \mathbb{R}^{M \times N}\) and noisy observations \(\boldsymbol{A} \boldsymbol{s} \sim \mathcal{N}(\boldsymbol{b}, \sigma^2 \boldsymbol{I})\), define \(\boldsymbol{\hat{Q}} := \boldsymbol{Q} + \sigma^{-2} \boldsymbol{A}^\top \boldsymbol{A}\). The posterior \(\boldsymbol{s}^* := (\boldsymbol{s} \mid \boldsymbol{A}\boldsymbol{s} = \boldsymbol{b})\) is \[\boldsymbol{s}^* \sim \mathcal{N}\bigl(\boldsymbol{\hat{Q}}^{-1}(\boldsymbol{Q} \boldsymbol{\mu}+ \sigma^{-2} \boldsymbol{A}^\top \boldsymbol{b}), \;\; \boldsymbol{\hat{Q}}^{-1}\bigr). \label{eq:gaussian-conditioning-general}\tag{5}\]

Returning to our example: We can e.g.set \(\boldsymbol{A} = \begin{bmatrix}0 & 1\end{bmatrix}\) (the row vector selecting \(u''(0.5)\) from \(\boldsymbol{s}_{\text{example}}\)) and \(\boldsymbol{b} = 0\) in 5 to obtain a posterior of \(\boldsymbol{s}_{\text{example}}\) that reflects \(u''(0.5) \approx 0\). We could extend our state by many more, spread-out points at which we evaluate the second derivative (so-called collocation points), which would allow us to approximate a solution to the 1D Laplace equation \[u''(x) = 0. \quad (x \in (0, 1))\] Physics-informed GP regression [10] develops this idea systematically: boundary conditions, PDE discretizations and sensor measurements are all just sources of information (in the sense of [9]), and conditioning yields a posterior that unifies them.

2.2 Finite Volume Discretization↩︎

The Finite Volume Method (FVM) discretizes 1 by integrating over control volumes \(\Omega_1, \dots, \Omega_N\): \[\frac{\partial}{\partial t} \int_{\Omega_i} u \, dV + \int_{\Omega_i} \nabla \cdot F(u) \, dV = \int_{\Omega_i} s \, dV . \label{eq:conservation-law-integrated}\tag{6}\] The divergence theorem yields: \[\frac{\partial}{\partial t} \textcolor{cbblue}{\int_{\Omega_i} u \, dV} + \textcolor{cborange}{\oint_{\partial\Omega_i} F(u) \cdot \mathbf{n} \, dS} = \textcolor{cbteal}{\int_{\Omega_i} s \, dV} . \label{eq:conservation-law}\tag{7}\]

Classical numerical FVM schemes represent the solution in terms of one or multiple nodes per cell, and then approximate these integrals numerically. In the simplest setting, a piecewise constant representation is used. The compute nodes of FVM then are the cell midpoints \(\bar{u}_i\), and the integral over each cell is approximated by the midpoint rule \[\textcolor{cbblue}{\int_{\Omega_i} u \, dV} \approx \lvert \Omega_i \rvert \bar{u}_i .\] The integral of the source function \(s\) as well as the surface integral of the flux \(F\) are similarly approximated through numerical quadrature. The latter requires a scheme to extrapolate or interpolate \(F\) at the cell surface from the compute nodes, for which many different flux schemes exist.

From the discretization of 1 to the use of quadrature and flux schemes, FVM makes various approximations. Next, we will derive a GP-based FVM approach that naturally expresses these approximation errors as uncertainty.

2.3 FVM as Bayesian Inference↩︎

In the following, we express the components in 7 in terms of linear functionals. We then obtain a joint Gaussian prior through the push-forward of GP priors under these linear functionals. Then fulfilling 7 reduces to (non-)linear conditioning of this prior.

Integration is linear. Recall that integration is a linear operation, and thus \(\textcolor{cbblue}{\ell_{vol_i}[u] := \int_{\Omega_i} u \, dV}\) is a linear functional. One might similarly define \(\textcolor{cborange}{\ell_{\text{flux}_i}[u] := \oint_{\partial\Omega_i} F(u) \cdot \mathbf{n} \, dS}\), as done by [13]. But this is only a linear functional if \(F\) is linear, which is not the case for nonlinear PDEs. Thus, the flux term warrants a different approach.

Flux term. Consider as example the flux for the inviscid Burgers’ equation: \(F_{\text{Burgers}}(u) = \frac{1}{2} u^2\). In one spatial dimension, the cells are intervals \(\Omega_i = [b_i, b_{i+1}]\), and so \[\textcolor{cborange}{\oint_{\partial\Omega_i} F(u) \cdot \mathbf{n} \, dS} = F(b_{i+1}) - F(b_i) .\] Thus, we only need linear functionals for the components that make up \(F\) at the cell boundaries. Concretely, for \(F_{\text{Burgers}}\), we may use \(\ell_{\text{boundary}_i}[u] = u(b_i)\).

In 2D and 3D, we need to compute actual line and surface integrals of the flux. We propose a two-stage approach to this. First, we model the flux field with its own prior \(F(u) \sim \mathcal{GP}(\mu_F, k_F)\), treated initially as a quantity independent of \(u\). The deterministic relationship \(F = F(u)\) imposed by the PDE is then enforced through nonlinear conditioning of pointwise evaluations of \(F\) on pointwise evaluations of \(u\) (detailed below), yielding a posterior that respects the actual flux–solution coupling. Afterwards, the flux term is directly a linear functional in \(F(u)\): \(\textcolor{cborange}{\tilde{\ell}_{\text{flux}_i}[F(u)] := \oint_{\partial\Omega_i} F(u) \cdot \mathbf{n} \, dS}\). This amounts to per-cell Bayesian quadrature [14].

Source term. Some commonly presented conservation laws have constant (often zero) source functions, in which case the treatment of the source term in 7 is trivial. If the source function is non-constant (or even unknown), we may employ the same method as for the flux term: Model the source with a GP prior \(s \sim \mathcal{GP}(\mu_s, k_s)\), incorporate point observations of the source (if available), and finally define linear functionals \(\textcolor{cbteal}{\ell_{\text{source}_i}[s] := \int_{\Omega_i} s \, dV}\).

Additional linear functionals are used to embed information on top of the FVM constraints:

• Initial and boundary conditions provide point observations of the solution \(u\).

• Point observations of the flux may be indirectly obtained through observations of \(u\).

• Point observations of the source are either available or the object of interest (in an inverse problem setting).

We add these linear functionals to the joint state. For notational simplicity, we write \(\boldsymbol{\ell}_u : U \rightarrow \mathbb{R}^{N + M_u}\) for the vector-valued functional stacking the \(N\) FVM cell-integral functionals \(\ell_{vol_i}\) together with the \(M_u\) functionals encoding additional information about \(u\); the stacked flux functional \(\boldsymbol{\ell}_F\) (collecting the per-cell \(\tilde{\ell}_{\text{flux}_i}\)) and stacked source functional \(\boldsymbol{\ell}_s\) (collecting \(\ell_{\text{source}_i}\)) are defined analogously.

Putting it all together. We have independent GP priors for the solution, flux, and source. This results in the joint state \[\boldsymbol{s} := \begin{bmatrix} \textcolor{cbblue}{\boldsymbol{\ell}_u[u]} \\ \textcolor{cborange}{\boldsymbol{\ell}_F[F(u)]} \\ \textcolor{cbteal}{\boldsymbol{\ell}_s[s]} \end{bmatrix} \sim \mathcal{N}\left( \begin{bmatrix} \textcolor{cbblue}{\boldsymbol{\mu}_u} \\ \textcolor{cborange}{\boldsymbol{\mu}_F} \\ \textcolor{cbteal}{\boldsymbol{\mu}_s} \end{bmatrix}, \begin{bmatrix} \textcolor{cbblue}{\boldsymbol{\Sigma}_u} & & \\ & \textcolor{cborange}{\boldsymbol{\Sigma}_F} & \\ & & \textcolor{cbteal}{\boldsymbol{\Sigma}_s} \end{bmatrix} \right), \label{eq:joint-state}\tag{8}\] where each block is given by 4 .

Now, we can express our various sources of information in terms of conditioning. Initial and boundary conditions, source observations, and even the FVM constraints in 7 are all linear combinations of state components. We obtain a posterior conditioned on these observations through one application of 5 . Though their priors are independent, this conditioning couples all three quantities in the posterior.

The only possible source of nonlinearity is the flux \(F(u)\): although \(F\) has its own GP prior, the flux–solution map \(F = F(u)\) itself is a deterministic nonlinear relationship that still needs to be enforced. Concretely, the joint state contains pointwise evaluations of both \(F(u)\) and \(u\) at points \(\boldsymbol{x} := (x_1, \dots, x_k)^\top\). We want to enforce a nonlinear relationship \(f(u(\boldsymbol{x}), F(u)(\boldsymbol{x})) = \boldsymbol{0}\). To do so, we find the maximum a posteriori (MAP) estimate of the joint state \(\boldsymbol{s}\) (with the constraint locations \(\boldsymbol{x}\) fixed) under the likelihood \[\boldsymbol{y} \mid \boldsymbol{s} \sim \mathcal{N}(f(u(\boldsymbol{x}), F(u)(\boldsymbol{x})), \sigma^2 \boldsymbol{I}),\] with observations \(\boldsymbol{y} = 0\). Gaussian likelihoods on zero-valued observations of this form have appeared in related probabilistic modeling contexts as virtual likelihoods or virtual observables for softly enforcing physical constraints [15]. We find the MAP estimate via Gauss-Newton optimization and then use a Laplace approximation to obtain a Gaussian posterior.

2.4 Implications↩︎

Unified inference. Our method forms a joint prior over solution, flux and source. This enables us to unify different sources of information (i.e. conservation law, boundary conditions, sensor measurements) in one common language. In addition, we clearly express our structural assumptions about each term through the prior. The result is a framework that handles both forward and inverse problems natively by simply changing what is observed vs.inferred:

The FVM constraints create coupling between the three quantities that propagate the different information sources. This concept of maintaining a joint distribution over the solution trajectory and latent force(s) has been explored in the context of ODE solvers by [16].

Differences to classical FVM. Where classical FVM uses numerical quadrature, we implicitly employ Bayesian quadrature and thus propagate quadrature error to the joint uncertainty. For classical FVM, there are schemes to achieve higher convergence orders for suitable problems. In our framework, as will be demonstrated in 4.3, this may be tweaked effortlessly through the smoothness of the kernel. For classical FVM, there is a whole zoo of flux schemes to arrive at the flux values from the solution function values. In our method, there is no separate flux scheme: The flux term is automatically derived from conditioning the flux prior on observations of the solution.

Computational bottleneck. The key challenge of our method is that it quickly produces a large total state dimension \(N\), and Gaussian inference naively scales as \(\mathcal{O}(N^3)\) compute. Hence, 3 will explore in detail techniques to implement our conceptual framework efficiently.

3 Scalable Inference↩︎

For a scalable implementation of the concepts explored in 2, we propose an approach based on a structured prior (3.1), sparse inference algorithms ([sec:sparse-information] [sec:kl-cholesky]) and a careful treatment of time (3.4).

3.1 A Structured Prior↩︎

We employ a separable prior with tensor product structure. For the spatial dimensions, we use a tensor product of Matérn kernels: \[k_{\text{space}}(\boldsymbol{x}, \boldsymbol{x}') = \prod_{i=1}^{d} k_{\text{Mat}}^{(\nu_i, \lambda_i)}(x_i, x_i'), \label{eq:kernel-space}\tag{9}\] where \(k_{\text{Mat}}^{(\nu, \lambda)}\) denotes the Matérn kernel with smoothness \(\nu\) and lengthscale \(\lambda\). This structure enables closed-form evaluation of integrals [17] and partial derivatives needed for 4 , and allows independent control of smoothness in each coordinate. To obtain closed-form kernel integrals, we use rectangular FVM cells. We deliberately use Matérn rather than e.g.squared-exponential (SE) kernels: the sparse Cholesky approximation in 3.3 relies on the screening effect, which holds for finitely smooth kernels but is known to degrade for the infinite smoothness of the SE kernel [18].

For the temporal dimension, we use an Integrated Wiener Process (IWP) of order \(q\) [19]. The IWP encodes that the \(q\)-th time derivative is Brownian motion, naturally capturing temporal evolution. A side effect is that the time derivative \(\frac{\partial}{\partial t}\) needed for 7 is built into the state representation. Importantly, the IWP is non-stationary, so conservation properties are determined entirely by the FVM constraints. A stationary prior (e.g., temporal Matérn) reverts toward its mean, which cancels the numerical conservation of energy of FVM.

With this, the full spacetime kernel is \[k\bigl((t, \boldsymbol{x}), (t', \boldsymbol{x}')\bigr) = k_{\text{IWP}}^{(q)}(t, t') \cdot k_{\text{space}}(\boldsymbol{x}, \boldsymbol{x}'). \label{eq:kernel-spacetime}\tag{10}\] Note that separability is a property of the prior, not of the posterior. The FVM constraints in 7 couple space and time through conditioning, so the posterior is generally non-separable — reflecting the actual spacetime structure imposed by the PDE. A similar assumption is made in high-dimensional probabilistic ODE solvers [20], where a Kronecker-structured (i.e.separable) prior across state dimensions is coupled through conditioning on ODE residuals.

3.2 Sparse Information Form Inference↩︎

Inference in our framework repeatedly applies 5 . Naively, this involves \(\mathcal{O}(N^2)\) memory and \(\mathcal{O}(N^3)\) compute. Recently, approaches based on sparse precision matrices have successfully alleviated this issue [12], [21].

Assume we discretize \(T\) into \(\tilde{T} := [t_1, \dots, t_{N_t}]\) and use the same spatial state \(s\) across all time points. Then from 10 the precision matrix of the joint distribution across all spacetime components factorizes into \[\boldsymbol{Q_{\tilde{T}, s}} = \boldsymbol{Q_{\tilde{T}}} \otimes \boldsymbol{Q_{s}}, \label{eq:spacetime-precision}\tag{11}\] with temporal precision \(\boldsymbol{Q_{\tilde{T}}}\) and spatial precision \(\boldsymbol{Q_s}\). \(\boldsymbol{Q_{\tilde{T}}}\) is naturally sparse since the IWP is Markovian (7.1). 3.3 covers a sparse approximation to \(\boldsymbol{Q_s}\).

Once we have a sparse approximation to 11 , we use sparse Cholesky factorizations for mean updates, sampling and to compute marginal variances. Memory and compute requirements depend on the sparsity pattern but generally scale much more favorably than dense equivalents. [21] discuss this inference scheme in detail.

3.3 Sparse Cholesky Approximation↩︎

Approximating \(\boldsymbol{Q_s}\) with a sparse matrix \(\boldsymbol{\tilde{Q}_s}\) is known in the literature as a Vecchia approximation [22]. [18] frame the problem as finding a sparse Cholesky factor \(\boldsymbol{\tilde{L}}\) such that the Kullback-Leibler (KL) divergence between \(\mathcal{N}(\boldsymbol{0}, \boldsymbol{Q_s}^{-1})\) and \(\mathcal{N}(\boldsymbol{0}, \boldsymbol{\tilde{L}}^{-\top} \boldsymbol{\tilde{L}}^{-1})\) is minimal. They find a closed-form expression for the optimal \(\boldsymbol{\tilde{L}}\), where each column can be computed in parallel at cost that depends cubically on the sparsity pattern of the column. Thus, there is an inherent tradeoff for the sparsity pattern between accuracy and computational efficiency.

To motivate this idea, one can show that for \(\boldsymbol{x} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{L}^{-\top} \boldsymbol{L}^{-1})\), it holds: \[\mathbb{E}[x_i \mid x_{i+1}, \dots, x_N] = \mu_i - \sum_{j > i} \frac{L_{ji}}{L_{ii}} (x_j - \mu_j) \label{eq:conditional-regression}\tag{12}\] The key insight is that the \(L_{ji}\) (\(j > i\)) form regression coefficients. For strategic node orderings and kernels with suitable properties, this results in so-called screening effects [23]. This means that close points “screen out” distant points, causing their regression coefficients (and thus the \(L_{ji}\)) in 12 to be close to zero. In this setting, Vecchia approximations achieve high accuracy.

Node ordering. Achieving strong screening effects requires a good ordering of the nodes \(x_i\). In a pure regression setting, [18] propose an optimal reverse maximin ordering. In a setting that mixes point evaluations with derivative evaluations, [12] observe that ordering point evaluations before derivative evaluations results in stronger screening effects. Following this intuition of ordering “coarse” nodes before “fine nodes”, we propose to order integrals before point evaluations, and point evaluations before derivatives. Within each category, we employ a reverse maximin ordering. See 7.3 for details and empirical validation.

2 summarizes our approach based on the KL-minimizing sparse Cholesky factorization [18]. Their algorithm includes a parameter \(\rho \in \mathbb{R}^+\) which tunes the sparsity, allowing for an explicit tradeoff between accuracy and computation time. As in the original work, we use supernodes to benefit from BLAS-3 and LAPACK calls. Additionally, due to the block diagonal structure in 8 , we may apply this approximation to each block instead of the entire matrix at once.

3.4 An Alternative Treatment of Time↩︎

Working with a big joint distribution over all spacetime observations via 11 is principled, but expensive. For inference, we need to compute sparse Cholesky factors of updated precision matrices, which induce numerical fill-in that quickly renders this approach prohibitively expensive. In the following, we discuss alternative strategies.

Filtering and smoothing. Probabilistic ODE solvers [24] use Kalman filter and smoother methods to achieve linear-in-time runtime complexity. Unfortunately, the Kalman filter prediction step produces a dense precision matrix even from a sparse prior. As such, it again induces a prohibitive cubic scaling in the spatial dimension.

Marginal Moment-Matching. We propose an alternative that approximates the filtering distributions by matching the marginal moments.

Given a current state distribution \(\boldsymbol{s}_t \sim \mathcal{N}(\boldsymbol{\mu}_t, \boldsymbol{\Sigma}_t)\) and transition distribution \(\boldsymbol{s}_{t+1} \mid \boldsymbol{s}_t \sim \mathcal{N}(\boldsymbol{A} \boldsymbol{s}_t, \boldsymbol{\Sigma}_{\varepsilon})\), the joint distribution is \((\boldsymbol{s}_t, \boldsymbol{s}_{t+1})^\top \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) with mean \(\boldsymbol{\mu} = (\boldsymbol{\mu}_t, \boldsymbol{A}\boldsymbol{\mu}_t)^\top\) and \[\boldsymbol{\Sigma}^{-1} = \begin{pmatrix} \boldsymbol{\Sigma}_t^{-1} + \boldsymbol{A}^\top \boldsymbol{\Sigma}_\varepsilon^{-1} \boldsymbol{A} & -\boldsymbol{A}^\top \boldsymbol{\Sigma}_\varepsilon^{-1} \\ -\boldsymbol{\Sigma}_\varepsilon^{-1} \boldsymbol{A} & \boldsymbol{\Sigma}_\varepsilon^{-1} \end{pmatrix}.\] In our setting, \(\boldsymbol{\Sigma}_t^{-1}\) is the current sparse precision matrix, and \(\boldsymbol{A}\) and \(\boldsymbol{\Sigma}_\varepsilon^{-1}\) result from the state-space representation of 10 together with the sparse approximation \(\boldsymbol{\tilde{Q}_s}\). Thus, \(\boldsymbol{\Sigma}^{-1}\) is sparse as a composition of sparse matrices. This enables us to efficiently perform inference on this two-step joint distribution.

From the conditioned joint distribution, we extract the mean \(\boldsymbol{\tilde{\mu}}_{t+1}\) and marginal variances \(\boldsymbol{\tilde{\sigma}}_{t+1}^2\) of the next state \(\boldsymbol{s}_{t+1}\) via selected inversion [25]. We then construct the approximate filtering distribution for the next time step as \[\boldsymbol{s}_{t+1} \sim \mathcal{N}\left(\boldsymbol{\tilde{\mu}}_{t+1}, \left(\boldsymbol{\tilde{Q}}_s + \mathrm{diag}(\boldsymbol{\tilde{\sigma}}_{t+1}^{-2})\right)^{-1}\right).\] This preserves the marginal means exactly while combining the prior precision structure with the inverse marginal variances. We repeat this procedure for each time step. In practice, we use a Crank-Nicolson discretization for the FVM constraints (7.2); the single Gauss–Newton step used for nonlinear constraints there is validated empirically in 12.

This approach is linear in time and leaves the spatial sparsity pattern unchanged. However, it only provides an approximation to the true filtering distributions; we give a formal error analysis in 13.

3.5 Complexity↩︎

The sparse Cholesky approximation of [18] achieves near-linear complexity: \(\mathcal{O}(N \log(N/\varepsilon)^d)\) space and \(\mathcal{O}(N \log(N/\varepsilon)^{2d})\) time for \(\varepsilon\)-accuracy in \(d\) dimensions. However, this requires access to covariance entries, which are cheap for the prior (kernel evaluations) but expensive for the posterior (requires inverting the precision). Thus, we use the sparse Cholesky approximation to obtain a sparse prior precision, but subsequent conditioning and inference use standard sparse Cholesky factorization [26], where fill-in is governed by the sparsity pattern. For 2D grids with nested dissection ordering, this yields \(\mathcal{O}(N^{3/2})\) factorization time and \(\mathcal{O}(N \log N)\) memory [27]. We validate this scaling empirically in 10 on the source-identification problem of 4.1. We emphasize that the functional ordering of 2 (integrals \(\rightarrow\) evaluations \(\rightarrow\) derivatives) and the nested-dissection ordering serve distinct purposes: the former is a node ordering for the KL-minimizing sparse approximation that produces \(\boldsymbol{\tilde{Q}_s}\), whereas nested dissection is a fill-reducing permutation for factorizing the resulting sparse matrix. The \(\mathcal{O}(N^{3/2})\) bound depends on the bounded-degree planar structure of \(\boldsymbol{\tilde{Q}_s}\) on a 2D mesh, not on the functional ordering used to construct it. Time-stepping via marginal moment-matching (3.4) is linear in \(N_t\), giving overall complexity \(\mathcal{O}(N_t \cdot N_s^{3/2})\) for 2D problems.

4 Experiments↩︎

We evaluate four questions: Can our method solve ill-posed inverse problems with meaningful uncertainty estimates? Does it match classical FVM accuracy on forward problems? What is its convergence order? Does it scale to complex systems?

Our implementation is in Julia.¹ Hardware details and additional experimental parameters are provided in 8.

4.1 Source Identification in Advection-Diffusion↩︎

In this experiment, we model the transport of contaminants in groundwater: A contaminant enters an aquifer and is spread downstream via the flow of water. In steady-state, this may be modeled by the advection-diffusion equation \[\nabla \cdot (\mathbf{v}c - D \nabla c) = s, \label{eq:advection-diffusion}\tag{13}\] where \(c\) is the contaminant concentration, \(s\) is the contaminant source field, \(\mathbf{v}\) the velocity field, and \(D\) the dispersion coefficient. The inverse problem in this setting is source identification: Given measurements of the concentration \(c(x, y)\) at sparse monitoring wells, recover the source field \(s(x, y)\). This problem is ill-posed: Distinct source field configurations can yield similar downstream observations.

1 depicts the concrete setup. The ground-truth source \(s\) (1e) consists of two bell-shaped regions where the contaminant enters. The contaminant then flows towards the right (\(\boldsymbol{v} = (1 \, 0)^{\top}\)) and diffuses (\(D = 0.05\)), resulting in a steady-state concentration \(c\) (1a). This concentration is measured at 12 monitoring wells with additive Gaussian noise (\(\sigma = 0.01\)). For further details, refer to 8.1.

Our method approaches this task through a joint distribution over concentration, source, and flux, as in 8 . As this is a steady-state problem (i.e. \(\frac{\partial c}{\partial t} = 0\)), we use purely spatial priors. The FVM constraints then automatically introduce coupling between the source and the concentration measurements. Posterior samples correspond to inferred possible forms of the source, while marginal means and variances yield a point estimate and confidence intervals.

Baseline. We compare against a physics-informed neural network [5] baseline. It uses a dual-network architecture with separate MLPs for the concentration (4 layers × 128 units) and source field (2 layers × 64 units), trained jointly on PDE residuals, boundary conditions, and observation data. For a fair comparison, we follow state-of-the-art training practices [28], including Fourier feature embeddings and learning rate scheduling.

Point estimates. 1 shows qualitative results of both methods. The posterior mean of our method correctly reconstructs the overall shape of both \(c\) and \(s\) (1c, g) and achieves a source RMSE of 0.44. By contrast, the PINN solution exhibits noticeable artifacts at a source RMSE of 0.76. Our method recovers the source field in 2 seconds, while the PINN baseline requires approximately 5 minutes — a 150x speedup.

Marginal variances. Our method provides uncertainty quantification on top of the point estimate — 1d and 1h show the marginal standard deviations. For both \(s\) and \(c\), the uncertainty shrinks near the measurements. In addition, the uncertainty correctly reflects the fact that our measurements provide more information about \(c\) — indeed, the uncertainty about \(s\) is smeared due to the flow of water. The concentration uncertainty further collapses at the left boundary, which prescribes \(c = 0\) (1d).

Samples. Our method’s posterior not only provides marginal statistics, but also enables us to sample possible solution pairs of \((s, c)\). The top row of 3 shows four samples of \(s\) for our experiment. To further highlight that these samples are indeed physically meaningful, we repeat our experiment with three added observations and sample again using the same seeds. This is shown in the bottom row of 3. The result is that for all samples, the top left corner indeed adapts and shifts the source mass, matching the ground-truth more closely.

Takeaway. The ill-posedness of source identification makes any estimate dependent on prior assumptions. Our method makes these explicit through the kernel; the PINN encodes them implicitly in the architecture and training procedure. In this case, a Matérn prior results in a better point estimate, and we additionally obtain confidence intervals and samples. We achieve all of this at a fraction of the cost of a PINN due to our proximity to classical numerical methods. We study the calibration of the posterior on this source identification problem in 9, and 11 extends the same framework to a genuinely nonlinear inverse problem (Burgers source identification with a softplus link).

4.2 Accuracy and Runtime on Forward Problems↩︎

We benchmark our method on a forward problem: the 2D viscous Burgers’ equation with bell-shaped initial conditions. We compare four methods: classical Crank-Nicolson FVM, a collocation-based sparse GP approach as in [12], and our GP-FVM in both full and minimal state configurations. The minimal configuration uses 8 state components per cell (values, derivatives, integrals), while the full configuration additionally tracks flux quantities. For scalability in the time dimension, all GP approaches use our moment-matching technique. For each method, we vary the grid resolution \(N \in \{8, 12, 16, 20, 24\}\) and average results over 5 random problem instances.

4 shows the results. Our GP-FVM (minimal) achieves accuracy within a factor of 2–3 of the classical method while providing uncertainty quantification. GP-Collocation, despite similar computational cost, performs an order of magnitude worse — highlighting the importance of conservative discretization. The full GP-FVM configuration achieves similar accuracy to the minimal variant but at substantially higher computational cost due to the larger state dimension. Full experimental specifications (problem setup, reference solver, kernel and lengthscale choices) are given in 8.2.

4.3 Convergence with Kernel Smoothness↩︎

In classical FVM, higher-order polynomial reconstructions yield faster convergence rates. In our method, the analogous mechanism is kernel smoothness: Matérn-\(\nu\) kernels with higher \(\nu\) yield smoother GP samples and faster convergence.

5 demonstrates this on a 1D Poisson problem (8.3). Panel (a) shows dense GP-FVM achieving rates close to the theoretical rate for Matérn-based kernel interpolation \(\mathcal{O}(h^{\nu+\frac{1}{2}})\) [29]. Panel (b) shows sparse GP-FVM with fixed \(\rho=5\) (16% fill at \(N=160\)). Matérn-3/2 maintains its convergence rate perfectly. For smoother kernels, errors plateau or increase at fine resolutions, where screening effects weaken and sparse approximation struggles to capture long-range correlations. Note however that this occurs in a regime of already low \(L^2\) error.

4.4 2D Shallow Water Equations↩︎

To demonstrate scalability to complex nonlinear systems, we simulate tsunami propagation using the shallow water equations (8.4). This is a system of three coupled conservation laws governing water height \(h\) and momentum \((hu, hv)\), with nonlinear flux terms that are quadratic in the conserved variables. The domain is 100 km \(\times\) 100 km with sloping bathymetry varying from 50 m depth at the shore to 500 m offshore, discretized on a \(31 \times 31\) grid. The spatially varying bathymetry introduces source terms coupling the momentum equations to the bottom topography.

6 shows the evolution of surface elevation \(\eta = h - b\) as a tsunami wave propagates toward shore. The method correctly captures wave propagation and shoaling—the steepening of waves as they enter shallow water—demonstrating physically plausible behavior without requiring problem-specific numerical treatment. This experiment shows that our framework extends naturally from scalar PDEs to complex systems of equations.

The joint state dimension per timestep is 28,566, comprising cell integrals, boundary values, and derivatives for all three conserved quantities. Our moment-matching scheme (3.4) requires approximately 80 seconds per timestep, of which roughly 60 seconds is spent computing marginal variances via selected inversion. While slower than classical solvers, this cost enables uncertainty propagation through the nonlinear dynamics within a Bayesian framework.

5 Related Work↩︎

GP-based PDE simulators. Probabilistic numerical methods cast numerical computation as Bayesian inference, providing a coherent probabilistic treatment of numerical uncertainty [9]. For PDEs, [10] show that physics-informed GP regression strictly generalizes classical methods of weighted residuals, including collocation, finite volume, and Galerkin methods. [13] develop a probabilistic FVM using integral functionals, but their method is restricted to linear PDEs and their low-rank approach struggles to accurately quantify uncertainties for large-scale problems.

Scalable inference. [12] apply sparse KL-minimizing Cholesky factorizations to solve nonlinear PDEs via collocation. [21] obtain scalability through the SPDE representation of Matérn priors in a Finite Element Method approach. By contrast, our work focuses on FVM, which is of central interest in computational fluid dynamics due to its conservative properties. By extending sparse KL-minimizing Cholesky to integral functionals, we frame FVM through a scalable Bayesian inference mechanism. Neither of these related works provides a practical answer to scalability along the time dimension.

6 Discussion↩︎

We frame the finite volume method as nonlinear conditioning under a structured GP prior. An extension of Vecchia approximations to integrals, as well as a novel moment-matching approach to time propagation, enables our method to circumvent the prohibitive cubic scaling of GPs. As a result, our method solves inverse problems quickly and under uncertainty and it scales to high-dimensional nonlinear simulations, while inheriting the accuracy and order-of-accuracy structure of classical FVM — with kernel smoothness playing the role of polynomial reconstruction order.

At the same time, our moment-matching approach incurs substantial fill-in in the temporal propagation step, which slows the underlying linear algebra and restricts practical scalability to 2D problems. The spatial sparse Cholesky scales well (\(\mathcal{O}(N_s^{3/2})\) in 2D, \(\mathcal{O}(N_s^{2})\) in 3D), so the bottleneck lies in how time is handled rather than in the spatial machinery. This points toward future work on filtering approaches that preserve sparsity across time — for which our scalable spatial framework is precisely the substrate needed.

A second limitation concerns the expressiveness of the prior. For inverse problems with positivity-enforcing link functions (e.g.the softplus prior \(s = \log(1 + e^g)\) on the source field), a stationary Matérn prior on the latent \(g\) has constant marginal variance, and the saturating link compresses negative values of \(g\) to near-zero \(s\). In regions the posterior identifies as “no source present” (\(g \ll 0\)), the saturating link collapses the marginal variance of \(s\), so the model underestimates its uncertainty about whether a weak source might be present. This is the locally overconfident source posterior visible where the softplus saturates in our Burgers experiment (11). Non-stationary kernels and sparsity-promoting priors (e.g.horseshoe-type) would relax this constraint; both remain compatible with the sparse Cholesky framework underlying our method.

Finally, our method relies on spatial Matérn priors and hyperrectangular FVM cells for closed-form computations of kernel integrals. Numerical integration would likely lift these restrictions, though we leave a thorough evaluation to future work.

6.0.0.1 Scope and outlook.

Our contribution is twofold: a scalable methodology for FVM-based probabilistic inference under structured GP priors, and empirical evidence — forward accuracy on benchmark PDEs, \(\mathcal{O}(N_s^{3/2})\) runtime scaling, and posterior coverage ranging from near-nominal to conservative in a calibration study — that the resulting posteriors are quantitatively meaningful. The Laplace approximation underlying our nonlinear inversion is of course only Gaussian-around-the-MAP; where richer non-Gaussian uncertainty is required, our scalable posterior provides a natural starting point and preconditioner for downstream samplers (e.g.as the initial state and mass matrix for Hamiltonian Monte Carlo, or as a proposal distribution), and therefore complements rather than competes with Markov chain Monte Carlo-based approaches. Beyond the engineering value of accelerated inference, the broader case for probabilistic numerical methods is that they place numerical and inferential uncertainty on a single footing: the same posterior that quantifies measurement noise and prior assumptions also quantifies discretization uncertainty, enabling coherent fusion of sparse data with mechanistic knowledge. Making this combination computationally practical for nonlinear conservation laws — the dominant PDE class in the sciences — is what our method delivers.

Acknowledgements↩︎

The authors gratefully acknowledge co-funding by the European Union (ERC, ANUBIS, 101123955). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. PH is supported by the DFG through Project HE 7114/6-1 in SPP2298/2. PH is a member of the Machine Learning Cluster of Excellence, funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC number 2064/1 – Project number 390727645. The authors also gratefully acknowledge the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ:01IS18039A); and funds from the Ministry of Science, Research and Arts of the State of Baden-Württemberg. The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting TW.

Impact Statement↩︎

This paper develops scalable Bayesian inference for nonlinear conservation laws, with the explicit goal of providing uncertainty quantification alongside numerical PDE solutions. The intended downstream uses include scientific decision-making under uncertainty — environmental monitoring, geophysical inversion, contaminant transport, and similar settings. A specific risk we wish to flag: the calibration of our posterior depends on prior choices (e.g.the lengthscale, smoothness, and stationarity of the spatial kernel) that encode modeling assumptions which may not hold for a given real-world problem. Our calibration study in 9 demonstrates broadly reasonable calibration — near-nominal to conservative — on a synthetic benchmark, but we caution against treating posterior uncertainty as automatically trustworthy in high-stakes settings (e.g.safety-critical or regulatory applications) without independent validation. We view the method as one source of evidence to be combined with domain expertise, sensitivity analyses, and corroborating data.

7 Theoretical Details↩︎

7.1 Integrated Wiener Process and Temporal Precision↩︎

The Integrated Wiener Process (IWP) of order \(q\) is a Gauss-Markov process whose \(q\)-th derivative is Brownian motion. This Markov property induces a block tridiagonal structure in the joint precision matrix over discrete time points.

7.1.0.1 State-space formulation.

The IWP(\(q\)) can be written as a linear stochastic differential equation (SDE). Define the state vector \(\boldsymbol{x}(t) := [u(t), u'(t), \ldots, u^{(q-1)}(t)]^\top \in \mathbb{R}^q\), collecting the function value and its first \(q-1\) derivatives. The SDE is \[d\boldsymbol{x}(t) = \boldsymbol{F} \boldsymbol{x}(t) \, dt + \boldsymbol{B} \, dW(t), \label{eq:iwp-sde}\tag{14}\] where \(W(t)\) is a Wiener process, \(\boldsymbol{F} \in \mathbb{R}^{q \times q}\) is the shift matrix with ones on the superdiagonal and zeros elsewhere, and \(\boldsymbol{B} = [0, \ldots, 0, \sigma]^\top \in \mathbb{R}^q\) drives only the highest derivative.

7.1.0.2 Discretization.

Discretizing 14 at time points \(t_1, \ldots, t_{N_t}\) with step size \(h = t_{i+1} - t_i\) yields the transition \[\boldsymbol{x}_{i+1} = \boldsymbol{A} \boldsymbol{x}_i + \boldsymbol{\eta}_i, \quad \boldsymbol{\eta}_i \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_\eta), \label{eq:iwp-transition}\tag{15}\] where \(\boldsymbol{A} = \exp(\boldsymbol{F} h)\) is the state transition matrix. We have the closed forms [19]: \[\begin{align} A_{ij} & = \frac{h^{j-i}}{(j-i)!} \quad \text{for } j \geq i, \quad A_{ij} = 0 \quad \text{otherwise}, \\ (\Sigma_\eta)_{ij} & = \frac{\sigma^2 h^{2q+1-i-j}}{(q-i)!(q-j)!(2q+1-i-j)}. \end{align}\]

7.1.0.3 Block tridiagonal precision.

The Markov property in 15 implies that the joint distribution \(p(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_{N_t})\) factorizes as \[p(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_{N_t}) = p(\boldsymbol{x}_1) \prod_{i=1}^{N_t - 1} p(\boldsymbol{x}_{i+1} \mid \boldsymbol{x}_i).\] Each transition density is Gaussian: \[p(\boldsymbol{x}_{i+1} \mid \boldsymbol{x}_i) \propto \exp\!\left( -\tfrac{1}{2} (\boldsymbol{x}_{i+1} - \boldsymbol{A} \boldsymbol{x}_i)^\top \boldsymbol{\Pi} (\boldsymbol{x}_{i+1} - \boldsymbol{A} \boldsymbol{x}_i) \right),\] where \(\boldsymbol{\Pi} := \boldsymbol{\Sigma}_\eta^{-1}\) is the precision of the process noise. Expanding and collecting quadratic terms across all time steps, the joint precision matrix takes the block tridiagonal form \[\boldsymbol{Q}_{\tilde{T}} = \begin{bmatrix} \boldsymbol{\Pi}_0 + \boldsymbol{A}^\top \boldsymbol{\Pi} \boldsymbol{A} & -\boldsymbol{A}^\top \boldsymbol{\Pi} & & \\ -\boldsymbol{\Pi} \boldsymbol{A} & \boldsymbol{\Pi} + \boldsymbol{A}^\top \boldsymbol{\Pi} \boldsymbol{A} & -\boldsymbol{A}^\top \boldsymbol{\Pi} & \\ & \ddots & \ddots & \ddots \\ & & -\boldsymbol{\Pi} \boldsymbol{A} & \boldsymbol{\Pi} \end{bmatrix}, \label{eq:iwp-precision}\tag{16}\] where \(\boldsymbol{\Pi}_0\) is the precision of the initial state \(\boldsymbol{x}_1\). Each block is \(q \times q\), so for \(N_t\) time steps the full matrix is \(q N_t \times q N_t\) with only \(\mathcal{O}(q^2 N_t)\) nonzero entries. This sparsity is what enables efficient temporal inference.

7.2 Crank-Nicolson Time Discretization↩︎

For time-dependent problems, we enforce the FVM constraints using Crank-Nicolson discretization, which is second-order accurate in time and unconditionally stable for linear advection-diffusion problems.

7.2.0.1 Two-step joint distribution.

At each time step, we form the joint distribution over consecutive states \((\boldsymbol{s}_{t-1}, \boldsymbol{s}_t)\) as described in 3.4. This joint distribution incorporates both the temporal prior (IWP transition) and the current spatial precision.

7.2.0.2 Crank-Nicolson constraint.

The FVM constraint 7 is enforced by averaging all quantities at consecutive time levels: \[\label{eq:cn95averaged} \frac{1}{2}\left( \dot{\bar{u}}_i^{(t-1)} + \dot{\bar{u}}_i^{(t)} \right) + \frac{1}{2}\left( \mathcal{F}_i^{(t-1)} + \mathcal{F}_i^{(t)} \right) = \frac{1}{2}\left( s_i^{(t-1)} + s_i^{(t)} \right),\tag{17}\] where \(\dot{\bar{u}}_i\) denotes the time derivative of the cell integral (which is explicitly available in the IWP state), \(\mathcal{F}_i\) the net flux through cell boundaries, and \(s_i\) the source integral.

7.2.0.3 Conditioning.

We enforce this constraint on the joint distribution \((\boldsymbol{s}_{t-1}, \boldsymbol{s}_t)\) via the nonlinear conditioning procedure in 2.3. For linear PDEs, this reduces to a single linear constraint. For nonlinear PDEs, inspired by the Extended Kalman Filter (EKF), we take a single Gauss-Newton step rather than iterating to convergence. This is justified by the small time step assumption: with sufficiently small \(\Delta t\), the nonlinearity within one step is mild, making a single linearization accurate. Any residual linearization error is corrected at subsequent time steps rather than accumulating, and the computational savings are substantial; we validate this empirically in 12. After conditioning, we extract the marginal distribution of \(\boldsymbol{s}_t\) and propagate to the next step.

7.2.0.4 Laplace approximation error.

The Laplace approximation used in the per-step conditioning is exact for linear constraints and approximate for nonlinear ones. Combining the averaged conservation law 17 with the cell-wise PDE \(\dot{\bar{u}}_i = s_i - \mathcal{F}_i\) and rearranging implicitly for \(\boldsymbol{s}_t\), the constraint function has the structure \[\label{eq:cn95constraint95struct} g(\boldsymbol{s}_t) = \boldsymbol{s}_t - \tfrac{\Delta t}{2}\, f(\boldsymbol{s}_t) - (\text{known from } \boldsymbol{s}_{t-1}),\tag{18}\] i.e.identity plus an \(O(\Delta t)\) nonlinear term. A short computation shows that all derivatives of \(g\) beyond the first inherit the \(\Delta t/2\) prefactor: \(D^k g = O(\Delta t)\) for \(k \geq 2\), whenever the flux \(f\) is sufficiently smooth. Heuristically, a perturbative (Edgeworth-style) expansion of the negative log-posterior around the MAP then gives the following per-step rate.

7.3 Node Ordering for Mixed Functionals↩︎

When the state vector contains multiple types of linear functionals (point evaluations, derivatives, integrals), the sparse Cholesky approximation quality depends on how these blocks are ordered. We propose ordering functionals from finest to coarsest scale, placing integrals in the rightmost columns of the Cholesky factor.

7.3.0.1 Functional classification.

We classify linear functionals into a hierarchy based on the spatial scale they capture:

• Integrals (coarsest): Cell averages \(\iint f \, dx\,dy\) capture global, smoothed behavior

• Face integrals: Mixed tensor products such as \(\int f(x, \cdot) \, dy\) at fixed \(x\)

• Evaluations: Point values \(f(x_i)\) capture local behavior

• Derivatives (finest): Local slopes \(\partial f / \partial x\) capture fine-scale variation

Our Julia implementation automatically classifies linear functionals via multiple dispatch on their types.

7.3.0.2 Ordering algorithm.

We order functional blocks from finest to coarsest (derivatives \(\to\) evaluations \(\to\) integrals), placing the coarsest observations in the rightmost columns of the Cholesky factor. Within each block, we apply the reverse maximin ordering of [18].

7.3.0.3 Why integrals should be coarsest.

In the sparse Cholesky factorization, each row \(i\) conditions on a subset of later rows \(j > i\). The screening effect [18] implies that coarse-scale measurements effectively “screen out” correlations, reducing the conditioning set sizes needed for accurate approximation.

Integrals are smoothing operators that average over spatial regions, making them natural coarse-scale measurements. When integrals occupy the rightmost (coarsest) positions, all finer-scale functionals to their left benefit from this strong screening effect, requiring smaller conditioning sets to achieve the same approximation quality.

Additionally, there is a geometric asymmetry: point evaluations are more densely packed than integral cells covering the same domain. When integrals (in the middle) condition on evaluations (coarsest), each integral’s neighborhood contains many evaluation points, leading to large conditioning sets. Conversely, when evaluations condition on integrals, each evaluation finds fewer integral neighbors. This asymmetry means that for the same sparsity parameter \(\rho\), placing integrals coarsest yields lower fill-in.

7.3.0.4 Empirical validation.

We compare two orderings on a 2D test problem: a \(35 \times 35\) grid with point evaluations, derivative sums (\(\partial_x + \partial_y\)), and cell integrals under a Matérn-3/2 kernel. For each ordering, we sweep the sparsity parameter \(\rho \in \{1, 1.5, 2, \ldots, 10\}\) and measure the KL divergence from the exact posterior.

7 shows that the “integrals first” ordering Pareto-dominates “evaluations first”: for comparable fill-in levels, it achieves lower approximation error. This validates placing integrals as the coarsest block: both the screening effect and the geometric asymmetry favor this ordering.

7.3.0.5 Visualizing the exact Cholesky factor.

8 shows the exact precision Cholesky factor \(\boldsymbol{L}\) (where \(\boldsymbol{\Theta} = \boldsymbol{L}\boldsymbol{L}^\top\)) for a smaller \(12 \times 12\) grid. The color scale shows \(\log_{10}|L_{ij}|\), revealing that many remote entries are already close to zero in the exact factorization. This decay justifies the sparse approximation: entries below a threshold can be dropped with minimal loss in accuracy. The dashed lines mark boundaries between functional blocks (derivatives, evaluations and integrals).

8 Experimental Details↩︎

All GP-based experiments were run on an Apple M1 Max with 32 GB RAM. The PINN baseline was run on an NVIDIA GeForce RTX 2080 Ti.

Our Julia implementation builds on GaussianMarkovRandomFields.jl [30] for information form inference. This in turn builds on sparse linear algebra routines from SuiteSparse, specifically CHOLMOD [26]. The sparse Jacobians required for Gauss-Newton steps are computed via sparse automatic differentiation, as in [31].

8.1 Source Identification in Advection-Diffusion↩︎

This section provides implementation details for the source identification experiment in 4.1.

8.1.0.1 Problem specification.

The domain is the unit square \([0,1]^2\) discretized on a uniform \(36 \times 36\) node grid (35 \(\times\) 35 cells). The ground-truth source field consists of two Gaussian sources: \[\begin{align} s_1(x,y) & = 4.0 \cdot \exp\!\left( -\frac{(x - 0.25)^2 + (y - 0.3)^2}{2 \cdot 0.10^2} \right), \\ s_2(x,y) & = 3.0 \cdot \exp\!\left( -\frac{(x - 0.35)^2 + (y - 0.7)^2}{2 \cdot 0.12^2} \right). \end{align}\] The physics parameters are velocity \(\boldsymbol{v} = (1, 0)^\top\) and diffusion coefficient \(D = 0.05\), yielding a Péclet number of \(\text{Pe} = v_x L / D = 20\) (advection-dominated). Boundary conditions are: Dirichlet \(c = 0\) at the left boundary (clean inflow), homogeneous Neumann at the right boundary (free outflow), and homogeneous Neumann at the top and bottom boundaries (impermeable walls).

Ground-truth concentration data is generated via a forward FVM solve with upwind advection. Observations are taken at 12 monitoring well locations spread across the domain (1), with additive Gaussian noise \(\sigma = 0.01\).

8.1.0.2 GP-FVM implementation.

We use Matérn-5/2 kernels (\(\nu = 2\)) with 2D tensor product structure. The lengthscale is set to \(\ell = 5 \Delta x \approx 0.14\), chosen to resolve approximately 5 grid cells. The sparsity parameter is \(\rho = 2.0\) for the concentration field and \(\rho = 4.0\) for the source field. The joint state vector comprises:

• Concentration: point evaluations \(c(x_i, y_j)\), vertical face integrals \(\int c \, dy\), horizontal face integrals \(\int c \, dx\), derivative face integrals \(\int \partial_x c \, dy\) and \(\int \partial_y c \, dx\)

• Source: point evaluations \(s(x_i, y_j)\) and cell integrals \(\iint_{\text{cell}} s \, dx\,dy\)

All kernel integrals are computed in closed form using the Matérn polynomial structure (no quadrature). FVM and boundary condition constraints are enforced with precision \(1/\sigma_{\text{constraint}}^2\) where \(\sigma_{\text{constraint}} = 10^{-5}\).

8.1.0.3 PINN baseline.

We implement the PINN baseline using the JAX-based jaxpi framework [28]. The architecture uses two separate MLPs:

• Concentration network: 4 layers \(\times\) 128 units, tanh activation, with Fourier feature embeddings (scale \(= 1.0\), dimension \(= 128\))

• Source network: 2 layers \(\times\) 64 units, tanh activation, no Fourier features

We train with the Adam optimizer (\(\beta_1 = 0.9\), \(\beta_2 = 0.999\), initial learning rate \(10^{-3}\)) using exponential decay (rate \(0.9\) every 2000 steps) for 30,000 iterations. We use fixed loss weights: PDE residual \(= 1.0\), boundary conditions \(= 1.0\), observation data \(= 1000.0\). We sample 4,096 interior collocation points and 100 boundary points per edge (400 total) at each iteration. All experiments use random seed 42. We also tested a larger source network (4 layers × 128 units with Fourier features), but found no improvement in reconstruction accuracy.

8.2 2D Burgers Benchmark↩︎

This section provides details for the accuracy and runtime comparison in 4.2.

8.2.0.1 Problem specification.

We solve the 2D viscous Burgers equation in conservative form on the unit square \([0,1]^2\) with periodic boundary conditions: \[\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}\left(\frac{1}{2}u^2\right) + \frac{\partial}{\partial y}\left(\frac{1}{2}u^2\right) = \nu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right),\] with viscosity \(\nu = 0.02\) and end time \(T = 0.5\). Initial conditions are Gaussian bumps with randomly sampled parameters: width \(\sigma \in [0.2, 0.3]\), amplitude \(A \in [0.3, 0.6]\), and center \((x_0, y_0) \in [0.2, 0.5]^2\).

8.2.0.2 Reference solution.

We compute ground-truth solutions using high-resolution classical FVM on a \(128 \times 128\) grid with RK4 time integration and a stable CFL condition (\(\Delta t = 0.25 \Delta x^2\)). Errors are computed as \(L^2\) norm against this reference, interpolated to the test grid.

8.2.0.3 Methods compared.

• Classical Crank-Nicolson FVM: Standard implicit scheme with \(N^2\) unknowns.

• GP-Collocation: Sparse GP approach following [12], using point evaluations and derivatives. State dimension: \(8N^2\).

• GP-FVM (minimal): Our method with point evaluations, derivatives, and cell integrals. State dimension: \(8N^2\).

• GP-FVM (full): Adds flux-related quantities. State dimension: \(26N^2\).

All methods use Crank-Nicolson time discretization (7.2). The GP methods additionally use Matérn-5/2 kernels with lengthscale \(\ell = 0.3\) and our moment-matching scheme (3.4) for temporal propagation.

8.2.0.4 Experimental parameters.

Grid sizes: \(N \in \{8, 12, 16, 20, 24\}\). Number of time steps scales with grid size: \(n_{\text{steps}} \approx 1.5 N\), so \(\Delta t = T / n_{\text{steps}}\). We evaluate on five random problem instances.

8.3 Poisson Convergence Study↩︎

This section provides details for the convergence experiment in 4.3.

8.3.0.1 Problem specification.

We solve the 1D Poisson equation on the unit interval: \[-u''(x) = f(x), \quad x \in (0, 1), \quad u(0) = u(1) = 0.\] The source is \(f(x) = \sin(2\pi x)\), with exact solution \(u(x) = \sin(2\pi x) / (4\pi^2)\).

8.3.0.2 FVM formulation.

Integrating over cell \([x_i, x_{i+1}]\) and applying the fundamental theorem of calculus: \[u'(x_i) - u'(x_{i+1}) = \int_{x_i}^{x_{i+1}} f(x) \, dx.\] This links derivatives of \(u\) at cell boundaries to integrals of the source \(f\).

8.3.0.3 GP-FVM implementation.

We place independent GP priors on \(u\) and \(f\), both with Matérn-\(\nu\) kernels and lengthscale \(\ell = 0.15\). The joint state for \(u\) includes point evaluations \(u(x_i)\), derivatives \(u'(x_i)\), and cell integrals \(\int u \, dx\). The joint state for \(f\) includes point evaluations \(f(x_i)\) and cell integrals \(\int f \, dx\). We prescribe the source function \(f\) at grid points, enforce the FVM constraints, and apply Dirichlet boundary conditions. All constraints use precision \(\sigma_{\text{constraint}}^{-2}\) with \(\sigma_{\text{constraint}} = 10^{-6}\).

8.3.0.4 Experimental parameters.

Grid sizes: \(N \in \{10, 20, 40, 80, 160\}\). Kernel smoothness: \(\nu \in \{3/2, 5/2, 7/2\}\). For dense GP-FVM, we use \(\rho = 60\) (essentially exact). For sparse GP-FVM, we use fixed \(\rho = 5\), yielding approximately 16% fill-in at \(N = 160\).

8.4 Shallow Water Equations↩︎

The nonlinear shallow water equations describe the evolution of water height \(h\) and momentum \((hu, hv)\) over a domain with bathymetry \(b(x,y)\): \[\begin{align} \frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} + \frac{\partial (hv)}{\partial y} & = 0, \tag{20} \\ \frac{\partial (hu)}{\partial t} + \frac{\partial}{\partial x}\left(\frac{(hu)^2}{h} + \frac{1}{2}gh^2\right) + \frac{\partial}{\partial y}\left(\frac{(hu)(hv)}{h}\right) & = -gh\frac{\partial b}{\partial x}, \tag{21} \\ \frac{\partial (hv)}{\partial t} + \frac{\partial}{\partial x}\left(\frac{(hu)(hv)}{h}\right) + \frac{\partial}{\partial y}\left(\frac{(hv)^2}{h} + \frac{1}{2}gh^2\right) & = -gh\frac{\partial b}{\partial y}, \tag{22} \end{align}\] where \(g\) is gravitational acceleration. The right-hand side terms represent pressure forces from sloping bathymetry.

8.4.0.1 Experimental setup.

We simulate tsunami propagation over a 100 km \(\times\) 100 km domain with linearly sloping bathymetry: depth varies from 50 m at the shore (\(x=0\)) to 500 m offshore (\(x=L_x\)). The initial condition is a Gaussian perturbation of amplitude 10 m centered at \((70\,\text{km}, 50\,\text{km})\), representing a tsunami source, with zero initial velocity. We do not include Coriolis forces or bottom friction.

9 Calibration Study↩︎

Posterior standard deviation contracting near measurement locations is a generic property of GP regression; on its own, it does not establish that the posterior is calibrated, i.e.that credible intervals achieve their nominal coverage. We assess calibration empirically via a multi-instance study on the source identification problem of 4.1.

9.0.0.1 Data generating process.

We generate 200 random source identification instances. Each instance has two Gaussian source bumps at random locations in \([0.2, 0.8]^2\) with random amplitudes in \([-3, 3]\) and widths in \([0.05, 0.15]\). For each instance, we solve the forward advection-diffusion PDE to obtain the true concentration, sample 10 noisy observations (\(\sigma = 0.1\)) at fixed locations, and run the GP-FVM solver on a \(31 \times 31\) grid.

9.0.0.2 Output scale calibration.

The GP prior’s output scale \(\sigma^2\) controls the overall magnitude of the posterior uncertainty. We calibrate it via a two-pass procedure: (i) run 10 calibration instances with an initial \(\sigma^2\), compute the mean squared marginal \(z\)-score \(\bar{z}^2 = \mathbb{E}[(f^\dagger - \mu)^2 / \sigma_\text{post}^2]\); (ii) set \(\sigma^2_\text{cal} = \sigma^2_\text{init} \cdot \bar{z}^2\) so that \(z\)-scores have unit variance on average. We calibrate the source and concentration fields independently, obtaining \(\sigma^2_s = 0.154\) and \(\sigma^2_c = 0.347\).

9.0.0.3 Results.

9 shows the calibration plot (nominal vs.empirical coverage) across all 200 instances. The concentration field (blue) tracks the diagonal closely across most of the range, with mild deviations at the extremes: empirical coverage is somewhat above nominal at low levels (e.g.\(59\%\) at \(50\%\) nominal) and somewhat below at high levels (\(88\%\) at \(95\%\) nominal). The source field (red) is broadly conservative — empirical coverage substantially exceeds nominal across the low and middle range (e.g.\(80\%\) at \(50\%\) nominal, \(93\%\) at \(90\%\) nominal) — and crosses the diagonal near the \(95\%\) level.

10 shows the spatial coverage map at the 95% level, averaged across all instances. The source field achieves near-nominal coverage at the domain boundaries (where the Dirichlet BC constrains the solution) and is slightly overconfident in the interior where the source bumps are located. This spatial pattern reflects the prior–truth mismatch: the stationary Matérn prior spreads uncertainty uniformly, while the actual error is concentrated at source locations.

9.0.0.4 Discussion.

The remaining miscalibration has two distinct sources:

1. Prior–truth mismatch in the source field. The stationary Matérn prior spreads uncertainty uniformly across the domain, while the actual reconstruction errors are concentrated at source locations. This produces the consistent over-coverage of the source curve at low and middle nominal levels: posterior credible intervals are wider than needed across most of the domain. Sparsity-promoting hierarchical GMRF priors (via local variance parameters with an appropriate hyperprior) or log-Gaussian priors for positivity would improve calibration while preserving the sparse Cholesky structure.

2. Tail behavior of the concentration field. The output-scale calibration matches the mean squared \(z\)-score to unity, ensuring the right average posterior spread but not the full shape of the empirical error distribution. The mild over-coverage at low nominal levels combined with mild under-coverage at high nominal levels is consistent with empirical errors that are heavier-tailed than the assumed Gaussian posterior. Heavier-tailed likelihood models or non-Gaussian posterior approximations could address this within the same sparse Cholesky framework.

Both limitations have concrete paths to improvement within the framework and do not require architectural changes.

10 Scalability Study↩︎

To complement the asymptotic complexity analysis of 3.5, we report wall-clock timings and reconstruction errors for the source identification problem of 4.1 across a range of grid resolutions.

10.0.0.1 Setup.

We run the 2D steady advection-diffusion source identification problem from 4.1 at grid sizes \(N = 11\) to \(N = 101\) (\(N \times N\) grid points, corresponding to \(N_s = 782\) to \(N_s = 70{,}802\) total degrees of freedom). We use physics-motivated lengthscales (\(\ell_c = 0.2 \approx \sqrt{D/v_x}\), \(\ell_s = 0.12\)) that are fixed across all grid sizes, and output scales \(\sigma^2_s = 0.154\), \(\sigma^2_c = 0.347\) taken from the calibration study (9). Timing is the best of 3 runs per grid size.

10.0.0.2 Timing.

11 shows wall-clock time versus total DOF (\(N_s\)) on a log-log scale. The total time tracks the theoretical \(\mathcal{O}(N_s^{3/2})\) reference line closely (empirical exponent \(\approx 1.35\)), confirming the complexity analysis in 3.5. The dominant cost at large \(N_s\) is posterior statistics extraction (selected inversion on the sparse Cholesky factor). The method solves the inverse problem with \(\sim 70{,}000\) DOF in under 4 minutes on a single CPU core.

10.0.0.3 Error.

2 shows that solution error is stable across grid sizes: source RMSE and concentration RMSE plateau beyond \(N = 31\), consistent with the observation-limited regime (10 noisy measurements constrain the inverse problem, not the grid resolution). The source location error remains below one grid cell for all \(N \geq 21\).

The error plateau confirms that the method does not degrade at larger grid sizes — the accuracy is limited by the observation information, not the discretization. The stable fill-in percentages (decreasing from \(37.5\%\) at \(N = 11\) to \(0.6\%\) at \(N = 101\) for the concentration field) demonstrate that the sparse Cholesky approximation becomes increasingly efficient at larger scales.

11 Nonlinear Inverse Problem: Burgers Source Identification↩︎

We demonstrate our framework on a genuinely nonlinear inverse problem: inferring an unknown spatial source field \(s(x,y) \geq 0\) from noisy observations of the solution \(u(x,y,t)\) to the 2D viscous Burgers equation, \[\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}\!\left(\frac{u^2}{2}\right) + \frac{\partial}{\partial y}\!\left(\frac{u^2}{2}\right) = \nu \Delta u + s(x,y), \label{eq:burgers95source}\tag{23}\] on \([0,1]^2 \times [0, 0.5]\) with \(\nu = 0.01\), homogeneous Dirichlet boundary conditions, and Gaussian bump initial condition. The source \(s\) is constant in time and couples all time steps.

11.0.0.1 Setup.

The source consists of two Gaussian bumps with amplitudes \(\approx 1\) and \(\approx 3.5\) and widths \(\approx 0.12\)–\(0.14\). We observe \(u\) at 50 random spatial locations at 5 out of 10 time steps (250 observations total, noise \(\sigma = 0.05\)). The source is never observed directly — it must be inferred entirely through the nonlinear PDE.

11.0.0.2 Prior and solver.

We place a Matérn-5/2 GP prior on \(u\) with IWP(1) temporal dynamics and a separate Matérn-5/2 GP prior on \(g\), where \(s = \mathrm{softplus}(g) = \log(1 + e^g) \geq 0\). The softplus link enforces positivity while keeping bounded Jacobian entries (unlike \(\exp\)), which improves the quality of the Gauss–Newton approximation. The joint spacetime GMRF has \(\sim 22{,}000\) degrees of freedom (\(16 \times 16\) spatial grid, 11 time steps). The Gauss–Newton iteration converges in 14 iterations (3), with Newton decrement dropping over 10 orders of magnitude.

11.0.0.3 Results.

12 shows the solution at five time steps (IC, two unobserved, two observed). The GP-FVM posterior mean tracks the nonlinear advection–diffusion dynamics closely. The posterior standard deviation (bottom row) is small near observation locations (white dots) and larger between observed time steps, reflecting the expected information pattern.

13 shows the inferred source field. The posterior mean identifies both source locations and captures the dominant bump’s amplitude. The posterior standard deviation (in \(g\)-space, interpretable as multiplicative uncertainty) is largest where the source is active, reflecting the difficulty of the inverse problem.

11.0.0.4 Discussion.

This experiment demonstrates two properties of the framework: (i) it handles a genuinely nonlinear inverse problem (Burgers flux \(F = u^2/2\)), and (ii) the Gauss–Newton iteration converges reliably to a good solution in \(\sim 10\)–\(14\) iterations. The softplus source prior demonstrates that nonlinear link functions (and the associated positivity constraints) integrate naturally into the GP-FVM framework — the same machinery handles both the PDE nonlinearity and the prior nonlinearity.

The remaining locally underestimated uncertainty in the source posterior (visible where the softplus is saturated) is discussed in detail in 6.

12 Gauss–Newton Ablation↩︎

The sequential (EKF-style) solver of 3.4 performs a single Gauss–Newton step per time step. This is a strong assumption, which we validate empirically in two complementary settings: a forward-simulation ablation showing that additional iterations do not improve accuracy, and a convergence study confirming that the joint solver used for inverse problems requires multiple iterations but reaches them reliably.

12.0.0.1 Ablation: sequential solver.

We run the 1D viscous Burgers forward problem with the EKF sequential solver, varying the number of Gauss–Newton iterations per time step from 1 (standard EKF) to 10. 14 shows the result across 5 initial conditions (sine and Gaussian families) at \(N = 100\) grid points: additional iterations provide no measurable improvement in solution accuracy (\(< 0.1\%\) change) while adding \(\sim 10\%\) computational overhead.

This validates the single-step assumption: in the small-\(\Delta t\) regime, the per-step nonlinearity is a mild perturbation, and one linearization suffices. The Gauss–Newton iteration converges in a single step because the prior from the previous time step already provides an excellent linearization point.

12.0.0.2 Convergence: joint solver.

For the nonlinear inverse problem (11), where we solve the full spacetime system jointly, multiple Gauss–Newton iterations are needed — the problem cannot be decomposed into small-\(\Delta t\) steps. The iteration converges reliably in \(\sim 10\)–\(14\) steps (3), with the Newton decrement dropping over 10 orders of magnitude.

13 Moment-Matching Error Analysis↩︎

We provide a theoretical analysis of the moment-matching approximation introduced in 3.4. We frame it within the assumed density filtering framework, derive a per-step projection error decomposition that explicitly accounts for the diagonal mismatch arising from the closed-form variance rule, and bound the accumulated KL divergence from the exact filter via a contraction argument. The analysis shows that the moment-matching error decomposes into an off-diagonal contribution that vanishes as \(\Delta t \to 0\) and a diagonal mismatch floor controlled by the spatial prior’s conditioning.

13.1 Setup and Notation↩︎

Standing notation.

• \(\Sigma_t^+\): marginal filtering covariance of \(x_t\) after conditioning at step \(t\).

• \(\Sigma_t^-\): marginal predictive covariance before conditioning.

• \(A = A(\Delta t)\): IWP transition matrix; \(\Sigma_\eta = \Sigma_\eta(\Delta t)\): process noise covariance.

• \(\widetilde{\Sigma}_t = (Q_s + \widetilde{D}_t)^{-1}\): moment-matched covariance, with \(\widetilde{D}_t = \mathop{\mathrm{diag}}(\sigma_t^{-2})\).

• \(\sigma_{t,i}^2 := [\Sigma_t^+]_{ii}\): the true marginal variances after conditioning.

• \(\Delta_t := G_t S_t^{-1} G_t^\top \succeq 0\): covariance removed by FVM conditioning.

• \(\lambda_{\min}(\cdot)\), \(\lambda_{\max}(\cdot)\): smallest and largest eigenvalues.

• \(\kappa(M) := \lambda_{\max}(M)/\lambda_{\min}(M)\): condition number.

• \(N\): state dimension; \(N_s\): number of spatial cells; \(m\): observation dimension per step.

Key spectral parameters. \[\label{eq:spectral} \lambda_{\min}(Q_s) =: \alpha > 0, \qquad \lambda_{\max}(Q_s) =: \beta, \qquad \kappa_s := \beta/\alpha.\tag{25}\] The filtering precision satisfies \(Q_s + \widetilde{D}_t \succeq Q_s \succeq \alpha I\), giving the covariance upper bound \[\label{eq:cov95upper} \widetilde{\Sigma}_t \preceq Q_s^{-1} \preceq \alpha^{-1} I.\tag{26}\] We assume throughout that \(\lVert \widetilde{D}_t \rVert_{\mathrm{op}} \leq \delta_{\max}\) uniformly in \(t\) for some finite \(\delta_{\max} > 0\); this is a mild filter-stability condition, satisfied whenever the moment-matched marginals are bounded away from zero. Under this assumption the covariance lower bound is \[\label{eq:cov95lower} \widetilde{\Sigma}_t \succeq (\beta + \delta_{\max})^{-1} I.\tag{27}\]

13.2 The Diagonal Mismatch↩︎

The algorithm stores \(\sigma_{t,i}^2 := [\Sigma_t^+]_{ii}\) and sets \(\widetilde{D}_t = \mathop{\mathrm{diag}}(1/\sigma_{t,i}^2)\). If \(Q_s = 0\), then \(\widetilde{\Sigma}_t = \widetilde{D}_t^{-1}\) and \([\widetilde{\Sigma}_t]_{ii} = \sigma_{t,i}^2\) exactly. But for \(Q_s \neq 0\), the Schur complement formula gives \[\label{eq:schur} [(Q_s + D)^{-1}]_{ii} = \frac{1}{d_i + [Q_s]_{ii} - \mathbf{q}_i^\top(Q_{s,-i} + D_{-i})^{-1}\mathbf{q}_i},\tag{28}\] where \(\mathbf{q}_i\) is the \(i\)-th row of \(Q_s\) with the diagonal removed, and \(Q_{s,-i}\), \(D_{-i}\) are the corresponding submatrices. The correction term \(c_i := \mathbf{q}_i^\top(Q_{s,-i} + D_{-i})^{-1}\mathbf{q}_i \geq 0\) makes the actual marginal variance differ from \(1/(d_i + [Q_s]_{ii})\), shifting it away from the target \(\sigma_{t,i}^2 = 1/d_i\).

13.3 Per-Step Update↩︎

The prediction and update follow the standard Kalman recursion, with the approximate filter propagating \(\widetilde{\Sigma}_{t-1}\) (the moment-matched covariance) rather than \(\Sigma_t^+\): \[\begin{align} \Sigma_t^- &= A\,\widetilde{\Sigma}_{t-1}\,A^\top + \Sigma_\eta, \tag{31}\\ \Sigma_t^+ &= \Sigma_t^- - \Delta_t, \qquad \Delta_t := G_t S_t^{-1} G_t^\top \succeq 0. \tag{32} \end{align}\]

13.4 Per-Step Projection Error↩︎

The covariance error \(E_t := \Sigma_t^+ - \widetilde{\Sigma}_t\) has two contributions: (a) off-diagonal discrepancy from projecting onto the sparse precision family; and (b) diagonal discrepancy from the approximate moment-matching rule.

13.5 Total Projection Error↩︎

13.6 Accumulated Error via Contraction↩︎

References↩︎

---

1. Code: https://github.com/timweiland/GPFiniteVolume.jl↩︎