Structure-Preserving Neural Surrogates with Tractable Uncertainty Quantification

# Structure-Preserving Neural Surrogates with Tractable Uncertainty Quantification ††thanks: \fundingThis material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under award numbers DE-SC0024563 and DE-SC0023163.
Source: [https://arxiv.org/html/2606.11650](https://arxiv.org/html/2606.11650)
\\newsiamremark

remarkRemark\\newsiamremarkhypothesisHypothesis\\newsiamthmclaimClaim\\newsiamremarkfactFact\\headersStructure\-Preserving Neural Surrogates with Tractable UQZhang et al\.\\externaldocument\[\]\[nocite\]ex\_supplement

Handi ZhangApplied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, PA, USA \(\)\.Adrienne M\. ProppInstitute for Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA \(\)\.Brooks KinchHouman OwhadiDepartment of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA, USA \(\)\.Nathaniel TraskMechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA, USA \(\)\.

###### Abstract

Recent advances in scientific machine learning provide a means of near\-real\-time solution to partial differential equations \(PDEs\), but lack the theoretical underpinnings of conventional simulators that support contemporary verification and validation\. In this work, we construct data\-driven reduced\-order models that serve as structure\-preserving, real\-time surrogates\. Remarkably, the exterior calculus that imposes physical conservation structure also exposes topological structure that we use to build a Gaussian process \(GP\) representation of uncertainty in state\-flux relationships, ultimately yielding a Dirichlet\-to\-Neumann map for quantities of interest with closed\-form expressions for posterior uncertainty\. We specifically propose structure\-preservingH​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}subspaces of conventional Raviart–Thomas andd​g​P0dgP\_\{0\}elements prescribed by a lightweight transformer\. Reduced\-order dynamics consistent with this subspace are learned by posing a conservation law in which a GP describes the fluxes between volumes\. This work hinges on a novel interface between mixed FEM spaces and GP regression; when training is posed as the optimal recovery problem \(ORP\), the resulting GP regression can be written as an optimization problem with equality constraints that impose a conservation structure, amenable to a fast Schur\-complement training strategy\. The trained model can then be solved in real time with closed\-form estimators for boundary fluxes driven by prescribed Dirichlet data\. The paper includes RKHS posterior error bounds for linear functionals to support uncertainty quantification, as well as numerical experiments demonstrating the accuracy of the posterior distribution as a surrogate for error estimation\.

###### keywords:

Scientific machine learning, Whitney forms, finite element exterior calculus, optimal recovery, uncertainty quantification

\{MSCcodes\}

65N30, 81Q30, 68T07, 60G15, 90C70

## 1Problem overview and relation to the literature

Neural operators and other machine learning surrogates are emerging as practical alternatives to classical simulation due to their computational efficiency and ability to generalize across problem instances\. However, the black\-box nature of these methods remains a major impediment to adoption in scientific and engineering applications and creates challenges for uncertainty quantification \(UQ\), where errors can propagate through the learned solution operator and further affect downstream predictions without clear traceability\[abdar2021review,NAP29212\]\. Moreover, for systems governed by partial differential equations \(PDEs\), the learned surrogate models should also respect underlying physical constraints, such as conservation laws and initial and boundary conditions\.

Motivated by this gap, we propose a surrogate framework for PDE\-governed systems that supports fast simulation while providing tractable posterior uncertainty\. Consider the following abstract problem, which we later elaborate in[Section4](https://arxiv.org/html/2606.11650#S4)\. Letu∈Qu\\in Qbe a field onΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}\(naturally,Q⊂L2​\(Ω\)Q\\subset L^\{2\}\(\\Omega\)\), assumed to satisfy an unknown conservation law∇⋅F​\(u\)=f\\nabla\\cdot F\(u\)=f, whereF∈VF\\in Vis a given flux function \(naturally,V⊂H​\(div;Ω\)V\\subset H\(\\mathrm\{div\};\\Omega\)\)\. We assume access to data consisting of sampled state and flux fieldsuuandFF\. Additionally, we allow for parametric dependence on a conditioning variableZZ, which may represent constitutive relationships, geometry, or other factors, yielding the training set

𝒟N=\{\(z\(n\),udata\(n\),Fdata\(n\)\)\}n=1N\.\\mathcal\{D\}\_\{N\}=\\left\\\{\\left\(z^\{\(n\)\},u^\{\(n\)\}\_\{\\text\{data\}\},F^\{\(n\)\}\_\{\\text\{data\}\}\\right\)\\right\\\}\_\{n=1\}^\{N\}\.
Our goal is to identify a functional form forFFthat is consistent with the available data, supports tractable posterior estimation, and generalizes to boundary data not seen during training\. This learning problem must be posed in an appropriate discrete setting so that the resulting surrogate preserves the structure of the underlying PDE\. We assume that the flux can be decomposed in the form

\(1\)𝐅=−ϵ​∇u\+𝒩​\[u\],∇⋅𝐅=f,\\mathbf\{F\}=\-\\epsilon\\nabla u\+\\mathcal\{N\}\[u\],\\qquad\\nabla\\cdot\\mathbf\{F\}=f,consisting of a diffusion term withϵ\>0\\epsilon\>0for numerical stability and a nonlinear correction𝒩\\mathcal\{N\}that identifies a state\-to\-flux map from data\. In previous work we demonstrated how neural operators can serve as𝒩\\mathcal\{N\}\[kinch2025structure\]; here we instead treat𝒩\\mathcal\{N\}as a Gaussian process \(GP\)\. Classical GPs are most tractable in small\-data, low\-dimensional settings, and a primary contribution of this work is that the reduced graph representation developed below brings the otherwise infinite\-dimensional state\-to\-flux learning problem into precisely such a regime\.

![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/roadmap.png)Figure 1:Roadmap to structure\-preserving surrogates with quantified uncertainty\.The transformer learnsH​\(div\)H\(\\mathrm\{div\}\)\-conforming bases and constructs a conservative coarse graph adapted to the conditioningZZ\(Boxes①\-③\)\. The graph edges carry GP models of the state\-to\-flux laws while conservation is imposed as an exact divergence constraint \(Box④\) and the resulting reduced model serves as a Dirichlet\-to\-Neumann surrogate with posterior error estimates for boundary fluxes \(Box⑤\)\.The proposed method separates the learning task into two parts: a data\-driven model mapping the fine\-scale space to a coarse\-scale space \(P1\), and a stochastic model for flux correction that supports uncertainty quantification \(P2\)\. Because the full construction is technical, we first summarize the salient features of each\.

H​\(g​r​a​d\)\{H\(grad\)\}H​\(c​u​r​l\)\{H\(curl\)\}H​\(d​i​v\)\{H\(div\)\}L2\{L^\{2\}\}Λ0\{\\Lambda\_\{0\}\}Λ1\{\\Lambda\_\{1\}\}Λ2\{\\Lambda\_\{2\}\}Λ3\{\\Lambda\_\{3\}\}Vhc\{V\_\{h\}^\{c\}\}Qhc\{Q\_\{h\}^\{c\}\}dddddd∇\\nabla∇×\\nabla\\times∇⋅\\nabla\\cdot∇⋅\\nabla\\cdotrVr\_\{V\}rQr\_\{Q\}H​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}subcomplexFigure 2:Construction ofH​\(div\)H\(\\mathrm\{div\}\)\-conforming reduced subspace\.We design a transformer that outputs a subspace of the de Rham complex appropriate for strongly imposing conservation laws\. While previous work builds “bottom\-up” dualΛ0/Λ1\\Lambda\_\{0\}/\\Lambda\_\{1\}reduced subcomplexes\[actor2024data\], this work provides a “top\-down” primal subcomplex onΛd/Λd−1\\Lambda\_\{d\}/\\Lambda\_\{d\-1\}\. A transformer prescribes the restriction maprQr\_\{Q\}conditioned on Z, and we provide a compatible construction ofrVr\_\{V\}that defines coarsenedR​T​0−d​g​P0RT0\-dgP\_\{0\}subspaces which preserve surjectivity of the divergence\.InP1, we work with low\-order Whitney forms\[arnold2010finite\]that provide conforming finite element spaces encoding the topological and cohomological properties tied to conservation structure\. Here the Raviart–Thomas and discontinuous\-piecewise\-constant \(RT0\\mathrm\{RT\}\_\{0\}/d​g​P0dgP\_\{0\}\) pair is naturally conforming, withRT0⊂H​\(div;Ω\)\\mathrm\{RT\}\_\{0\}\\subset H\(\\mathrm\{div\};\\Omega\)andd​g​P0⊂L2​\(Ω\)dgP\_\{0\}\\subset L^\{2\}\(\\Omega\), and has the surjectivity property∇⋅:RT0→dgP0\\nabla\\cdot:\\mathrm\{RT\}\_\{0\}\\rightarrow dgP\_\{0\}\. Because this space interpolates cell\-based scalar degrees of freedom and facet\-based flux degrees of freedom, it provides a discrete divergence theorem, allowing the discrete div/grad matrices to be interpreted as adjacency matrices between cells and facets\. This dual interpretation of div/grad is the crucial ingredient for Bayesian analysis on fields: in previous work\[owhadi2022computational,propp2026discovery\]we showed how to use optimal recovery to learn circuit models with tractable posteriors, and this linkage lets us extend the analysis to finite element fields\.

The goal is to identify reduced\-order spacesQhc​\(z;θ\)⊂QhfQ\_\{h\}^\{c\}\(z;\\theta\)\\subset Q^\{f\}\_\{h\}andVhc​\(z;θ\)⊂VhfV\_\{h\}^\{c\}\(z;\\theta\)\\subset V^\{f\}\_\{h\}, where the superscripts⋅f\\cdot^\{f\}and⋅c\\cdot^\{c\}denote fine and coarse spaces and the subscript⋅h\\cdot\_\{h\}denotes discretization\. To preserve exterior calculus structure in the reduced space, we require a construction in which the restrictionsrQ:Qhf→Qhcr\_\{Q\}:Q^\{f\}\_\{h\}\\rightarrow Q^\{c\}\_\{h\}andrV:Vhf→Vhcr\_\{V\}:V^\{f\}\_\{h\}\\rightarrow V^\{c\}\_\{h\}commute with the divergence operator, so that∇⋅Vhc⊆Qhc\\nabla\\cdot V\_\{h\}^\{c\}\\subseteq Q\_\{h\}^\{c\}\. To achieve this, we first design a transformer that evaluatesrQr\_\{Q\}by outputting a coarsening matrixWW, so thatQhc=span\{∑aWi​aχa\}iQ\_\{h\}^\{c\}=\\operatorname\{span\}\\left\\\{\\sum\_\{a\}W\_\{ia\}\\chi\_\{a\}\\right\\\}\_\{i\}, where theχa\\chi\_\{a\}are the indicator functions over the cells spanningQhfQ\_\{h\}^\{f\}\. We then designVhcV\_\{h\}^\{c\}by identifying coarsened degrees of freedom associated with coarse facets shared between coarse cells\. This novel construction represents a “top\-down” coarsening of the\(Λd/Λd−1\)\(\\Lambda\_\{d\}/\\Lambda\_\{d\-1\}\)subcomplex, in contrast to the “bottom\-up” coarsening of\(Λ0/Λ1\)\(\\Lambda\_\{0\}/\\Lambda\_\{1\}\)developed previously\[actor2024data\]\.P1therefore produces a reducedH​\(div\)H\(\\mathrm\{div\}\)\-conforming finite element space that can be conditioned onZZ\.

InP2, we pose the discrete representation of the physics as anoptimal recovery problem\[owhadi2022computational\]\. We identifyF∈VhcF\\in V\_\{h\}^\{c\}andu∈Qhcu\\in Q\_\{h\}^\{c\}by their basis coefficientsF^\\widehat\{F\}andu^\\widehat\{u\}, where each coarse flux degree of freedomF^i​j\\hat\{F\}\_\{ij\}is associated with the coarse oriented boundary shared by the coarse cellsu^i\\hat\{u\}\_\{i\}andu^j\\hat\{u\}\_\{j\}in the graph induced byP1\. We model the nonlinearity edgewise by a Gaussian process,𝒩=𝒢​𝒫​\(ui​j\)\\mathcal\{N\}=\\mathcal\{GP\}\(u\_\{ij\}\), whereui​ju\_\{ij\}concatenates the statesu^i,u^j\\hat\{u\}\_\{i\},\\hat\{u\}\_\{j\}with the additional metric features needed to generalize the flux law across the mesh\.

In traditional GP regression, the posterior distribution is derived from the joint normal distribution via a Schur complement\. Its mean may equivalently be obtained by solving the optimal recovery problem: given a kernelKKwith reproducing kernel Hilbert space \(RKHS\)ℋK\\mathcal\{H\}\_\{K\}andNNnoisy data pairs\(𝐗,𝐘\)=\{\(𝐱i,yi\)\}i=1N\(\\mathbf\{X\},\\mathbf\{Y\}\)=\\\{\(\\mathbf\{x\}\_\{i\},y\_\{i\}\)\\\}\_\{i=1\}^\{N\}with noise varianceσε2\\sigma\_\{\\varepsilon\}^\{2\}, the minimizer of the RKHS norm penalized by data misfit

\(2\)f^=arg​minf∈ℋK⁡‖f‖ℋK2\+1σε2​∑i=1N\(f​\(𝐱i\)−yi\)2,\\hat\{f\}\\;=\\;\\operatorname\*\{arg\\,min\}\_\{f\\in\\mathcal\{H\}\_\{K\}\}\\;\\\|f\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}\\;\+\\;\\frac\{1\}\{\\sigma\_\{\\varepsilon\}^\{2\}\}\\sum\_\{i=1\}^\{N\}\\bigl\(f\(\\mathbf\{x\}\_\{i\}\)\-y\_\{i\}\\bigr\)^\{2\},recovers the conventional GP estimatorf^​\(⋅\)=K​\(⋅,𝐗\)​\(K​\(𝐗,𝐗\)\+σε2​I\)−1​𝐘\\hat\{f\}\(\\cdot\)=K\(\\cdot,\\mathbf\{X\}\)\\bigl\(K\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\\bigr\)^\{\-1\}\\mathbf\{Y\}\.

Recasting the conventional GP problem in this manner allows us to incorporate equality constraints directly\. In the graph interpretation ofP1, fluxes overVhcV\_\{h\}^\{c\}are edge currents between nodes associated with states inQhcQ\_\{h\}^\{c\}\. On each edge we cast an optimal recovery problem mapping state to flux, subject to the equality constraint that the conservation law holds at each node\. The connectivity of the coarse bases imposes sparsity on this dense graph through a metric weighting, and after training we are left with a discrete boundary value problem whose boundary fluxes are encoded by GPs\. Because this is a constrained optimization problem, its KKT conditions expose a saddle\-point structure that we exploit to derive a fast optimizer \([Section4\.2](https://arxiv.org/html/2606.11650#S4.SS2)\)\.

In concert,P1andP2amount to agraph discovery problem: by concurrently training the transformer and solving the optimal recovery problem, we interpret the identification of reduced div/grad finite element spaces as the identification of a dense graph, sparsified by the connectivity of the bases, that encodes the flow of conserved quantities through the system \([Figure1](https://arxiv.org/html/2606.11650#S1.F1)\)\. While circuit analogies are commonplace in engineering \(e\.g\., compact models for semiconductor devices\[pmlr\-v107\-aadithya20a,fan2023two\], hydraulic circuits in hydrodynamics\[9385620,vacca2021hydraulic\], or thermal circuits for heat transfer\[wang2017microscale\]\), they are typically built from simplified analytic solutions and empirical curve\-fitting through a slow iterative process\. In contrast, our approach can be carried out autonomously to obtain rapid, uncertainty\-quantified surrogates that preserve a conservative input/output relationship\. We demonstrate this by considering two representative examples\. Advection\-diffusion on a triangular mesh of Philadelphia’s Liberty Bell serves to illustrate how real\-time surrogates with quantified uncertainty can be constructed on complex geometries; in this setting, we condition on the direction of advection as an example of how to construct parametric models\. We finally construct a digital twin of a semiconductor device \(specifically, ap−np\-ndiode\)\. Training data can be constructed by TCAD simulations solving the device drift\-diffusion equations\[musson2022charon\], providing examples of internal device transport and how it defines the voltage\-current relationship governing the device\. This problem specifically shows how the UQ developed here identifies the range of inputs over which the learned surrogate can be trustworthy\.

### 1\.1Relation to literature

The construction above brings together three goals usually pursued separately: fast reduced\-order simulation, exact enforcement of physical structure, and uncertainty quantification\. We review each in turn\. Classical discretizations such as finite element methods yield accurate, verifiable predictions but at substantial cost, motivating a large body of scientific machine learning \(SciML\) that trades rigor for speed\. Physics\-informed neural networks \(PINNs\) embed the governing PDEs in the training loss\[karniadakis2021physics,yu2022gradient\]; neural operators such as DeepONet\[deeponetNatureML,pideeponet\]and Fourier neural operators\[li2020fourier\]learn maps between function spaces; and data\-driven reduced\-order models \(ROMs\) accelerate parametric simulation by working in a low\-dimensional subspace\[fresca2022pod,jung2025accelerating,kapteyn2022data\]\. These surrogates can be fast, but typically forfeit the structural guarantees and error control of the solvers they replace\. A related line of structure\-preserving model\-reduction methods addresses part of this gap by seeking reduced spaces that retain conservation, stability, passivity, or compatibility properties of the full\-order discretization\[benner2015survey,carlberg2018conservative,quarteroni2015reduced\]\. Our construction is close in spirit to this line of work, but learns the compatible reducedH​\(div\)−L2H\(\\mathrm\{div\}\)\-L^\{2\}complex from data rather than selecting it through a fixed projection, and couples it to a GP optimal\-recovery formulation for posterior uncertainty\.

Enforcing physical structure reliably is the central difficulty, particularly under limited data and complex geometry\. Most physics\-informed approaches impose conservation through soft penalties on the PDE residual or boundary conditions\[chen2024physics,jiao2024solving\]; these do not guarantee exact conservation, and the resulting violations can accumulate, amplify errors, and destabilize training on stiff or multiscale problems\[bonfanti2024challenges,krishnapriyan2021characterizing\]\. Structure\-preserving methods instead build the constraints into the function spaces themselves\. The mixed finite element theory underlying Raviart\-Thomas and discontinuous Galerkin pairs provides the classical foundation for suchH​\(div\)H\(\\mathrm\{div\}\)\-conforming flux approximations and compatibleL2L^\{2\}scalar spaces; see, for example,\[boffi2013mixed\]\. Finite element exterior calculus \(FEEC\) provides a canonical framework for this approach, choosing spaces and operators that form a discrete de Rham complex and so inherit the topological identities underlying conservation\[arnold2010finite,Arnold\_Falk\_Winther\_2006\]\. Data\-driven exterior calculus extends this to learned, graph\-based operators while preserving exact\-sequence compatibility\[trask2022enforcing\]\. Our coarsenedRT0\\mathrm\{RT\}\_\{0\}/d​g​P0dgP\_\{0\}complex is designed to inherit these guarantees by construction\.

Finally, trustworthy surrogates require calibrated uncertainty estimates that can be propagated to downstream quantities of interest\[abdar2021review,psaros2023uncertainty,xu2021accurate\]\. UQ for neural PDE surrogates is often pursued through ensembles, Bayesian neural networks, dropout, latent\-variable models, or post\-hoc calibration; these approaches can be effective but usually do not yield the closed\-form functional posterior bounds available in GP/RKHS optimal recovery\[song2026structure,yang2019adversarial\]\. Gaussian processes provide a natural alternative because conditioning a GP prior on data yields a posterior mean and covariance in closed\-form\. The same estimator arises from the optimal recovery problem, which characterizes the minimum\-norm RKHS interpolant of noisy data\[chen2021solving,owhadi2022computational,propp2026discovery\]and, as we exploit below, accommodates hard linear constraints\. We build most directly on the graph\-based optimal recovery of Dirichlet\-to\-Neumann maps in\[propp2026discovery\], extending it from fixed graphs to the learned, conditioned finite element complexes produced byP1\.

### 1\.2Main contributions

[Figure1](https://arxiv.org/html/2606.11650#S1.F1)summarizes the resulting framework with structure preservation and uncertainty quantification\. Our specific contributions are:

- •a*top\-down*coarsening of the\(Λd/Λd−1\)\(\\Lambda\_\{d\}/\\Lambda\_\{d\-1\}\)subcomplex, in which a lightweight transformer emits a coarsening matrixWWthat defines a reduced,H​\(div\)H\(\\mathrm\{div\}\)\-conformingRT0\\mathrm\{RT\}\_\{0\}/d​g​P0dgP\_\{0\}pair conditioned on the problem parametersZZ, in contrast to the bottom\-up\(Λ0/Λ1\)\(\\Lambda\_\{0\}/\\Lambda\_\{1\}\)coarsening of\[actor2024data\];
- •a structure\-preserving optimal recovery formulation in which an edgewise Gaussian process models the state\-to\-flux law and conservation is imposed exactly through linear equality constraints, yielding a saddle\-point KKT system with a fast Schur\-complement solve;
- •a joint “graph discovery” training procedure that learns the reduced complex and the flux Gaussian process concurrently, producing a sparse circuit model of the conserved dynamics that generalizes across geometries and boundary data;
- •a closed\-form RKHS posterior error bound for linear functionals of the flux, in particular the boundary fluxes defining the Dirichlet\-to\-Neumann map, validated on complex geometries and semiconductor device equations\.

The remainder of the paper is organized as follows\.[Section2](https://arxiv.org/html/2606.11650#S2)reviews the necessary mathematical background for the proposed framework\.[Section2\.1](https://arxiv.org/html/2606.11650#S2.SS1)reviews the finite element exterior calculus preliminaries for theH​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}mixed setting and[Section2\.2](https://arxiv.org/html/2606.11650#S2.SS2)specifies the finite element space with conformingRT0\\mathrm\{RT\}\_\{0\}/d​g​P0dgP\_\{0\}discretizations\.[Section2\.3](https://arxiv.org/html/2606.11650#S2.SS3)explains the optimal recovery problem for Gaussian processes\.[Section3](https://arxiv.org/html/2606.11650#S3)develops the first component \(P1\) by introducing the data\-driven PoU Whitney construction \([Section3\.1](https://arxiv.org/html/2606.11650#S3.SS1)\) and the induced coarse operators \([Section3\.2](https://arxiv.org/html/2606.11650#S3.SS2)\) for structure preservation\. We show that under this construction, properties such as surjectivity and flux balance still hold in the reduced space \([Section3\.3](https://arxiv.org/html/2606.11650#S3.SS3)\)\. In addition, we provide a graph interpretation of the learned reduced spaces \([Section3\.4](https://arxiv.org/html/2606.11650#S3.SS4)\), which supports the subsequent optimal recovery formulation\.[Section4](https://arxiv.org/html/2606.11650#S4)develops the second component \(P2\)\. We cast the GP optimal recovery formulation on the coarse spaces subject to equality constraints that guarantee the exact enforcement of conservation laws \([Section4\.1](https://arxiv.org/html/2606.11650#S4.SS1)\)\. In addition to a fast Schur\-complement solve for the KKT system, we propose a bilevel training procedure for efficient training with nested coarse variables and models \([Section4\.2](https://arxiv.org/html/2606.11650#S4.SS2)\) and incorporate geometric embedding \([Section4\.3](https://arxiv.org/html/2606.11650#S4.SS3)\)\. We also provide the theoretical analysis for the RKHS posterior error bound in[Section4\.4](https://arxiv.org/html/2606.11650#S4.SS4)\.[Section5](https://arxiv.org/html/2606.11650#S5)reports numerical experiments on representative examples, a complex\-geometry advection\-diffusion problem, and a semiconductor diode example to validate the proposed framework\.

## 2Mathematical preliminaries

In this section we gather the necessary background in both finite element exterior calculus \(FEEC\) and the optimal recovery formalisms\. For further background, please see\[arnold2010finite,Arnold\_Falk\_Winther\_2006,micchelli1977survey\]and\[trask2022enforcing\]\.

### 2\.1Finite element exterior calculus \(FEEC\)

Exterior calculus generalizes classical calculus to differential forms of higher degree on differentiable manifolds and provides a robust mathematical framework for analyzing PDEs on different geometries\. Finite element exterior calculus extends the exterior calculus to discrete finite element meshes using chain complexes and preserves crucial topological and geometric structures underlying the governing PDEs\[kinch2025structure,trask2022enforcing\]\.

One advantage of FEEC is that it unifies the main differential operators, i\.e\., grad, div, and curl, with the de Rham complex\. LetΛ​\(Ω\)\\Lambda\(\\Omega\)denote the exterior algebra with exterior product∧\\wedgeandΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}\. We then have the de Rham complex:

\(3\)0→Λ0​\(Ω\)→d0Λ1​\(Ω\)→d1⋯→dk−1Λk​\(Ω\)→0,0\\rightarrow\\Lambda^\{0\}\(\\Omega\)\\xrightarrow\{d^\{0\}\}\\Lambda^\{1\}\(\\Omega\)\\xrightarrow\{d^\{1\}\}\\cdots\\xrightarrow\{d^\{k\-1\}\}\\Lambda^\{k\}\(\\Omega\)\\rightarrow 0,whered:Λk−1​\(Ω\)→Λk​\(Ω\)d:\\Lambda^\{k\-1\}\(\\Omega\)\\rightarrow\\Lambda^\{k\}\(\\Omega\)is the exterior derivative operator\. ForΩ⊂ℝ3\\Omega\\subset\\mathbb\{R\}^\{3\}, the de Rham complex becomes:

\(4\)0→C∞​\(Ω\)→grad\[C∞​\(Ω\)\]3→curl\[C∞​\(Ω\)\]3→divC∞​\(Ω\)→0\.0\\rightarrow C^\{\\infty\}\(\\Omega\)\\xrightarrow\{\\text\{grad\}\}\[C^\{\\infty\}\(\\Omega\)\]^\{3\}\\xrightarrow\{\\text\{curl\}\}\[C^\{\\infty\}\(\\Omega\)\]^\{3\}\\xrightarrow\{\\text\{div\}\}C^\{\\infty\}\(\\Omega\)\\rightarrow 0\.If we further consider Sobolev spaces with homogeneous Dirichlet boundary conditions, we can derive the FEEC specialization primal and dual cochain complex:

\(5\)0⟶H0​\(grad,Ω\)⇌−gradgradH0​\(curl,Ω\)⇌curlcurlH0​\(div,Ω\)⇌−divdivL2​\(Ω\)⟶0\.0\\longrightarrow H\_\{0\}\(\\mathrm\{grad\},\\Omega\)\\xrightleftharpoons\[\\,\-\\mathrm\{grad\}\\,\]\{\\ \\mathrm\{grad\}\\ \}H\_\{0\}\(\\mathrm\{curl\},\\Omega\)\\xrightleftharpoons\[\\,\\mathrm\{curl\}\\,\]\{\\ \\mathrm\{curl\}\\ \}H\_\{0\}\(\\mathrm\{div\},\\Omega\)\\xrightleftharpoons\[\\,\-\\mathrm\{div\}\\,\]\{\\ \\mathrm\{div\}\\ \}L^\{2\}\(\\Omega\)\\longrightarrow 0\.In this work, instead of working on “bottom\-up” subcomplexes, we work on the terminal mapH​\(div,Ω\)→∇⋅L2​\(Ω\)H\(\\mathrm\{div\},\\Omega\)\\xrightarrow\{\\nabla\\cdot\}L^\{2\}\(\\Omega\)because this is a natural algebraic mechanism through which flux balance is imposed\. This is illustrated in the commuting diagram in[Figure2](https://arxiv.org/html/2606.11650#S1.F2)\.

### 2\.2H​\(div\)H\(\\mathrm\{div\}\)space and lowest\-order Raviart–Thomas element

Since one core objective of this work is to guarantee the flux conservation law, instead of working on the entire de Rham complex, we focus on its tail, i\.e\.,H​\(div,Ω\)→divL2​\(Ω\)H\(\\mathrm\{div\},\\Omega\)\\xrightarrow\{\\text\{div\}\}L^\{2\}\(\\Omega\)\. This segment of the complex supports flux conservation in a natural way through the mixed formulation of elliptic problems\. To make this connection explicit, consider the general second\-order elliptic equation in a bounded domainΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}:

\(6\)−∇⋅\(𝐊​\(x\)​∇u\)=fin​Ω,\-\\nabla\\cdot\\bigl\(\\mathbf\{K\}\(x\)\\nabla u\\bigr\)=f\\qquad\\text\{in \}\\Omega,where𝐊​\(x\)∈ℝd×d\\mathbf\{K\}\(x\)\\in\\mathbb\{R\}^\{d\\times d\}is a symmetric and uniformly positive definite diffusion tensor\. The boundary∂Ω\\partial\\Omegais decomposed into disjoint parts∂Ω=ΓD∪ΓN\\partial\\Omega=\\Gamma\_\{D\}\\cup\\Gamma\_\{N\}on which Dirichlet and Neumann boundary conditions are defined, respectively\. Then by introducing the flux variable𝐅≔−𝐊​\(x\)​∇u\\mathbf\{F\}\\coloneqq\-\\mathbf\{K\}\(x\)\\nabla u, the problem can be rewritten as the first\-order system:

\(7\)\{𝐊​\(x\)−1​𝐅\+∇u=0,∇⋅𝐅=f,in​Ω\.\\begin\{cases\}\\mathbf\{K\}\(x\)^\{\-1\}\\mathbf\{F\}\+\\nabla u&=0,\\\\ \\nabla\\cdot\\mathbf\{F\}&=f,\\end\{cases\}\\qquad\\text\{in \}\\Omega\.In this setting, the natural space for the flux variable is:

\(8\)H​\(div,Ω\)≔\{𝐯∈\[L2​\(Ω\)\]d\|∇⋅𝐯∈L2​\(Ω\)\}\.H\(\\mathrm\{div\},\\Omega\)\\coloneqq\\left\\\{\\mathbf\{v\}\\in\[L^\{2\}\(\\Omega\)\]^\{d\}\\;\\middle\|\\;\\nabla\\cdot\\mathbf\{v\}\\in L^\{2\}\(\\Omega\)\\right\\\}\.This space ensures that the divergence operator is well\-defined in a weak sense and the normal components𝐯⋅𝐧\\mathbf\{v\}\\cdot\\mathbf\{n\}are well\-defined on element interfaces, thus naturally making it the appropriate space for modeling local flux conservation\. The corresponding mixed variational formulation is to find\(𝐅,u\)∈H​\(div,Ω\)×L2​\(Ω\)\(\\mathbf\{F\},u\)\\in H\(\\mathrm\{div\},\\Omega\)\\times L^\{2\}\(\\Omega\)such that

\(9\)∫Ω𝐊​\(x\)−1​𝐅⋅𝐯​𝑑x−∫Ωu​\(∇⋅𝐯\)​𝑑x=−∫ΓDuD​\(𝐯⋅𝐧\)​𝑑s,∀𝐯∈H​\(div,Ω\),∫Ω\(∇⋅𝐅\)​q​𝑑x=∫Ωf​q​𝑑x,∀q∈L2​\(Ω\)\.\\begin\{split\}\\int\_\{\\Omega\}\\mathbf\{K\}\(x\)^\{\-1\}\\mathbf\{F\}\\cdot\\mathbf\{v\}\\,dx\-\\int\_\{\\Omega\}u\\,\(\\nabla\\cdot\\mathbf\{v\}\)\\,dx&=\-\\int\_\{\\Gamma\_\{D\}\}u\_\{D\}\\,\(\\mathbf\{v\}\\cdot\\mathbf\{n\}\)\\,ds,\\quad\\forall\\,\\mathbf\{v\}\\in H\(\\mathrm\{div\},\\Omega\),\\\\ \\int\_\{\\Omega\}\(\\nabla\\cdot\\mathbf\{F\}\)\\,q\\,dx&=\\int\_\{\\Omega\}f\\,q\\,dx,\\quad\\forall\\,q\\in L^\{2\}\(\\Omega\)\.\\end\{split\}A conforming discretization of this mixed formulation requires a finite element flux space inH​\(div,Ω\)H\(\\mathrm\{div\},\\Omega\)and a compatible scalar space inL2​\(Ω\)L^\{2\}\(\\Omega\)\. In this work, we consider the Raviart–Thomas \(RT\) elements for the flux variable and piecewise constants for the scalar variable\. Let𝒯h\\mathcal\{T\}\_\{h\}be a simplicial mesh ofΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}\. On each elementK∈𝒯hK\\in\\mathcal\{T\}\_\{h\}, the Raviart–Thomas space of orderk≥0k\\geq 0is defined by

\(10\)RTk​\(K\)≔Pk​\(K\)d\+𝐱​Pk​\(K\),\\mathrm\{RT\}\_\{k\}\(K\)\\coloneqq P\_\{k\}\(K\)^\{d\}\+\\mathbf\{x\}\\,P\_\{k\}\(K\),wherePk​\(K\)\{P\}\_\{k\}\(K\)denotes the space of polynomials of degree at mostkkonKK\. The global Raviart–Thomas space is defined as

\(11\)RTk​\(𝒯h\)≔\{𝐯∈H​\(div,Ω\)\|𝐯\|K∈RTk​\(K\),∀K∈𝒯h\},\\mathrm\{RT\}\_\{k\}\(\\mathcal\{T\}\_\{h\}\)\\coloneqq\\left\\\{\\mathbf\{v\}\\in H\(\\mathrm\{div\},\\Omega\)\\;\\middle\|\\;\\mathbf\{v\}\|\_\{K\}\\in\\mathrm\{RT\}\_\{k\}\(K\),\\;\\forall K\\in\\mathcal\{T\}\_\{h\}\\right\\\},which enforces the continuity of normal components across element interfaces\.

Fork=0k=0, the*lowest–order Raviart–Thomas space*onKKis defined as

\(12\)RT0​\(K\)=P0​\(K\)d\+x​P0​\(K\)\.\\text\{RT\}\_\{0\}\(K\)\\;=\\;P\_\{0\}\(K\)^\{d\}\\;\+\\;x\\,P\_\{0\}\(K\)\.Equivalently, any𝐯∈RT0​\(K\)\\mathbf\{v\}\\in\\mathrm\{RT\}\_\{0\}\(K\)can be written as𝐯​\(x\)=𝐚\+p​x\\mathbf\{v\}\(x\)=\\mathbf\{a\}\+p\\,xwith constant vector𝐚∈ℝd\\mathbf\{a\}\\in\\mathbb\{R\}^\{d\}and constant scalarp∈ℝp\\in\\mathbb\{R\}\. The degrees of freedom of this lowest\-order RT0space are defined by the normal fluxes across element faces:

\(13\)∫e\(𝐯⋅𝐧\)​q​𝑑s,∀e⊂∂K\.\\int\_\{e\}\(\\mathbf\{v\}\\cdot\\mathbf\{n\}\)\\,q\\,ds,\\qquad\\forall e\\subset\\partial K\.Thus the degree of freedom admits a direct interpretation in terms of cochains\. In particular, the RT space ensures that the discrete divergence operator is compatible with the underlying differential complex, and the\(RT0,P0\)\(\\mathrm\{RT\}\_\{0\},\\mathrm\{P\}\_\{0\}\)can form a meaningful discretization of theH​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}mixed formulation\. Within the FEEC framework, this construction corresponds to the final part of the de Rham complex where the scalar field is associated withdd\-forms and the flux corresponds to\(d−1\)\(d\-1\)\-forms\. The gradient and divergence operators are unified through the exterior derivative and the codifferential\. The resulting compatible structure provides the algebraic foundation for conservation and motivates the learnable coarse space construction below\.

### 2\.3Optimal recovery problem for Gaussian processes

Gaussian processes provide a probabilistic model for unknown functions\. Formally, a GP is a collection of random variables such that any finite collection has a joint Gaussian distribution\. Equivalently, a GP can be viewed as a distribution over functions: each draw from the process is a possible function, and for any finite set of input locations, the corresponding vector of function values is multivariate Gaussian\. This makes GPs particularly useful in settings where the uncertainty quantification is important, since conditioning on observed data yields both a posterior mean prediction and a posterior covariance that quantifies uncertainty at new inputs\.

Given data𝒟=\(𝐗,𝐘\)\\mathcal\{D\}=\(\\mathbf\{X\},\\mathbf\{Y\}\)withX=\[𝐱1,…,𝐱N\]⊤X=\[\\mathbf\{x\}\_\{1\},\\dots,\\mathbf\{x\}\_\{N\}\]^\{\\top\}and corresponding output𝐘=\[y1,…,yN\]⊤\\mathbf\{Y\}=\[y\_\{1\},\\dots,y\_\{N\}\]^\{\\top\}generated by an underlying functionyi=f​\(𝐱i\)y\_\{i\}=f\(\\mathbf\{x\}\_\{i\}\), a GP prior onffis fully specified by its mean functionm​\(x\)m\(x\)and covariance kernelK​\(x,x′\)K\(x,x^\{\\prime\}\):

\(14\)m​\(𝐱\)\\displaystyle m\(\\mathbf\{x\}\)=𝔼​\[f​\(𝐱\)\]\\displaystyle=\\mathbb\{E\}\[f\(\\mathbf\{x\}\)\]\(15\)K​\(𝐱,𝐱′\)\\displaystyle K\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)=𝔼​\[\(f​\(𝐱\)−m​\(𝐱\)\)​\(f​\(𝐱′\)−m​\(𝐱′\)\)\]\\displaystyle=\\mathbb\{E\}\[\(f\(\\mathbf\{x\}\)\-m\(\\mathbf\{x\}\)\)\(f\(\\mathbf\{x\}^\{\\prime\}\)\-m\(\\mathbf\{x\}^\{\\prime\}\)\)\]\(16\)f​\(𝐱\)\\displaystyle f\(\\mathbf\{x\}\)∼𝒢​𝒫​\(m​\(𝐱\),K​\(𝐱,𝐱′\)\)\\displaystyle\\sim\\mathcal\{GP\}\(m\(\\mathbf\{x\}\),K\(\\mathbf\{x\},\\mathbf\{x\}^\{\\prime\}\)\)In the noise\-free scenario, for test inputs𝐗∗=\(𝐱1∗,…,𝐱N∗∗\)T\\mathbf\{X\}\_\{\*\}=\(\\mathbf\{x\}\_\{1\}^\{\*\},\\ldots,\\mathbf\{x\}\_\{N\_\{\*\}\}^\{\*\}\)^\{T\}, the joint distribution of the observed training values𝐘=f​\(𝐗\)\\mathbf\{Y\}=f\(\\mathbf\{X\}\)and the unobserved test values𝐟∗=f​\(X∗\)\\mathbf\{f\_\{\*\}\}=f\(X\_\{\*\}\)is111In \([17](https://arxiv.org/html/2606.11650#S2.E17)\), we use a zero\-mean GP prior for simplicity\. This is a standard convention when the data have been centered or when the GP models a residual correction around a deterministic mean model\. A nonzero mean can be incorporated by applying the same formulas to the centered quantities obtained after subtracting the mean\. Under this zero\-mean prior, the joint distribution of the observed training outputs and the unobserved test values is also zero mean\.

\(17\)\[𝐘𝐟∗\]∼𝒩​\(𝟎,\[K​\(𝐗,𝐗\)K​\(𝐗,𝐗∗\)K​\(𝐗∗,𝐗\)K​\(𝐗∗,𝐗∗\)\]\)\.\\begin\{bmatrix\}\\mathbf\{Y\}\\\\ \\mathbf\{f\}\_\{\*\}\\end\{bmatrix\}\\sim\\mathcal\{N\}\\left\(\\mathbf\{0\},\\begin\{bmatrix\}K\(\\mathbf\{X\},\\mathbf\{X\}\)&K\(\\mathbf\{X\},\\mathbf\{X\}\_\{\*\}\)\\\\ K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)&K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{\*\}\)\\end\{bmatrix\}\\right\)\.Conditioning on this joint Gaussian distribution gives

\(18\)𝐟∗\|𝐗∗,𝐗,𝐘∼𝒩​\(μ∗,Σ∗\),\\mathbf\{f\}\_\{\*\}\|\\mathbf\{X\}\_\{\*\},\\mathbf\{X\},\\mathbf\{Y\}\\sim\\mathcal\{N\}\(\\mu\_\{\*\},\\Sigma\_\{\*\}\),where

\(19\)μ∗\\displaystyle\\mu\_\{\*\}=K​\(𝐗∗,𝐗\)​K​\(𝐗,𝐗\)−1​𝐘,\\displaystyle=K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)K\(\\mathbf\{X\},\\mathbf\{X\}\)^\{\-1\}\\mathbf\{Y\},\(20\)Σ∗\\displaystyle\\Sigma\_\{\*\}=K​\(𝐗∗,𝐗∗\)−K​\(𝐗∗,𝐗\)​K​\(𝐗,𝐗\)−1​K​\(𝐗,𝐗∗\)\.\\displaystyle=K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{\*\}\)\-K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)K\(\\mathbf\{X\},\\mathbf\{X\}\)^\{\-1\}K\(\\mathbf\{X\},\\mathbf\{X\}\_\{\*\}\)\.In most applications, however, the data are noisy\. Suppose𝐘=f​\(𝐗\)\+ε\\mathbf\{Y\}=f\(\\mathbf\{X\}\)\+\\varepsilon, whereε∼𝒩​\(0,σε2​I\)\\varepsilon\\sim\\mathcal\{N\}\(0,\\sigma\_\{\\varepsilon\}^\{2\}I\)is independent identically distributed \(i\.i\.d\.\) Gaussian noise\. Then the joint distribution in \([17](https://arxiv.org/html/2606.11650#S2.E17)\)\-\([18](https://arxiv.org/html/2606.11650#S2.E18)\) becomes

\(21\)\[𝐘𝐟∗\]\\displaystyle\\begin\{bmatrix\}\\mathbf\{Y\}\\\\ \\mathbf\{f\}\_\{\*\}\\end\{bmatrix\}∼𝒩​\(𝟎,\[K​\(𝐗,𝐗\)\+σε2​IK​\(𝐗,𝐗∗\)K​\(𝐗∗,𝐗\)K​\(𝐗∗,𝐗∗\)\]\),\\displaystyle\\sim\\mathcal\{N\}\\left\(\\mathbf\{0\},\\begin\{bmatrix\}K\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\sigma\_\{\\varepsilon\}^\{2\}I&K\(\\mathbf\{X\},\\mathbf\{X\}\_\{\*\}\)\\\\ K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)&K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{\*\}\)\\end\{bmatrix\}\\right\),𝐟∗\|𝐗∗,𝐗,𝐘\\displaystyle\\mathbf\{f\}\_\{\*\}\|\\mathbf\{X\}\_\{\*\},\\mathbf\{X\},\\mathbf\{Y\}∼𝒩\(K\(𝐗∗,𝐗\)\(K\(𝐗,𝐗\)\+σε2I\)−1𝐘,\\displaystyle\\sim\\mathcal\{N\}\\Big\(K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)\(K\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\)^\{\-1\}\\mathbf\{Y\},\(22\)K\(𝐗∗,𝐗∗\)−K\(𝐗∗,𝐗\)\(K\(𝐗,𝐗\)\+σε2I\)−1K\(𝐗,𝐗∗\)\)\.\\displaystyle\\qquad\\quad K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\_\{\*\}\)\-K\(\\mathbf\{X\}\_\{\*\},\\mathbf\{X\}\)\(K\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\)^\{\-1\}K\(\\mathbf\{X\},\\mathbf\{X\}\_\{\*\}\)\\Big\)\.The posterior covariance provides a meaningful probabilistic measure of uncertainty at the test points, which is useful for uncertainty\-aware scientific machine learning tasks\. Indeed, prior work has established that incorporating GPs within SciML frameworks provides rigorous uncertainty quantification for learning complex PDE\-governed systems\[brunton2024promising,chen2021solving,harkonen2023gaussian,shi2025survey\]\.

Given limited noisy data and some prior knowledge of the function space, the optimal recovery problem aims to learn an unknown function by minimizing the worst\-case error\. Mathematically, given a normed spaceFFand an unknown functionf∈Ff\\in F, the partial knowledge offfis defined through point evaluations:

\(23\)yi=ℓi​\(f\),i=1,…,m,\\displaystyle y\_\{i\}=\\ell\_\{i\}\(f\),\\quad i=1,\\dots,m,for some linear functionalsℓ1,…,ℓm∈F∗\\ell\_\{1\},\\dots,\\ell\_\{m\}\\in F^\{\*\}, whereF∗F^\{\*\}is the dual space ofFF\. The objective is to approximateffby some estimatef^∈F\\widehat\{f\}\\in Fsuch that the worst\-case error is minimized:

ℰ​\(f^\)≔inff^∈Fsupf∈Fℓi​\(f\)=yi‖f−f^‖‖f‖\.\\mathcal\{E\}\(\\widehat\{f\}\)\\coloneqq\\inf\_\{\\widehat\{f\}\\in F\}\\sup\_\{\\begin\{subarray\}\{c\}f\\in F\\\\ \\ell\_\{i\}\(f\)=y\_\{i\}\\end\{subarray\}\}\\frac\{\\\|f\-\\widehat\{f\}\\\|\}\{\\\|f\\\|\}\.LetK:𝒳×𝒳→ℝK:\\mathcal\{X\}\\times\\mathcal\{X\}\\rightarrow\\mathbb\{R\}be a symmetric positive definite bivariate kernel, andℋK\\mathcal\{H\}\_\{K\}denote its reproducing kernel Hilbert space \(RKHS\) with accompanying induced RKHS norm∥⋅∥K\\\|\\cdot\\\|\_\{K\}\. Given noisy observations\(𝐗,𝐘\)\(\\mathbf\{X\},\\mathbf\{Y\}\), the corresponding regularized optimal recovery problem is to find the functionf^∈ℋK\\hat\{f\}\\in\\mathcal\{H\}\_\{K\}that minimizes the sum of the squared RKHS norm and the scaled data\-misfit term:

𝒥​\(𝐗;𝐘\)=ming∈ℋK⁡‖g‖K2\+1σε2​‖g​\(𝐗\)−𝐘‖22\.\\mathcal\{J\}\(\\mathbf\{X\};\\mathbf\{Y\}\)=\\min\_\{g\\in\\mathcal\{H\}\_\{K\}\}\\\|g\\\|\_\{K\}^\{2\}\+\\frac\{1\}\{\\sigma\_\{\\varepsilon\}^\{2\}\}\\\|g\(\\mathbf\{X\}\)\-\\mathbf\{Y\}\\\|\_\{2\}^\{2\}\.By the Representer Theorem, the minimizerf^\\hat\{f\}has a finite\-dimensional representation in terms of kernel evaluations at the training inputs\. In particular,f^\\hat\{f\}evaluated at a new point can be written as

f^​\(⋅\)=K​\(⋅,𝐗\)​\(K​\(𝐗,𝐗\)\+σε2​I\)−1​𝐘\.\\hat\{f\}\(\\cdot\)=K\(\\cdot,\\mathbf\{X\}\)\(K\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\)^\{\-1\}\\mathbf\{Y\}\.This estimator coincides with the posterior mean of Gaussian process regression with covariance kernelKKand independent Gaussian observation noise of varianceσε2\\sigma\_\{\\varepsilon\}^\{2\}\. The corresponding minimum value of the objective𝒥​\(𝐗;𝐘\)\\mathcal\{J\}\(\\mathbf\{X\};\\mathbf\{Y\}\)is

𝒥​\(𝐗;𝐘\)=𝐘⊤​\(K​\(𝐗,𝐗\)\+σε2​I\)−1​𝐘,\\mathcal\{J\}\(\\mathbf\{X\};\\mathbf\{Y\}\)=\\mathbf\{Y\}^\{\\top\}\\left\(K\(\\mathbf\{X\},\\mathbf\{X\}\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\\right\)^\{\-1\}\\mathbf\{Y\},proved in\[propp2026discovery\]\.

## 3Method part 1: Data\-drivenH​\(div\)H\(\\mathrm\{div\}\)\-conforming reduced space construction

In this section, we construct the structure\-preserving reduced spaces used by the surrogate\. The goal is to learn the data\-driven Whitney forms from a trainable partition of unity \(PoU\) conditioned onZZ\. The PoU functions define coarse control volumes, represented by 0\-forms, while the associated 1\-forms encode generalized fluxes between them \([Section3\.1](https://arxiv.org/html/2606.11650#S3.SS1)\)\. This construction yields reduced scalar and flux spaces compatible with theH​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}structure introduced above, allowing the resultingH​\(div\)H\(\\mathrm\{div\}\)\-conforming finite element spaces to adapt to the conditioning variableZZ\([Sections3\.2](https://arxiv.org/html/2606.11650#S3.SS2)and[3\.3](https://arxiv.org/html/2606.11650#S3.SS3)\)\. Finally, we connect the induced coarse spaces to a graph interpretation in[Section3\.4](https://arxiv.org/html/2606.11650#S3.SS4)that will support the optimal recovery formulation\.

### 3\.1PoU Whitney forms and coarse spaces

We construct the learnable coarse spaces from the fine\-scaleRT0\\mathrm\{RT\}\_\{0\}/d​g​P0dgP\_\{0\}discretization\. Consider the fine\-scale piecewise\-constant scalar spaceQhfQ\_\{h\}^\{f\}and the corresponding fine\-scale lowest\-orderR​T0RT\_\{0\}flux spaceVhfV\_\{h\}^\{f\}defined as

\(24\)Qhf≔span⁡\{ϕaP0:a∈𝒞h\},Vhf≔span⁡\{ϕeR​T:e∈ℰh\},Q\_\{h\}^\{f\}\\coloneqq\\operatorname\{span\}\\\{\\phi\_\{a\}^\{P\_\{0\}\}:a\\in\\mathcal\{C\}\_\{h\}\\\},\\qquad V\_\{h\}^\{f\}\\coloneqq\\operatorname\{span\}\\\{\\phi\_\{e\}^\{RT\}:e\\in\\mathcal\{E\}\_\{h\}\\\},where𝒞h=\{1,…,Ncell\}\\mathcal\{C\}\_\{h\}=\\\{1,\\ldots,N\_\{\\mathrm\{cell\}\}\\\}denotes the set of fine cell indices andℰh\\mathcal\{E\}\_\{h\}is the set of edges for fineR​T0RT\_\{0\}flux degrees of freedom\. The basis\{ϕaP0\}\\\{\\phi\_\{a\}^\{P\_\{0\}\}\\\}is chosen so that

∑a∈𝒞hϕaP0​\(x\)=1\.\\sum\_\{a\\in\\mathcal\{C\}\_\{h\}\}\\phi\_\{a\}^\{P\_\{0\}\}\(x\)=1\.
In this work, we only use the 0\-form and 1\-form components needed for theH​\(div\)−L2H\(\\mathrm\{div\}\)\-L^\{2\}tail\. Furthermore, instead of directly constructing the Whitney forms using the classic barycentric interpolant, we employ neural networks to construct neural Whitney forms that still possess the desired PoU property\.

LetN0N\_\{0\}denote the number of partitions\. For a fixed conditioning variable, letW=\[Wi​a\]∈ℝN0×Nc​e​l​lW=\[W\_\{ia\}\]\\in\\mathbb\{R\}^\{N\_\{0\}\\times N\_\{cell\}\}be a learnable PoU weight matrix satisfying:

\(25\)Wi​a≥0,∑i=1N0Wi​a=1for every​a∈𝒞h\.W\_\{ia\}\\geq 0,\\qquad\\sum\_\{i=1\}^\{N\_\{0\}\}W\_\{ia\}=1\\quad\\text\{for every \}a\\in\\mathcal\{C\}\_\{h\}\.Then the coarse 0\-forms and 1\-forms are constructed as follows\.

###### Proposition 3\.1\(Coarse 0\-forms\)\.

For a fixed conditioning variable, letWWbe the learned PoU weights as in[Equation25](https://arxiv.org/html/2606.11650#S3.E25)\. The coarse 0\-form is defined as

\(26\)ψi0=∑a∈𝒞hWi​a​ϕaP0,i=1,…,N0,\\psi\_\{i\}^\{0\}=\\sum\_\{a\\in\\mathcal\{C\}\_\{h\}\}W\_\{ia\}\\phi\_\{a\}^\{P\_\{0\}\},\\qquad i=1,\\ldots,N\_\{0\},which forms a partition of unity:

\(27\)∑i=1N0ψi0=∑a∈𝒞h\(∑i=1N0Wi​a\)​ϕaP0=1\.\\sum\_\{i=1\}^\{N\_\{0\}\}\\psi\_\{i\}^\{0\}=\\sum\_\{a\\in\\mathcal\{C\}\_\{h\}\}\\left\(\\sum\_\{i=1\}^\{N\_\{0\}\}W\_\{ia\}\\right\)\\phi\_\{a\}^\{P\_\{0\}\}=1\.The corresponding coarseP0P\_\{0\}space is defined as:

\(28\)Qhc≔𝒲0=span​\{ψi0\}i⊂Qhf\.Q^\{c\}\_\{h\}\\coloneqq\\mathcal\{W\}^\{0\}=\\mathrm\{span\}\\\{\\psi^\{0\}\_\{i\}\\\}\_\{i\}\\subset Q\_\{h\}^\{f\}\.

Next, we construct the coarse 1\-forms that still follow the same PoU Whitney principle with a special design of separating the interior edges and boundary edges\. This separation allows us to maintain the boundary information more efficiently during the coarsening procedure, especially for scenarios where the geometry is complex or where boundary conditions are different for boundary subdomains\.

###### Proposition 3\.2\(Interior and boundary coarse 1\-forms\)\.

For a fixed conditioning variable, letWWbe the learned PoU weights as in[Equation25](https://arxiv.org/html/2606.11650#S3.E25)\. Lete∈ℰinte\\in\\mathcal\{E\}\_\{\\mathrm\{int\}\}denote the interior fine edges with two adjacent cellsKL​\(e\),KR​\(e\)K\_\{L\}\(e\),\\ K\_\{R\}\(e\)\. For each pair of coarse 0\-forms\(i,j\)\(i,j\)with1≤i<j≤N01\\leq i<j\\leq N\_\{0\}, the coarse interior 1\-form is defined as

\(29\)𝝍i​j1,int​\(x\)=∑e∈ℰint\(Wi,KL​\(e\)​Wj,KR​\(e\)−Wj,KL​\(e\)​Wi,KR​\(e\)\)​ϕe1​\(x\)\.\\boldsymbol\{\\psi\}\_\{ij\}^\{1,\\mathrm\{int\}\}\(x\)\\;=\\;\\sum\_\{e\\in\\mathcal\{E\}\_\{\\mathrm\{int\}\}\}\\Bigl\(W\_\{i,K\_\{\\mathrm\{L\}\}\(e\)\}\\,W\_\{j,K\_\{\\mathrm\{R\}\}\(e\)\}\-W\_\{j,K\_\{\\mathrm\{L\}\}\(e\)\}\\,W\_\{i,K\_\{\\mathrm\{R\}\}\(e\)\}\\Bigr\)\\,\\boldsymbol\{\\phi\}\_\{e\}^\{1\}\(x\)\.Let the boundary fine edges be partitioned into disjoint groups

ℰbc=ℰ1∪⋯∪ℰr,ℰγ∩ℰγ′=∅\(γ≠γ′\),\\mathcal\{E\}\_\{\\rm bc\}=\\mathcal\{E\}\_\{1\}\\cup\\cdots\\cup\\mathcal\{E\}\_\{r\},\\qquad\\mathcal\{E\}\_\{\\gamma\}\\cap\\mathcal\{E\}\_\{\\gamma^\{\\prime\}\}=\\emptyset\\quad\(\\gamma\\neq\\gamma^\{\\prime\}\),where each group may represent a distinct boundary component or boundary type\. Lete∈ℰγe\\in\\mathcal\{E\}\_\{\\mathrm\{\\gamma\}\}denote the boundary fine edges with one adjacent cellK​\(e\)K\(e\)\. Forγ=1,…,r\\gamma=1,\\dots,randi=1,…,N0i=1,\\dots,N\_\{0\}, the coarse boundary 1\-form is defined as

\(30\)𝝍i,γ1,bc​\(x\)=∑e∈ℰγWi,K​\(e\)​ϕe1​\(x\)\.\\boldsymbol\{\\psi\}\_\{i,\\gamma\}^\{1,\\mathrm\{bc\}\}\(x\)\\;=\\;\\sum\_\{e\\in\\mathcal\{E\}\_\{\\gamma\}\}W\_\{i,K\(e\)\}\\,\\boldsymbol\{\\phi\}\_\{e\}^\{1\}\(x\)\.Concatenating the interior and boundary 1\-forms, we have:

𝝍1=\[𝝍1,int𝝍1,bc\],\\boldsymbol\{\\psi\}^\{1\}=\[\\boldsymbol\{\\psi\}^\{1,\\mathrm\{int\}\}\\quad\\boldsymbol\{\\psi\}^\{1,\\mathrm\{bc\}\}\],whereN1=N1i​n​t\+N1b​c=\(N02\)\+r​N0N\_\{1\}=N\_\{1\}^\{int\}\+N\_\{1\}^\{bc\}=\\binom\{N\_\{0\}\}\{2\}\+rN\_\{0\}\. The corresponding coarseR​T0RT\_\{0\}space is

\(31\)Vhc≔𝒲1=span\{𝝍α1\}α=1N1⊂Vhf\.V\_\{h\}^\{c\}\\coloneqq\\mathcal\{W\}^\{1\}=\\operatorname\{span\}\\\{\\boldsymbol\{\\psi\}\_\{\\alpha\}^\{1\}\\\}\_\{\\alpha=1\}^\{N\_\{1\}\}\\subset V\_\{h\}^\{f\}\.

### 3\.2Coarse operators

For convenience, we hereafter denote by\{𝝍α1\}α=1N1\\\{\\boldsymbol\{\\psi\}\_\{\\alpha\}^\{1\}\\\}\_\{\\alpha=1\}^\{N\_\{1\}\}the full set of interior and boundary coarse 1\-form basis functions\. With the coarse 0\-forms\{ψi0\}i=1N0\\\{\\psi\_\{i\}^\{0\}\\\}\_\{i=1\}^\{N\_\{0\}\}and 1\-forms\{ψα1\}α=1N1\\\{\\psi\_\{\\alpha\}^\{1\}\\\}\_\{\\alpha=1\}^\{N\_\{1\}\}, we can further define the corresponding coarse operators\.

###### Proposition 3\.4\(Coarse operators\)\.

Fori,j=1,…,N0,i,j=1,\\dots,N\_\{0\},andα,β=1,…,N1\\alpha,\\beta=1,\\dots,N\_\{1\}, the coarse mass matrices and divergence matrix are defined by

\(𝐌0\)i​j=\(ψi0,ψj0\)Ω,\(𝐌1\)α​β=\(𝝍α1,𝝍β1\)Ω,𝐃i​α=\(∇⋅𝝍α1,ψi0\)Ω\.\\displaystyle\(\\mathbf\{M\}\_\{0\}\)\_\{ij\}=\(\\psi\_\{i\}^\{0\},\\ \\psi\_\{j\}^\{0\}\)\_\{\\Omega\},\\quad\(\\mathbf\{M\}\_\{1\}\)\_\{\\alpha\\beta\}=\(\\boldsymbol\{\\psi\}\_\{\\alpha\}^\{1\},\\ \\boldsymbol\{\\psi\}\_\{\\beta\}^\{1\}\)\_\{\\Omega\},\\quad\\mathbf\{D\}\_\{i\\alpha\}=\(\\nabla\\\!\\cdot\\\!\\boldsymbol\{\\psi\}\_\{\\alpha\}^\{1\},\\psi\_\{i\}^\{0\}\)\_\{\\Omega\}\.

The above construction for coarse bases and operators is purely algebraic and therefore remains valid for anyWWsatisfying the property in[Equation25](https://arxiv.org/html/2606.11650#S3.E25)\. For a conditioning variablezz, the neural coarsening map givesW=Wθ​\(z\)W=W\_\{\\theta\}\(z\)and the resulting conditional operators in reduced space areQhc​\(z;θ\),Vhc​\(z;θ\),Dθ​\(z\)Q\_\{h\}^\{c\}\(z;\\theta\),\\ V\_\{h\}^\{c\}\(z;\\theta\),\\ D\_\{\\theta\}\(z\)\. In practice, forNNtraining instances,WWcan be a tensor of shapeW∈ℝN×N0×NcellW\\in\\mathbb\{R\}^\{N\\times N\_\{0\}\\times N\_\{\\rm cell\}\}, whereW\(n\)∈ℝN0×NcellW^\{\(n\)\}\\in\\mathbb\{R\}^\{N\_\{0\}\\times N\_\{\\rm cell\}\}represents the conditional neural Whitney form for samplenngiven a conditioning variablez\(n\)z^\{\(n\)\}\. The choice of coarse forms may vary depending on different factors, such as the treatment of boundaries or the choices of finite element spaces, but the fundamental principle is that by introducing the trainableWW, we are able to adapt the PoU flexibly with respect to certain conditioning variables\. Moreover, the learned coarse spaces inherit theH​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}structure from the fine space with the coarse divergence operator as in the following diagram:

\(32\)H​\(div,Ω\)→∇⋅L2​\(Ω\)∪∪Vhc​\(z;θ\)→Dθ​\(z\)Qhc​\(z;θ\)\.\\begin\{array\}\[\]\{ccc\}H\(\\mathrm\{div\},\\Omega\)&\\xrightarrow\{\\;\\nabla\\\!\\cdot\\;\}&L^\{2\}\(\\Omega\)\\\\\[3\.99994pt\] \\cup&&\\cup\\\\\[\-1\.99997pt\] V\_\{h\}^\{c\}\(z;\\theta\)&\\xrightarrow\{\\;D\_\{\\theta\}\(z\)\\;\}&Q\_\{h\}^\{c\}\(z;\\theta\)\.\\end\{array\}
The following proposition provides a sufficient condition for the full row rank of the coarse divergence matrix\. This condition is a nondegeneracy condition of the learned coarse basis and can be checked from the learned PoU weights\.

###### Proposition 3\.5\(Rank condition\)\.

Let the boundary fine edges be partitioned as in[Proposition3\.2](https://arxiv.org/html/2606.11650#S3.Thmtheorem2)\. Assume that the fine boundary basisϕe1\\boldsymbol\{\\phi\}\_\{e\}^\{1\}for boundary edgee∈ℰγe\\in\\mathcal\{E\}\_\{\\gamma\}is oriented so thatβeγ≔\(∇⋅ϕe1,ϕK​\(e\)P0\)Ω\>0\\beta\_\{e\}^\{\\gamma\}\\coloneqq\(\\nabla\\\!\\cdot\\boldsymbol\{\\phi\}\_\{e\}^\{1\},\\phi\_\{K\(e\)\}^\{P\_\{0\}\}\)\_\{\\Omega\}\>0\. LetDγ∈ℝN0×N0D\_\{\\gamma\}\\in\\mathbb\{R\}^\{N\_\{0\}\\times N\_\{0\}\}denote the block of the coarse divergence matrixDDwhose columns correspond to the boundary basis functions associated withℰγ\\mathcal\{E\}\_\{\\gamma\}\. Then

\(Dγ\)i​α=∑e∈ℰγβeγ​Wα,K​\(e\)​Wi,K​\(e\)=\(Wγ​Bγ​Wγ⊤\)i​α,\(D\_\{\\gamma\}\)\_\{i\\alpha\}=\\sum\_\{e\\in\\mathcal\{E\}\_\{\\gamma\}\}\\beta\_\{e\}^\{\\gamma\}W\_\{\\alpha,K\(e\)\}W\_\{i,K\(e\)\}=\(W\_\{\\gamma\}B\_\{\\gamma\}W\_\{\\gamma\}^\{\\top\}\)\_\{i\\alpha\},where

\(Wγ\)i​e=Wi,K​\(e\),Bγ=diag⁡\{βeγ:e∈ℰγ\}\.\(W\_\{\\gamma\}\)\_\{ie\}=W\_\{i,K\(e\)\},\\qquad B\_\{\\gamma\}=\\operatorname\{diag\}\\\{\\beta\_\{e\}^\{\\gamma\}:e\\in\\mathcal\{E\}\_\{\\gamma\}\\\}\.If∑γ=1rWγ​Bγ​Wγ⊤\\sum\_\{\\gamma=1\}^\{r\}W\_\{\\gamma\}B\_\{\\gamma\}W\_\{\\gamma\}^\{\\top\}is positive definite, then the coarse divergence matrix has full row rank,rank​\(D\)=N0\\text\{rank\}\(D\)=N\_\{0\}\.

###### Proof 3\.6\.

See Appendix[A](https://arxiv.org/html/2606.11650#A1)\.

### 3\.3Surjectivity property underH​\(div\)H\(\\mathrm\{div\}\)setting and flux conservation

The divergence operator∇⋅\\nabla\\cdotis surjective from the Raviart–Thomas spaceRTk​\(K\)\\mathrm\{RT\}\_\{k\}\(K\)ontoPk​\(K\)P\_\{k\}\(K\)in the standard fine scale\[boffi2013mixed\]\. In particular, forRT0\\mathrm\{RT\}\_\{0\}/d​g​P0dgP\_\{0\}, the divergence represents the cellwise flux conservation\. Here, we extend this surjectivity property and conservation structure to coarse operators induced by the Whitney forms construction\.

###### Proposition 3\.7\.

Let𝒲1≔span​\{ψa1\}a=1N1,𝒲0≔span​\{ψi0\}i=1N0\\mathcal\{W\}^\{1\}\\coloneqq\\mathrm\{span\}\\\{\\psi^\{1\}\_\{a\}\\\}\_\{a=1\}^\{N\_\{1\}\},\\ \\mathcal\{W\}^\{0\}\\coloneqq\\mathrm\{span\}\\\{\\psi^\{0\}\_\{i\}\\\}\_\{i=1\}^\{N\_\{0\}\}be the coarse flux and scalar spaces following[Proposition3\.1](https://arxiv.org/html/2606.11650#S3.Thmtheorem1)and[Proposition3\.2](https://arxiv.org/html/2606.11650#S3.Thmtheorem2)\. Then:

\(33\)∇⋅𝒲1⊆𝒲0\.\\nabla\\\!\\cdot\\mathcal\{W\}^\{1\}\\subseteq\\mathcal\{W\}^\{0\}\.Moreover, if the rank condition in[Proposition3\.5](https://arxiv.org/html/2606.11650#S3.Thmtheorem5)holds, then the coarse divergence matrix has full row rank and the coarse divergence is surjective from𝒲1\\mathcal\{W\}^\{1\}onto𝒲0\\mathcal\{W\}^\{0\}\.

###### Proof 3\.8\.

Recall that we have fine\-scale basis functions\{ϕ1\}⊂R​T0\\\{\\boldsymbol\{\\phi\}^\{1\}\\\}\\subset RT\_\{0\}and\{ϕ0\}⊂P0\\\{\\phi^\{0\}\\\}\\subset P\_\{0\}\. Given the fine\-scale divergence identity for\(RT0,dgP0\)\(\\mathrm\{RT\}\_\{0\},\\mathrm\{dgP\}\_\{0\}\), we have

\(34\)∇⋅ϕeR​T0=ϕKR​\(e\)P0−ϕKL​\(e\)P0\.\\nabla\\\!\\cdot\\boldsymbol\{\\phi\}\_\{e\}^\{RT\_\{0\}\}=\\phi\_\{K\_\{R\}\(e\)\}^\{P\_\{0\}\}\-\\phi\_\{K\_\{L\}\(e\)\}^\{P\_\{0\}\}\.For interior edgee=\(i,j\)e=\(i,j\), we can derive the identity in the coarse scale by substituting the definition of coarse basis functions:

∇⋅ψi​j1\\displaystyle\\nabla\\\!\\cdot\\psi\_\{ij\}^\{1\}=∇⋅\(∑a,bWi​a​Wj​b​ϕa​bR​T0\)\\displaystyle=\\nabla\\\!\\cdot\\left\(\\sum\_\{a,b\}W\_\{ia\}W\_\{jb\}\\phi\_\{ab\}^\{RT\_\{0\}\}\\right\)=∑a,bWi​a​Wj​b​ϕbP0−∑a,bWi​a​Wj​b​ϕaP0\\displaystyle=\\sum\_\{a,b\}W\_\{ia\}W\_\{jb\}\\phi\_\{b\}^\{P\_\{0\}\}\-\\sum\_\{a,b\}W\_\{ia\}W\_\{jb\}\\phi\_\{a\}^\{P\_\{0\}\}=\(∑aWi​a\)​ψj0−\(∑bWj​b\)​ψi0\\displaystyle=\\left\(\\sum\_\{a\}W\_\{ia\}\\right\)\\psi\_\{j\}^\{0\}\-\\left\(\\sum\_\{b\}W\_\{jb\}\\right\)\\psi\_\{i\}^\{0\}=mi​ψj0−mj​ψi0∈𝒲0\.\\displaystyle=m\_\{i\}\\psi\_\{j\}^\{0\}\-m\_\{j\}\\psi\_\{i\}^\{0\}\\in\\mathcal\{W\}^\{0\}\.Similarly, by linearity, the boundary coarse 1\-form is also a linear combination of fineP0P\_\{0\}basis functions and belongs to the span of the coarse scalar space\. Therefore,

∇⋅𝒲1⊆𝒲0\.\\nabla\\\!\\cdot\\mathcal\{W\}^\{1\}\\subseteq\\mathcal\{W\}^\{0\}\.Next we prove the surjectivity of the coarse divergence operator\. Let

\(35\)𝒗=∑a=1N1F^a​ψa1∈𝒲1,q=∑i=1N0q^i​ψi0∈𝒲0\.\\boldsymbol\{v\}=\\sum\_\{a=1\}^\{N\_\{1\}\}\\widehat\{F\}\_\{a\}\\psi\_\{a\}^\{1\}\\in\\mathcal\{W\}^\{1\},\\qquad q=\\sum\_\{i=1\}^\{N\_\{0\}\}\\widehat\{q\}\_\{i\}\\psi\_\{i\}^\{0\}\\in\\mathcal\{W\}^\{0\}\.By linearity,

\(36\)∇⋅𝒗=∑a=1N1F^a​∇⋅ψa1∈𝒲0\.\\nabla\\cdot\\boldsymbol\{v\}=\\sum\_\{a=1\}^\{N\_\{1\}\}\\widehat\{F\}\_\{a\}\\nabla\\cdot\\psi\_\{a\}^\{1\}\\in\\mathcal\{W\}^\{0\}\.Then testing againstψi0\\psi\_\{i\}^\{0\}fori=1,…,N0i=1,\\dots,N\_\{0\}and expanding for coarse variables gives:

\(37\)\(∇⋅v,ψi0\)Ω\\displaystyle\(\\nabla\\cdot v,\\psi\_\{i\}^\{0\}\)\_\{\\Omega\}=\(∇⋅\(∑a=1N1F^a​𝝍a1\),ψi0\)Ω\\displaystyle=\\left\(\\nabla\\cdot\\left\(\\sum\_\{a=1\}^\{N\_\{1\}\}\\widehat\{F\}\_\{a\}\\boldsymbol\{\\psi\}\_\{a\}^\{1\}\\right\),\\psi\_\{i\}^\{0\}\\right\)\_\{\\Omega\}\(38\)=∑a=1N1F^a​\(∇⋅𝝍a1,ψi0\)Ω=∑a=1N1F^a​Di​a\\displaystyle=\\sum\_\{a=1\}^\{N\_\{1\}\}\\widehat\{F\}\_\{a\}\\left\(\\nabla\\cdot\\boldsymbol\{\\psi\}\_\{a\}^\{1\},\\psi\_\{i\}^\{0\}\\right\)\_\{\\Omega\}=\\sum\_\{a=1\}^\{N\_\{1\}\}\\widehat\{F\}\_\{a\}D\_\{ia\}\(39\)\(q,ψk0\)Ω\\displaystyle\(q,\\psi\_\{k\}^\{0\}\)\_\{\\Omega\}=∑i=1N0qi​\(ψi0,ψk0\)Ω=∑i=1N0qi​\(M0\)i​k,\\displaystyle=\\sum\_\{i=1\}^\{N\_\{0\}\}q\_\{i\}\(\\psi\_\{i\}^\{0\},\\psi\_\{k\}^\{0\}\)\_\{\\Omega\}=\\sum\_\{i=1\}^\{N\_\{0\}\}q\_\{i\}\(M\_\{0\}\)\_\{ik\},whereDDandMMare as defined in[Proposition3\.4](https://arxiv.org/html/2606.11650#S3.Thmtheorem4)\. Whenrank⁡\(D\)=N0\\operatorname\{rank\}\(D\)=N\_\{0\}, the mapD:ℝN1→ℝN0D:\\mathbb\{R\}^\{N\_\{1\}\}\\to\\mathbb\{R\}^\{N\_\{0\}\}is onto and thus for arbitraryqq, there exist coefficientsF^∈ℝN1\\widehat\{F\}\\in\\mathbb\{R\}^\{N\_\{1\}\}such thatD​F^=M0​qD\\widehat\{F\}=M\_\{0\}q\. For this choice ofF^\\widehat\{F\}, the corresponding𝐯∈𝒲1\\boldsymbol\{v\}\\in\\mathcal\{W\}^\{1\}satisfies

\(40\)\(∇⋅𝒗,ψk0\)Ω=\(q,ψk0\)Ω\.\(\\nabla\\\!\\cdot\\boldsymbol\{v\},\\psi\_\{k\}^\{0\}\)\_\{\\Omega\}=\(q,\\psi\_\{k\}^\{0\}\)\_\{\\Omega\}\.Thus we conclude that the coarse divergence operator is surjective\.

Moreover, since the coefficient in \([29](https://arxiv.org/html/2606.11650#S3.E29)\) is the discrete analogue of the classical Whitney exchange formλi​∇λj−λj​∇λi\\lambda\_\{i\}\\nabla\\lambda\_\{j\}\-\\lambda\_\{j\}\\nabla\\lambda\_\{i\}, the resulting antisymmetry property∇⋅ψi​j1=−∇⋅ψj​i1\\nabla\\\!\\cdot\\psi^\{1\}\_\{ij\}=\-\\,\\nabla\\\!\\cdot\\psi^\{1\}\_\{ji\}with respect to the exchange of indices ensures by construction that the 1\-forms encode equal and opposite fluxes exchanged between 0\-forms\. Summing over all indices cancels the internal fluxes, so that only the boundary contributions remain\. This recovers a discrete conservation law consistent with the divergence theorem and shows that the learned coarse construction preserves the flux balance\.

### 3\.4Graph calculus and circuit interpretation for FEEC

The preceding construction in[Section3\.1](https://arxiv.org/html/2606.11650#S3.SS1)and[Section3\.2](https://arxiv.org/html/2606.11650#S3.SS2)produces anH​\(div\)H\(\\mathrm\{div\}\)\-conforming reduced space where conservation is preserved by the coarse divergence operator\. To formulate the GP optimal recovery problem in the next section, we now reinterpret the reduced subcomplexes as a graph\-valued circuit model and connect coarse operators with the graph calculus\. This interpretation bridges the structure\-preserving FEEC with the graph\-based GP model used to learn the nonlinearity and perform uncertainty quantification later\. The viewpoint is related to the computational graph completion \(CGC\)\[owhadi2022computational\]problem, where the unknown functions and unobserved variables can be recovered on a given graph\. In contrast, the graph in this work is induced by the coarse finite element subcomplexes, which can be considered an autonomous identification of a dense graph\.

Let𝒢=\(𝒱,ℰ\)\\mathcal\{G\}=\(\\mathcal\{V\},\\mathcal\{E\}\)be a graph with\|𝒱\|=N0\\lvert\\mathcal\{V\}\\rvert=N\_\{0\}vertices and\|ℰ\|=N1\\lvert\\mathcal\{E\}\\rvert=N\_\{1\}edges\. The coarse spaces constructed via the partition of unity Whitney forms naturally define a graph where each vertexi∈𝒱i\\in\\mathcal\{V\}corresponds to a coarse 0\-formψi0\\psi^\{0\}\_\{i\}and each edgee∈ℰe\\in\\mathcal\{E\}corresponds to a coarse11\-formψe1\\psi^\{1\}\_\{e\}\. We define the cochain spaces as:

C0​\(𝒢\)≔ℝ\|V\|,C1​\(𝒢\)≔ℝ\|ℰ\|\.C^\{0\}\(\\mathcal\{G\}\)\\coloneqq\\mathbb\{R\}^\{\\lvert V\\rvert\},\\qquad C^\{1\}\(\\mathcal\{G\}\)\\coloneqq\\mathbb\{R\}^\{\\lvert\\mathcal\{E\}\\rvert\}\.Then the coarse variables admit the interpretation on vertices and edges, respectively:

𝐮^∈C0​\(𝒢\),𝐅^∈C1​\(𝒢\),\\widehat\{\\mathbf\{u\}\}\\in C^\{0\}\(\\mathcal\{G\}\),\\qquad\\widehat\{\\mathbf\{F\}\}\\in C^\{1\}\(\\mathcal\{G\}\),where𝐮^\\widehat\{\\mathbf\{u\}\}represents node states and𝐅^\\widehat\{\\mathbf\{F\}\}represents oriented fluxes or currents\.

In order to formulate the mapping between nodal values and edge values, we introduce the edge\-vertex signed incidence matrixB0∈ℝE×VB\_\{0\}\\in\\mathbb\{R\}^\{E\\times V\}defined by

\(41\)\(B0\)i,j=\{−1if edge​ei​leaves vertex​vj,1if edge​ei​enters vertex​vj,0otherwise\.\(B\_\{0\}\)\_\{i,j\}=\\begin\{cases\}\-1&\\text\{if edge \}e\_\{i\}\\text\{ leaves vertex \}v\_\{j\},\\\\ 1&\\text\{if edge \}e\_\{i\}\\text\{ enters vertex \}v\_\{j\},\\\\ 0&\\text\{otherwise\}\.\\end\{cases\}Subsequently for an oriented edgee=\(i,j\)e=\(i,j\), we can define the discrete graph gradient operator and graph divergence as:

\(Graph gradient\)\(δ0​𝐮^\)​\(e\)=u^j−u^i,\\displaystyle\\quad\(\\delta\_\{0\}\\widehat\{\\mathbf\{u\}\}\)\(e\)=\\widehat\{u\}\_\{j\}\-\\widehat\{u\}\_\{i\},\(Graph divergence\)\(δ0⊤​𝐅^\)​\(v\)=∑vi∼vF^​\(vi,v\)−∑v∼vjF^​\(v,vj\)\.\\displaystyle\\quad\(\\delta\_\{0\}^\{\\top\}\\widehat\{\\mathbf\{F\}\}\)\(v\)=\\sum\_\{v\_\{i\}\\sim v\}\\widehat\{F\}\(v\_\{i\},v\)\-\\sum\_\{v\\sim v\_\{j\}\}\\widehat\{F\}\(v,v\_\{j\}\)\.The discrete graph gradient, mapping the nodal values to edge values, captures the difference across edges\. Conversely, the graph divergence, mapping edge values back to values at nodes, measures the net flux entering or leaving the nodevvover its adjacent edges\. LetD𝒢∈ℝ\|𝒱\|×\|ℰ\|D\_\{\\mathcal\{G\}\}\\in\\mathbb\{R\}^\{\|\\mathcal\{V\}\|\\times\|\\mathcal\{E\}\|\}denote the graph divergence matrix, defined by

\(D𝒢​𝐅^\)i=∑e∈ℰDi​e​F^e\.\(D\_\{\\mathcal\{G\}\}\\widehat\{\\mathbf\{F\}\}\)\_\{i\}=\\sum\_\{e\\in\\mathcal\{E\}\}D\_\{ie\}\\widehat\{F\}\_\{e\}\.With the construction of coarse Whitney forms in[Section3\.2](https://arxiv.org/html/2606.11650#S3.SS2),DDhas the same incidence structure as the graph divergence operator, possibly weighted by the mass or basis normalization depending on the specific construction\. Thus, this gives the conservation law in the reduced spaces and supplies the finite\-dimensional variables on which the GP acts\.

## 4Method part 2: Optimal recovery of physics with quantified uncertainty

In this section, we develop a probabilistic representation of physical fields on the graph induced by the coarse spaces\. This provides a unified framework where the physical law acts on discrete structures while uncertainty is quantified in a functionally consistent manner\.

### 4\.1Structure\-preserving Optimal Recovery Problem

We now introduce a GP optimal recovery formulation for learning the nonlinear flux correction in the coarseH​\(div\)−L2H\(\\mathrm\{div\}\)\-L^\{2\}space, consistent with the FEEC structure established above\. For each coarse edgee=\(i,j\)∈ℰe=\(i,j\)\\in\\mathcal\{E\}, we can define the local input with different choices:

𝐮^e≔\{\(u^i,u^j\)∈ℝ2,\(endpoint representation\),δ0​u^e=u^j−u^i∈ℝ,\(edge gradient representation\)\.\\hat\{\\mathbf\{u\}\}\_\{e\}\\coloneqq\\begin\{cases\}\(\\hat\{u\}\_\{i\},\\hat\{u\}\_\{j\}\)\\in\\mathbb\{R\}^\{2\},&\\text\{\(endpoint representation\)\},\\\\\[4\.0pt\] \\delta\_\{0\}\\hat\{u\}\_\{e\}=\\hat\{u\}\_\{j\}\-\\hat\{u\}\_\{i\}\\in\\mathbb\{R\},&\\text\{\(edge gradient representation\)\}\.\\end\{cases\}For training instancesn=1,…,Nn=1,\\ldots,N, the desired GP input features become

xe\(n\)=\(𝐮^e\(n\),z\(n\),s\(n\)\)\.x\_\{e\}^\{\(n\)\}=\(\\widehat\{\\mathbf\{u\}\}\_\{e\}^\{\(n\)\},z^\{\(n\)\},s^\{\(n\)\}\)\.Here𝐮^e\(n\)\\widehat\{\\mathbf\{u\}\}\_\{e\}^\{\(n\)\}is the coarse state variable associated with edgeee,z\(n\)z^\{\(n\)\}denotes the conditioning variable, ands\(n\)s^\{\(n\)\}denotes the geometric embedding\.

In the independent kernel construction, for a fixed edgeee, we collect the training inputs asXe=\{xe\(n\)\}n=1NX\_\{e\}=\\\{x\_\{e\}^\{\(n\)\}\\\}\_\{n=1\}^\{N\}and the corresponding global kernel is the block\-diagonal matrix defined as

\(42\)Ke≔Ke​\(Xe,Xe\)∈ℝN×N,K=diag​\(K1,…,KN1\)∈ℝN1​N×N1​N\.K\_\{e\}\\coloneqq K\_\{e\}\(X\_\{e\},X\_\{e\}\)\\in\\mathbb\{R\}^\{N\\times N\},\\qquad K=\\text\{diag\}\\left\(K\_\{1\},\\ldots,K\_\{N\_\{1\}\}\\right\)\\in\\mathbb\{R\}^\{N\_\{1\}N\\times N\_\{1\}N\}\.The GP model on edgeeemaps local features to the flux on the edge:

\(43\)fe​\(Xe\)↦F^egp,𝐅^gp=vec​\(\[F^egp,\(n\)\]\)∈ℝN1​N,f\_\{e\}\(X\_\{e\}\)\\;\\mapsto\\;\\widehat\{F\}^\{\\rm gp\}\_\{e\},\\qquad\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}=\\text\{vec\}\\left\(\\left\[\\widehat\{F\}\_\{e\}^\{\\rm gp,\(n\)\}\\right\]\\right\)\\in\\mathbb\{R\}^\{N\_\{1\}N\},where eachfef\_\{e\}belongs to a reproducing kernel Hilbert spaceℋKe\\mathcal\{H\}\_\{K\_\{e\}\}and𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}is the vectorized global output over all edges\.

Instead of assigning independent GPs to each edge, we may also stack all edge samples into a global matrix

\(44\)X=\{xe\(n\):e=1,…,N1,n=1,…,N\}X=\\\{x\_\{e\}^\{\(n\)\}:e=1,\\ldots,N\_\{1\},\\ n=1,\\ldots,N\\\}and construct a single shared GP across all cochains:

K=\[kθ​\(xe\(n\),xe′\(m\)\)\]\(e,n\),\(e′,m\)∈ℝN1​N×N1​NK=\\left\[k\_\{\\theta\}\\\!\\left\(x\_\{e\}^\{\(n\)\},x\_\{e^\{\\prime\}\}^\{\(m\)\}\\right\)\\right\]\_\{\(e,n\),\(e^\{\\prime\},m\)\}\\in\\mathbb\{R\}^\{N\_\{1\}N\\times N\_\{1\}N\}wherekθk\_\{\\theta\}denotes the covariance kernel\. Since both cases use a large global matrix, for simplicity we denote the global matrix asKKin the following\.

Recall the governing equation for the physical law in[Equation1](https://arxiv.org/html/2606.11650#S1.E1)\. By restricting the PDE to data\-driven Whitney spaces𝒲0​\(Ω;θ\)\\mathcal\{W\}^\{0\}\(\\Omega;\\theta\)and𝒲1​\(Ω;θ\)\\mathcal\{W\}^\{1\}\(\\Omega;\\theta\), the goal is to find\(u,𝐅\)∈𝒲0​\(Ω;θ\)×𝒲1​\(Ω;θ\)\(u,\\mathbf\{F\}\)\\in\\mathcal\{W\}^\{0\}\(\\Omega;\\theta\)\\times\\mathcal\{W\}^\{1\}\(\\Omega;\\theta\)such that

\(45\)\(𝐅,𝐯\)Ω\\displaystyle\(\\mathbf\{F\},\\mathbf\{v\}\)\_\{\\Omega\}=\(u,∇⋅𝐯\)Ω\+\(𝒩​\[u\],𝐯\)Ω−⟨uD,𝐯⋅𝐧⟩ΓD,\\displaystyle=\(u,\\nabla\\cdot\\mathbf\{v\}\)\_\{\\Omega\}\+\(\\mathcal\{N\}\[u\],\\mathbf\{v\}\)\_\{\\Omega\}\-\\langle u\_\{D\},\\mathbf\{v\}\\cdot\\mathbf\{n\}\\rangle\_\{\\Gamma\_\{D\}\},𝐯∈𝒲1​\(Ω;θ\)\\displaystyle\\qquad\\mathbf\{v\}\\in\\mathcal\{W\}^\{1\}\(\\Omega;\\theta\)\(46\)\(∇⋅𝐅,q\)Ω\\displaystyle\(\\nabla\\cdot\\mathbf\{F\},q\)\_\{\\Omega\}=\(f,q\)Ω,\\displaystyle=\(f,q\)\_\{\\Omega\},q∈𝒲0​\(Ω;θ\),\\displaystyle\\qquad q\\in\\mathcal\{W\}^\{0\}\(\\Omega;\\theta\),whereffis the source term and⟨uD,𝐯⋅𝐧⟩ΓD=∫ΓDuD​\(𝐯⋅𝐧\)​𝑑s\\langle u\_\{D\},\\mathbf\{v\}\\cdot\\mathbf\{n\}\\rangle\_\{\\Gamma\_\{D\}\}\\;=\\;\\int\_\{\\Gamma\_\{D\}\}u\_\{D\}\\,\(\\mathbf\{v\}\\cdot\\mathbf\{n\}\)\\,dsis the weak Dirichlet boundary condition\. Recall that𝐮^\\widehat\{\\mathbf\{u\}\}denotes the coarse scalar coefficient\. Let𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}denote the nonlinear flux correction on the coarse edge and𝐅^\\widehat\{\\mathbf\{F\}\}denote the total coarse flux coefficient\. Using the FEEC structure, the problem can be further written as

\(47\)𝐌1​𝐅^\\displaystyle\\mathbf\{M\}\_\{1\}\\widehat\{\\mathbf\{F\}\}=D⊤​𝐮^\+𝐌1​𝐅^gp−𝐠^D,\\displaystyle=D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\+\\mathbf\{M\}\_\{1\}\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\-\\widehat\{\\mathbf\{g\}\}\_\{D\},\(48\)D​𝐅^\\displaystyle D\\widehat\{\\mathbf\{F\}\}=𝐟^,\\displaystyle=\\widehat\{\\mathbf\{f\}\},where𝐟^\\widehat\{\\mathbf\{f\}\}is the coefficient vector of the source term and𝐠^D\\widehat\{\\mathbf\{g\}\}\_\{D\}represents the Dirichlet boundary contribution with𝐠^D=∫ΓDuD​𝝍1⋅𝒏​𝑑s\\widehat\{\\mathbf\{g\}\}\_\{D\}=\\int\_\{\\Gamma\_\{D\}\}u\_\{D\}\\boldsymbol\{\\psi\}^\{1\}\\cdot\\boldsymbol\{n\}\\,ds\.𝐌1\\mathbf\{M\}\_\{1\}andDDare the discrete mass matrix and divergence operator as defined in[Proposition3\.4](https://arxiv.org/html/2606.11650#S3.Thmtheorem4)\. For simplicity, we letθ=\(θG​P,θW\)\\theta=\(\\theta\_\{GP\},\\theta\_\{W\}\)denote the set of trainable parameters for the GP model and neural Whitney forms\. Eliminating𝐅^\\widehat\{\\mathbf\{F\}\}gives the reduced residual:

\(49\)Rθ​\(𝐮^,𝐅^gp;z\):=D​𝐌1−1​D⊤​𝐮^\+D​𝐅^gp−𝐛^,with​𝐛^=𝐟^\+D​𝐌1−1​𝐠^D\.R\_\{\\theta\}\(\\widehat\{\\mathbf\{u\}\},\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\};z\):=D\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\+D\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\-\\widehat\{\\mathbf\{b\}\},\\qquad\\text\{with \}\\widehat\{\\mathbf\{b\}\}=\\widehat\{\\mathbf\{f\}\}\+D\\mathbf\{M\}\_\{1\}^\{\-1\}\\widehat\{\\mathbf\{g\}\}\_\{D\}\.
Finally, we arrive at the full optimization problem, combining the GP recovery and the data\-driven Whitney forms on the graph:

\(50\)minθ,𝐮^⁡min𝐅^gp\(𝐅^gp\)⊤​\(K​\(X,X\)\+σε2​I\)−1​𝐅^gp\+log​det\(K​\(X,X\)\+σε2​I\)\+∑n=1N‖𝐮data\(n\)−𝐮approx\(n\)‖2\+∑n=1N\|𝐅Γ,data\(n\)−𝐅Γ,approx\(n\)\|2s\.t\.ℛθ​\(𝐮^,𝐅^gp;z\(n\)\)=0,\\begin\{split\}\\min\_\{\\theta,\\widehat\{\\mathbf\{u\}\}\}\\min\_\{\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\}&\\quad\(\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\)^\{\\top\}\(K\(X,X\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\)^\{\-1\}\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\+\\log\\det\(K\(X,X\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\)\\\\ &\\quad\+\\sum\_\{n=1\}^\{N\}\\left\\\|\\mathbf\{u\}\_\{\\rm data\}^\{\(n\)\}\-\\mathbf\{u\}\_\{\\rm approx\}^\{\(n\)\}\\right\\\|^\{2\}\+\\sum\_\{n=1\}^\{N\}\\left\|\\mathbf\{F\}\_\{\\Gamma,\{\\rm data\}\}^\{\(n\)\}\-\\mathbf\{F\}\_\{\\Gamma,\\rm approx\}^\{\(n\)\}\\right\|^\{2\}\\\\ \\text\{s\.t\.\}\\qquad&\\mathcal\{R\}\_\{\\theta\}\\left\(\\widehat\{\\mathbf\{u\}\},\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\};z^\{\(n\)\}\\right\)=0,\\\\ \\end\{split\}where𝐮approx,𝐅approx\\mathbf\{u\}\_\{\\text\{approx\}\},\\ \\mathbf\{F\}\_\{\\text\{approx\}\}are the reconstruction of the solution field and flux back to the original fine\-scale space\. The first two terms promote smoothness in the RKHS norm and penalize overly complex models to avoid overfitting\. Moreover, as the recovery problem is now defined for the coarse coefficientsu^\\hat\{u\}and𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}, there is no guarantee that the solution and flux projected back to the original space can reconstruct the target solution𝐮data\\mathbf\{u\}\_\{\\text\{data\}\}and𝐅data\\mathbf\{F\}\_\{\\text\{data\}\}as desired\. Thus, we introduce two additional terms consisting of the mean squared error \(MSE\) reconstruction loss in the fine scale\.

The reconstruction of the solution in \([50](https://arxiv.org/html/2606.11650#S4.E50)\) can be directly obtained by mapping the coarse0\-cochain coefficients back to the fine scale with the neural Whitney bases:

\(51\)𝐮approx=𝐮^⊤​Ψ0\.\\mathbf\{u\}\_\{\\text\{approx\}\}=\\widehat\{\\mathbf\{u\}\}^\{\\top\}\\Psi\_\{0\}\.The total ground\-truth flux and reconstructed flux on the target boundaryΓ\\Gammafrom the PoU Whitney form and coarse coefficient are:

\(52\)𝐅Γ,data=∑e∈ΓFe,𝐅Γ,approx=∑i=1N​p​o​u∑e∈ΓWi​e​F^i,\\mathbf\{F\}\_\{\\Gamma,\\text\{data\}\}=\\sum\_\{e\\in\\Gamma\}F\_\{e\},\\qquad\\mathbf\{F\}\_\{\\Gamma,\\text\{approx\}\}=\\sum\_\{i=1\}^\{Npou\}\\sum\_\{e\\in\\Gamma\}W\_\{ie\}\\widehat\{F\}\_\{i\},whereWi​eW\_\{ie\}denotes the corresponding weight for theii\-th PoU and fine\-scale edgeeeand𝐅^=𝐌1−1​D⊤​u^\+𝐅^gp−𝐌1−1​𝐠D\\widehat\{\\mathbf\{F\}\}=\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\hat\{u\}\+\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\-\\mathbf\{M\}\_\{1\}^\{\-1\}\\mathbf\{g\}\_\{D\}\.

### 4\.2Bilevel optimization

The optimization problem is nested and strongly coupled because the GP inputs depend on the unknown coarse potentials𝐮^e\\widehat\{\\mathbf\{u\}\}\_\{e\}\. Therefore, we apply a bilevel optimization procedure to train the transformer for data\-driven Whitney forms and the GP model for flux prediction simultaneously\. The bilevel training algorithm is summarized in[Algorithm1](https://arxiv.org/html/2606.11650#alg1)\.

Algorithm 1Bilevel training for structure\-preserving optimal recovery0:

𝒟=\{\(z\(n\),udata\(n\),Fdata\(n\)\)\}n=1N\\mathcal\{D\}=\\\{\(z^\{\(n\)\},u^\{\(n\)\}\_\{\\text\{data\}\},F^\{\(n\)\}\_\{\\text\{data\}\}\)\\\}\_\{n=1\}^\{N\}; fine operators; PoU network

WθWW\_\{\\theta\_\{W\}\}; GP kernel

KθGPK\_\{\\theta\_\{\\mathrm\{GP\}\}\}
0:Optimizers: Adam for

θGP\\theta\_\{\\mathrm\{GP\}\}, Shampoo for

θW\\theta\_\{W\}
1:Initialize

𝐮^\\widehat\{\\mathbf\{u\}\},

θW\\theta\_\{W\}, and

θGP\\theta\_\{\\mathrm\{GP\}\}
2:forepoch

=1,…,T=1,\\ldots,Tdo

3:PoU forward:Compute

W←WθW​\(z\)W\\leftarrow W\_\{\\theta\_\{W\}\}\(z\)
4:Compute coarse bases

Ψ0,Ψ1\\Psi^\{0\},\\Psi^\{1\}, and operators

M1,D,RΓ,𝐛^M\_\{1\},D,R\_\{\\Gamma\},\\widehat\{\\mathbf\{b\}\}
5:Form coarse GP features

X=\{xe\(n\)\}e,n,K~=KθGP​\(X,X\)\+σε2​IX=\\\{x\_\{e\}^\{\(n\)\}\\\}\_\{e,n\},\\ \\widetilde\{K\}=K\_\{\\theta\_\{\\mathrm\{GP\}\}\}\(X,X\)\+\\sigma\_\{\\varepsilon\}^\{2\}I
6:Inner problem: Fix

𝐮^\\widehat\{\\mathbf\{u\}\},

θW\\theta\_\{W\},

θGP\\theta\_\{\\mathrm\{GP\}\}, solve the constrained recovery problem for

𝐅^gp,⋆\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\}
7:Update

𝐅^gp,⋆←SolveKKT⁡\(𝐊~,D,RΓ,𝐛^,FΓ,data\)\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\}\\leftarrow\\operatorname\{SolveKKT\}\\left\(\\widetilde\{\\mathbf\{K\}\},D,R\_\{\\Gamma\},\\widehat\{\\mathbf\{b\}\},F\_\{\\Gamma,\{\\rm data\}\}\\right\)
8:Update

θGP←Adam⁡\(θGP,∇θGPℒinner\)\\theta\_\{\\mathrm\{GP\}\}\\leftarrow\\operatorname\{Adam\}\\left\(\\theta\_\{\\mathrm\{GP\}\},\\nabla\_\{\\theta\_\{\\mathrm\{GP\}\}\}\\mathcal\{L\}\_\{\\mathrm\{inner\}\}\\right\)
9:Outer problem: Fix

𝐅^gp,⋆\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\},

θGP\\theta\_\{\\mathrm\{GP\}\}, solve for

𝐮^⋆\\widehat\{\\mathbf\{u\}\}^\{\\star\}with least squares method

10:Update

𝐮^⋆←lstsq⁡\(D​M1−1​D⊤,𝐛^−D​𝐅^gp,⋆\)\\widehat\{\\mathbf\{u\}\}^\{\\star\}\\leftarrow\\operatorname\{lstsq\}\\left\(DM\_\{1\}^\{\-1\}D^\{\\top\},\\,\\widehat\{\\mathbf\{b\}\}\-D\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\}\\right\)
11:Update

θW←Shampoo⁡\(θW,∇θWℒouter\)\\theta\_\{W\}\\leftarrow\\operatorname\{Shampoo\}\\left\(\\theta\_\{W\},\\nabla\_\{\\theta\_\{W\}\}\\mathcal\{L\}\_\{\\mathrm\{outer\}\}\\right\)
12:endfor

13:

14:returntrained

θW⋆,θGP⋆,𝐮^⋆\\theta\_\{W\}^\{\\star\},\\theta\_\{\\mathrm\{GP\}\}^\{\\star\},\\widehat\{\\mathbf\{u\}\}^\{\\star\}, and

𝐅^⋆\\widehat\{\\mathbf\{F\}\}^\{\\star\}\.

For the inner loop, we first fix𝐮^\\widehat\{\\mathbf\{u\}\}andθW\\theta\_\{W\}and solve for the closed\-form solution𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}via the Karush–Kuhn–Tucker \(KKT\) system\. Since the flux observations are defined on fine boundary edges, for simplicity, we introduce the coarse\-to\-fine reconstruction operatorRΓ​\(W\):ℝN1→ℝR\_\{\\Gamma\}\(W\):\\mathbb\{R\}^\{N\_\{1\}\}\\to\\mathbb\{R\}that maps the coarse flux degrees of freedom \(𝐅^\\hat\{\\mathbf\{F\}\}\) back to the total flux through the boundary regionΓ\\Gamma\. Substituting𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}gives the inner optimization problem as

\(53\)min𝐅^gp\\displaystyle\\min\_\{\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\}\\;‖𝐅Γ,data−RΓ​\(W\)​\(𝐌1−1​D⊤​𝐮^\+𝐅^gp−𝐌1−1​𝐠D\)‖22\\displaystyle\\left\\\|\\mathbf\{F\}\_\{\\Gamma,\\text\{data\}\}\-R\_\{\\Gamma\}\(W\)\\left\(\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\+\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\-\\mathbf\{M\}\_\{1\}^\{\-1\}\\mathbf\{g\}\_\{D\}\\right\)\\right\\\|\_\{2\}^\{2\}\(54\)\+\(𝐅^gp\)⊤​\(K\+σε2​I\)−1​𝐅^gp\+log​det\(K\+σ2​I\)\\displaystyle\\quad\+\(\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\)^\{\\top\}\(K\+\\sigma\_\{\\varepsilon\}^\{2\}I\)^\{\-1\}\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\+\\log\\det\(K\+\\sigma^\{2\}I\)\(55\)s\.t\.D​𝐅^gp=𝐛−D​𝐌1−1​D⊤​𝐮^\.\\displaystyle D\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}=\\mathbf\{b\}\-D\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\.Since the termlog​det\(K\+σε2​I\)\\log\\det\(K\+\\sigma\_\{\\varepsilon\}^\{2\}I\)does not depend on𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}, it can be omitted when solving for𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\. Introducing the Lagrange multiplierλ\\lambda, the optimality system becomes:

\[K~−1\+RΓ⊤​RΓD⊤D0\]​\[𝐅^gpλ\]=\[RΓ⊤​\(𝐅Γ,data−RΓ​\(𝐌1−1​D⊤​𝐮^−𝐌1−1​𝐠D\)\)𝐛−D​𝐌1−1​D⊤​𝐮^\]\.\\begin\{bmatrix\}\\tilde\{K\}^\{\-1\}\+R\_\{\\Gamma\}^\{\\top\}R\_\{\\Gamma\}&D^\{\\top\}\\\\ D&0\\end\{bmatrix\}\\begin\{bmatrix\}\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}\\\\ \\lambda\\end\{bmatrix\}=\\begin\{bmatrix\}R\_\{\\Gamma\}^\{\\top\}\\Bigl\(\\mathbf\{F\}\_\{\\Gamma,\\text\{data\}\}\-R\_\{\\Gamma\}\\left\(\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\-\\mathbf\{M\}\_\{1\}^\{\-1\}\\mathbf\{g\}\_\{D\}\\right\)\\Bigr\)\\\\ \\mathbf\{b\}\-D\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\\end\{bmatrix\}\.whereK~=K\+σ2​I\\tilde\{K\}=K\+\\sigma^\{2\}I\. Solving by the Schur complement, we obtain the closed\-form solution for optimal𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}at the current stage:

\(56\)λ\\displaystyle\\lambda=\(D​\(K~−1\+RΓ⊤​RΓ\)−1​D⊤\)−1​\(D​\(K~−1\+RΓ⊤​RΓ\)−1​rt−rb\),\\displaystyle=\\bigl\(D\\left\(\\tilde\{K\}^\{\-1\}\+R\_\{\\Gamma\}^\{\\top\}R\_\{\\Gamma\}\\right\)^\{\-1\}D^\{\\top\}\\bigr\)^\{\-1\}\\bigl\(D\\left\(\\tilde\{K\}^\{\-1\}\+R\_\{\\Gamma\}^\{\\top\}R\_\{\\Gamma\}\\right\)^\{\-1\}r\_\{t\}\-r\_\{b\}\\bigr\),\(57\)𝐅^gp,⋆\\displaystyle\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\}=\(K~−1\+RΓ⊤​RΓ\)−1​\(rt−D⊤​λ\),\\displaystyle=\\left\(\\tilde\{K\}^\{\-1\}\+R\_\{\\Gamma\}^\{\\top\}R\_\{\\Gamma\}\\right\)^\{\-1\}\\bigl\(r\_\{t\}\-D^\{\\top\}\\lambda\\bigr\),wherert=RΓ​\(W\)⊤​\(𝐅Γ,data−RΓ​\(𝐌1−1​D⊤​𝐮^−𝐌1−1​𝐠D\)\)r\_\{t\}=R\_\{\\Gamma\}\(W\)^\{\\top\}\\Bigl\(\\mathbf\{F\}\_\{\\Gamma,\\text\{data\}\}\-R\_\{\\Gamma\}\\bigl\(\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\-\\mathbf\{M\}\_\{1\}^\{\-1\}\\mathbf\{g\}\_\{D\}\\bigr\)\\Bigr\)andrb=𝐛−D​𝐌1−1​D⊤​𝐮^\.r\_\{b\}=\\mathbf\{b\}\-D\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\.

Then with the optimal𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}, we optimize the parameters for GP kernels and noise with the Adam optimizer\. The choice of kernels varies for different cases\. For example, for the independent GP, we choose the RBF kernel matrixKi​jK\_\{ij\}for a given edgeei​je\_\{ij\}, defined as

Ki​j​\(ei​j\)=exp⁡\(−‖𝐮i−𝐮j‖22​lei​j2\),K\_\{ij\}\(e\_\{ij\}\)=\\exp\\left\(\-\\frac\{\\\|\\mathbf\{u\}\_\{i\}\-\\mathbf\{u\}\_\{j\}\\\|^\{2\}\}\{2l\_\{e\_\{ij\}\}^\{2\}\}\\right\),where the length scalellandσ\\sigmaare the trainable parameters, i\.e\.,θG​P=\{l,σ\}\\theta\_\{GP\}=\\\{l,\\sigma\\\}\. Specifically, in numerical experiments, we choose the same initialization on all edges\. For the shared GP scenario, we consider the Automatic Relevance Determination \(ARD\) Matérn kernel\[Rasmussen2006\]withν=32\\nu=\\tfrac\{3\}\{2\}:

Kν=3/2​\(x,x′\)=σ2​\(1\+3​r\)​exp⁡\(−3​r\),r2=∑d\(xd−xd′\)2ℓd2\.K\_\{\\nu=3/2\}\(x,x^\{\\prime\}\)=\\sigma^\{2\}\\left\(1\+\\sqrt\{3\}r\\right\)\\exp\(\-\\sqrt\{3\}r\),\\quad r^\{2\}=\\sum\_\{d\}\\frac\{\(x\_\{d\}\-x^\{\\prime\}\_\{d\}\)^\{2\}\}\{\\ell\_\{d\}^\{2\}\}\.Here,ℓd\>0\\ell\_\{d\}\>0is a trainable length scale associated with thedd\-th input feature, where a smallℓd\\ell\_\{d\}indicates that the predicted flux varies strongly in thedd\-th feature direction\. This ARD parameterization allows the GP to have anisotropic scaling across different physical or geometric features\.

For the outer loop, with𝐅^gp\\widehat\{\\mathbf\{F\}\}^\{\\rm gp\}andθG​P\\theta\_\{GP\}fixed, we optimize over𝐮^\\widehat\{\\mathbf\{u\}\}andθW\\theta\_\{W\}with the Shampoo optimizer\[gupta2018shampoo\]by minimizing:

\(58\)ℒo​u​t​e​r=‖𝐅Γ,data−RΓ​\(𝐌1−1​D⊤​𝐮^\+𝐅^gp,⋆−𝐌1−1​𝐠D\)‖2\+‖𝐮data−𝐮^⊤​Ψ0‖2s\.t\.D​𝐌1−1​D⊤​𝐮^\+D​𝐅^gp,⋆=𝐛^\.\\begin\{split\}\\mathcal\{L\}\_\{outer\}&=\\left\\\|\\mathbf\{F\}\_\{\\Gamma,\\text\{data\}\}\-R\_\{\\Gamma\}\\left\(\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\+\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\}\-\\mathbf\{M\}\_\{1\}^\{\-1\}\\mathbf\{g\}\_\{D\}\\right\)\\right\\\|^\{2\}\\\\ &\\qquad\+\\left\\\|\\mathbf\{u\}\_\{\\text\{data\}\}\-\\widehat\{\\mathbf\{u\}\}^\{\\top\}\\Psi\_\{0\}\\right\\\|^\{2\}\\\\ \\text\{s\.t\.\}&\\quad D\\mathbf\{M\}\_\{1\}^\{\-1\}D^\{\\top\}\\widehat\{\\mathbf\{u\}\}\+D\\widehat\{\\mathbf\{F\}\}^\{\\rm gp,\\star\}=\\hat\{\\mathbf\{b\}\}\.\\\\ \\end\{split\}Essentially, we cast the optimal recovery problem on the reduced space subject to the equality constraint that enforces the conservation law exactly\. Under this setting, we have a constrained optimization problem where the KKT conditions yield a saddle\-point structure that we can use to derive a fast solution with the Schur complement\. The bilevel training is not strictly necessary, but can be useful for efficient training when variables are coupled and different model components, here the GP and transformer, are trained jointly\.

### 4\.3Geometric embedding forH​\(div\)H\(\\mathrm\{div\}\)coarse representations

To incorporate geometric information about the coarse degrees of freedom into the GP, we introduce a geometric embedding based on domain integrals of the coarse Whitney basis functions\. This embedding is designed to encode both the coarse variables in 0\- and 1\-forms in a manner consistent with the underlying FEEC structure\. Because scalarP0P\_\{0\}coefficients represent coarse cellwise quantities whileH​\(div\)H\(\\mathrm\{div\}\)coefficients represent flux degrees of freedom, we use different geometric summaries for 0\-forms and 1\-forms\.

First, we consider the coarseP0P\_\{0\}geometric embedding

###### Proposition 4\.1\.

Let\{ψi0\}i=1N0\\\{\\psi^\{0\}\_\{i\}\\\}\_\{i=1\}^\{N\_\{0\}\}denote the coarse 0\-form basis functions constructed as in[Proposition3\.1](https://arxiv.org/html/2606.11650#S3.Thmtheorem1)\. Fori=1,…,N0i=1,\\ldots,N\_\{0\}, substituting the definition ofψi0\\psi^\{0\}\_\{i\}gives the explicit representation for the geometric embedding on theiith PoU as

\(59\)si≔∫Ωψi0​\(x\)​𝑑x=∑a∈𝒞hWi​a​∫ΩϕaP0​\(x\)​𝑑x\.s\_\{i\}\\coloneqq\\int\_\{\\Omega\}\\psi\_\{i\}^\{0\}\(x\)\\,dx=\\sum\_\{a\\in\\mathcal\{C\}\_\{h\}\}W\_\{ia\}\\int\_\{\\Omega\}\\phi\_\{a\}^\{P\_\{0\}\}\(x\)\\,dx\.

This embedding defines a positive coarse area scalar that represents the effective measure of the coarse region associated with theii\-th PoU component and does not encode any directional information\. In terms of the GP input,sis\_\{i\}should be considered the geometric information for coarse nodeii\.

We next define the coarseR​T0RT\_\{0\}geometric embedding that accounts for both direction and magnitude\.

###### Proposition 4\.2\.

Let\{ψk1\}k=1N1\\\{\\psi^\{1\}\_\{k\}\\\}\_\{k=1\}^\{N\_\{1\}\}denote the coarseR​T0RT\_\{0\}basis functions constructed as in[Proposition3\.2](https://arxiv.org/html/2606.11650#S3.Thmtheorem2)\. We define the vector\-valued geometric embedding associated with thekk\-th coarse flux degree of freedom as

sk≔∫Ωψk1​\(x\)​𝑑x∈ℝ2\.s\_\{k\}\\coloneqq\\int\_\{\\Omega\}\\psi^\{1\}\_\{k\}\(x\)\\,dx\\;\\in\\;\\mathbb\{R\}^\{2\}\.
Substituting the expansion ofψk1\\psi^\{1\}\_\{k\}, we obtain:

sk=∑eWk​K​\(e\)​∫ΩφeR​T0​\(x\)​𝑑x\.s\_\{k\}=\\sum\_\{e\}W\_\{kK\(e\)\}\\int\_\{\\Omega\}\\varphi^\{RT\_\{0\}\}\_\{e\}\(x\)\\,dx\.
For each fine edgeeeshared by adjacent cellsKL​\(e\)K\_\{L\}\(e\)andKR​\(e\)K\_\{R\}\(e\), theR​T0RT\_\{0\}basis function satisfies:

∫KL​\(e\)∪KR​\(e\)φeR​T0​\(x\)​𝑑x=ce​n^e​\|e\|,\\int\_\{K\_\{L\}\(e\)\\cup K\_\{R\}\(e\)\}\\varphi^\{RT\_\{0\}\}\_\{e\}\(x\)\\,dx=c\_\{e\}\\,\\hat\{n\}\_\{e\}\\,\|e\|,wheren^e\\hat\{n\}\_\{e\}is the unit normal vector associated with the edge orientation,\|e\|\|e\|is the edge length, andcec\_\{e\}depends on the reference element geometry\.

From the above expression,𝐬k\\mathbf\{s\}\_\{k\}is a weighted sum of edge\-normal vectorsn^e\\hat\{n\}\_\{e\}scaled by edge lengths\|e\|\|e\|, and it encodes both the net orientation and the effective magnitude of the coarse flux degree of freedomψk1\\psi^\{1\}\_\{k\}\. Specifically, its direction reflects the dominant orientation of the fine\-scale flux contributions aggregated intoψk1\\psi^\{1\}\_\{k\}, while its magnitude reflects the total weighted edge measure\. This provides a geometric characterization of the coarse flux degree of freedom\. Therefore, the above geometric embeddings are consistent with the structure ofH​\(div\)H\(\\mathrm\{div\}\)spaces and enable the GP to incorporate geometric and directional information without relying on explicit spatial coordinates\.

### 4\.4Posterior error bound for boundary flux functionals

The quantified uncertainty below is conditional on the learned representation from PoU Whitney forms and kernel hyperparameters\. Under this conditioning, the boundary flux is a linear functional of the GP output at the coarse level and admits a closed\-form posterior error bound\. Previous work establishes a pointwise posterior error bound under the assumption that the GP input locations are fixed and known exactly\[propp2026discovery\]\. In this work, we extend this result to weighted linear functionals of the GP posterior, which is the relevant setting for boundary flux quantities\.

###### Proposition 4\.3\.

LetRΓ​\(W\)R\_\{\\Gamma\}\(W\)denote the coarse\-to\-fine reconstruction operator determined byWWin the proposed method \([53](https://arxiv.org/html/2606.11650#S4.E53)\)\. We define the flux functional on a given boundaryΓ\\Gammaand its estimator by

\(60\)ℱΓ≔Lv​\(f\)≔v⊤​f​\(Z\)≔∑m=1Mvm​f​\(zm\),ℱ^Γ≔v⊤​f^​\(Z\),\\begin\{split\}\\mathcal\{F\}\_\{\\Gamma\}\\coloneqq L\_\{v\}\(f\)&\\coloneqq v^\{\\top\}f\(Z\)\\coloneqq\\sum\_\{m=1\}^\{M\}v\_\{m\}\\,f\(z\_\{m\}\),\\\\ \\widehat\{\\mathcal\{F\}\}\_\{\\Gamma\}&\\coloneqq v^\{\\top\}\\hat\{f\}\(Z\),\\end\{split\}where the weightv⊤v^\{\\top\}is the row vector ofRΓ​\(W\)R\_\{\\Gamma\}\(W\), i\.e\.,vm=\(RΓ​\(W\)\)m,:v\_\{m\}=\(R\_\{\\Gamma\}\(W\)\)\_\{m,:\}, and the error is defined by

\(61\)e≔ℱΓ−ℱ^Γ\.e\\coloneqq\\mathcal\{F\}\_\{\\Gamma\}\-\\widehat\{\\mathcal\{F\}\}\_\{\\Gamma\}\.

###### Theorem 4\.4\.

Let𝒳⊆ℝp\\mathcal\{X\}\\subseteq\\mathbb\{R\}^\{p\}be the Gaussian process input space, letK:𝒳×𝒳→ℝK:\\mathcal\{X\}\\times\\mathcal\{X\}\\to\\mathbb\{R\}be a positive definite kernel with associated RKHSℋK\\mathcal\{H\}\_\{K\}, and letf∈ℋKf\\in\\mathcal\{H\}\_\{K\}\. LetX=\(x1,…,xN\)∈𝒳N,Z=\(z1,…,zM\)∈𝒳M,X=\(x\_\{1\},\\ldots,x\_\{N\}\)\\in\\mathcal\{X\}^\{N\},\\,Z=\(z\_\{1\},\\ldots,z\_\{M\}\)\\in\\mathcal\{X\}^\{M\},wherexi,zm∈ℝpx\_\{i\},z\_\{m\}\\in\\mathbb\{R\}^\{p\}and assume we have noisy data𝒟=\(X,Y\)\\mathcal\{D\}=\(X,Y\)withY=f​\(X\)\+εY=f\(X\)\+\\varepsilonwhereε∼𝒩​\(0,σε2​I\)\\varepsilon\\sim\\mathcal\{N\}\(0,\\sigma\_\{\\varepsilon\}^\{2\}I\)\. Letf^∈ℋK\\hat\{f\}\\in\\mathcal\{H\}\_\{K\}denote the minimizer of the optimal recovery problem:

f^​\(⋅\)=K​\(⋅,X\)​\(K​\(X,X\)\+σε2​I\)−1​Y\.\\hat\{f\}\(\\cdot\)=K\(\\cdot,X\)\\bigl\(K\(X,X\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\\bigr\)^\{\-1\}Y\.Define the posterior covariance onZZby

ΣZ≔K​\(Z,Z\)−K​\(Z,X\)​\(K​\(X,X\)\+σε2​I\)−1​K​\(X,Z\)\.\\Sigma\_\{Z\}\\coloneqq K\(Z,Z\)\-K\(Z,X\)\\Bigl\(K\(X,X\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\\Bigr\)^\{\-1\}K\(X,Z\)\.Forv∈ℝMv\\in\\mathbb\{R\}^\{M\}, consider the weighted functionalsℱΓ,ℱ^Γ\\mathcal\{F\}\_\{\\Gamma\},\\,\\hat\{\\mathcal\{F\}\}\_\{\\Gamma\}as defined in[Proposition4\.3](https://arxiv.org/html/2606.11650#S4.Thmtheorem3)\. Then we have the bounded mean\-squared error satisfying:

\(62\)𝔼ε​\[e2\]≤‖f‖ℋK2​v⊤​ΣZ​v\+σε2​‖v⊤​KZ​X​A−1‖22\.\\mathbb\{E\}\_\{\\varepsilon\}\[e^\{2\}\]\\leq\\\|f\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}v^\{\\top\}\\Sigma\_\{Z\}v\+\\sigma\_\{\\varepsilon\}^\{2\}\\bigl\\\|v^\{\\top\}K\_\{ZX\}A^\{\-1\}\\bigr\\\|\_\{2\}^\{2\}\.If we consider the diagonal bound, then \([62](https://arxiv.org/html/2606.11650#S4.E62)\) can be written as

\(63\)𝔼ε​\[e2\]≤‖f‖ℋK2​\(∑m=1M\|vm\|​σ​\(zm\)\)2\+σε2​‖∑m=1Mvm​K​\(zm,X\)​\(KX​X\+σε2​I\)−1‖22,\\mathbb\{E\}\_\{\\varepsilon\}\[e^\{2\}\]\\leq\\\|f\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}\\left\(\\sum\_\{m=1\}^\{M\}\|v\_\{m\}\|\\,\\sigma\(z\_\{m\}\)\\right\)^\{2\}\+\\sigma\_\{\\varepsilon\}^\{2\}\\left\\\|\\sum\_\{m=1\}^\{M\}v\_\{m\}K\(z\_\{m\},X\)\(K\_\{XX\}\+\\sigma\_\{\\varepsilon\}^\{2\}I\)^\{\-1\}\\right\\\|\_\{2\}^\{2\},whereσ2​\(zm\)≔\(ΣZ\)m​m\\sigma^\{2\}\(z\_\{m\}\)\\coloneqq\(\\Sigma\_\{Z\}\)\_\{mm\}\.

###### Proof 4\.5\.

LetA=K​\(X,X\)\+σε2​IA=K\(X,X\)\+\\sigma\_\{\\varepsilon\}^\{2\}I\. For finite ordered setsZ=\(zi\)i=1MZ=\(z\_\{i\}\)\_\{i=1\}^\{M\}andX=\(xj\)j=1NX=\(x\_\{j\}\)\_\{j=1\}^\{N\}, define:

KZ​X≔\(K​\(zi,xj\)\)i,j\.K\_\{ZX\}\\coloneqq\\bigl\(K\(z\_\{i\},x\_\{j\}\)\\bigr\)\_\{i,j\}\.SinceY=f​\(X\)\+εY=f\(X\)\+\\varepsilon, we can decompose the error as

\(64\)e:=ℱ−ℱ^=v⊤​f​\(Z\)−v⊤​KZ​X​A−1​Y=v⊤​f​\(Z\)−v⊤​KZ​X​A−1​f​\(X\)⏟=⁣:edet,v​\(f\)−v⊤​KZ​X​A−1​ε⏟=⁣:enoise,v\.\\begin\{split\}e:&=\\mathcal\{F\}\-\\widehat\{\\mathcal\{F\}\}\\\\ &=v^\{\\top\}f\(Z\)\-v^\{\\top\}K\_\{ZX\}A^\{\-1\}Y\\\\ &=\\underbrace\{v^\{\\top\}f\(Z\)\-v^\{\\top\}K\_\{ZX\}A^\{\-1\}f\(X\)\}\_\{=:\\,e\_\{\\text\{det\},v\}\(f\)\}\-\\underbrace\{v^\{\\top\}K\_\{ZX\}A^\{\-1\}\\varepsilon\}\_\{=:\\,e\_\{\\text\{noise\},v\}\}\.\\end\{split\}We then evaluate the MSE:

\(65\)MSE​\(e\)\\displaystyle\\mathrm\{MSE\}\(e\)=𝔼ε​\[e2\]=𝔼ε​\[\(edet−enoise\)2\]\\displaystyle=\\mathbb\{E\}\_\{\\varepsilon\}\[e^\{2\}\]=\\mathbb\{E\}\_\{\\varepsilon\}\\bigl\[\(e\_\{\\det\}\-e\_\{\\text\{noise\}\}\)^\{2\}\\bigr\]\(66\)=edet2\+Var​\(enoise\)\.\\displaystyle=e\_\{\\det\}^\{2\}\+\\mathrm\{Var\}\(e\_\{\\text\{noise\}\}\)\.For the noise part, sinceε∼𝒩​\(0,σε2​IN\)\\varepsilon\\sim\\mathcal\{N\}\(0,\\sigma\_\{\\varepsilon\}^\{2\}I\_\{N\}\), we have:

\(67\)𝔼ε​\[enoise,v\]=0,Var⁡\[enoise,v\]=σε2​‖v⊤​KZ​X​A−1‖22\.\\mathbb\{E\}\_\{\\varepsilon\}\[e\_\{\\text\{noise\},v\}\]=0,\\qquad\\operatorname\{Var\}\[e\_\{\\text\{noise\},v\}\]=\\sigma\_\{\\varepsilon\}^\{2\}\\\|v^\{\\top\}K\_\{ZX\}A^\{\-1\}\\\|\_\{2\}^\{2\}\.By the reproducing property,

\(68\)f​\(zm\)=⟨f,K​\(zm,⋅\)⟩ℋK,f​\(xi\)=⟨f,K​\(xi,⋅\)⟩ℋK\.f\(z\_\{m\}\)=\\langle f,K\(z\_\{m\},\\cdot\)\\rangle\_\{\\mathcal\{H\}\_\{K\}\},\\qquad f\(x\_\{i\}\)=\\langle f,K\(x\_\{i\},\\cdot\)\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\.Therefore, we further write the deterministic error in \([64](https://arxiv.org/html/2606.11650#S4.E64)\) as

\(69\)edet,v​\(f\)=⟨f,∑m=1Mvm​\[K​\(zm,⋅\)−K​\(zm,X\)​A−1​K​\(X,⋅\)\]⟩ℋK=⟨f,rv⟩ℋK\.e\_\{\\text\{det\},v\}\(f\)=\\left\\langle f,\\sum\_\{m=1\}^\{M\}v\_\{m\}\\Bigl\[K\(z\_\{m\},\\cdot\)\-K\(z\_\{m\},X\)A^\{\-1\}K\(X,\\cdot\)\\Bigr\]\\right\\rangle\_\{\\mathcal\{H\}\_\{K\}\}=\\langle f,r\_\{v\}\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\.Thus

\(70\)𝔼ε​\[e2\]=⟨f,rv⟩ℋK2\+σε2​‖v⊤​KZ​X​A−1‖22\.\\mathbb\{E\}\_\{\\varepsilon\}\[e^\{2\}\]=\\langle f,r\_\{v\}\\rangle\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}\+\\sigma\_\{\\varepsilon\}^\{2\}\\bigl\\\|v^\{\\top\}K\_\{ZX\}A^\{\-1\}\\bigr\\\|\_\{2\}^\{2\}\.By applying the Cauchy–Schwarz inequality, we have:

\(71\)\|⟨f,rv⟩ℋK\|2≤‖f‖ℋK2​‖rv‖ℋK2\.\\displaystyle\|\\langle f,r\_\{v\}\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\|^\{2\}\\leq\\\|f\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}\\\|r\_\{v\}\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}\.We now compute‖rv‖ℋK2\\\|r\_\{v\}\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}\. Recall the reproducing property that for ordered setsZ=\(zi\)i=1MZ=\(z\_\{i\}\)\_\{i=1\}^\{M\}andX=\(xj\)j=1NX=\(x\_\{j\}\)\_\{j=1\}^\{N\}, and vectorsa∈ℝM,b∈ℝNa\\in\\mathbb\{R\}^\{M\},\\,b\\in\\mathbb\{R\}^\{N\},

\(72\)⟨a⊤​K​\(Z,⋅\),b⊤​K​\(X,⋅\)⟩ℋK\\displaystyle\\left\\langle a^\{\\top\}K\(Z,\\cdot\),b^\{\\top\}K\(X,\\cdot\)\\right\\rangle\_\{\\mathcal\{H\}\_\{K\}\}=∑i=1M∑j=1Nai​bj​⟨K​\(zi,⋅\),K​\(xj,⋅\)⟩ℋK\\displaystyle=\\sum\_\{i=1\}^\{M\}\\sum\_\{j=1\}^\{N\}a\_\{i\}b\_\{j\}\\langle K\(z\_\{i\},\\cdot\),K\(x\_\{j\},\\cdot\)\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\(73\)=∑i=1M∑j=1Nai​bj​K​\(zi,xj\)=a⊤​KZ​X​b\.\\displaystyle=\\sum\_\{i=1\}^\{M\}\\sum\_\{j=1\}^\{N\}a\_\{i\}b\_\{j\}K\(z\_\{i\},x\_\{j\}\)=a^\{\\top\}K\_\{ZX\}b\.Thenrv​\(⋅\)r\_\{v\}\(\\cdot\)can be equivalently written as

\(74\)rv​\(⋅\)=v⊤​K​\(Z,⋅\)−α⊤​K​\(X,⋅\),r\_\{v\}\(\\cdot\)=v^\{\\top\}K\(Z,\\cdot\)\-\\alpha^\{\\top\}K\(X,\\cdot\),where

α⊤\\displaystyle\\alpha^\{\\top\}=v⊤​KZ​X​A−1\\displaystyle=v^\{\\top\}K\_\{ZX\}A^\{\-1\}K​\(Z,⋅\)\\displaystyle K\(Z,\\cdot\)=\(K​\(z1,⋅\),…,K​\(zM,⋅\)\)⊤,\\displaystyle=\\bigl\(K\(z\_\{1\},\\cdot\),\\ldots,K\(z\_\{M\},\\cdot\)\\bigr\)^\{\\top\},K​\(X,⋅\)\\displaystyle K\(X,\\cdot\)=\(K​\(x1,⋅\),…,K​\(xN,⋅\)\)⊤\\displaystyle=\\bigl\(K\(x\_\{1\},\\cdot\),\\ldots,K\(x\_\{N\},\\cdot\)\\bigr\)^\{\\top\}Applying \([73](https://arxiv.org/html/2606.11650#S4.E73)\) yields:

\(75\)‖rv‖ℋK2\\displaystyle\\\|r\_\{v\}\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}=⟨v⊤​K​\(Z,⋅\)−α⊤​K​\(X,⋅\),v⊤​K​\(Z,⋅\)−α⊤​K​\(X,⋅\)⟩ℋK\\displaystyle=\\left\\langle v^\{\\top\}K\(Z,\\cdot\)\-\\alpha^\{\\top\}K\(X,\\cdot\),v^\{\\top\}K\(Z,\\cdot\)\-\\alpha^\{\\top\}K\(X,\\cdot\)\\right\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\(76\)=⟨v⊤​K​\(Z,⋅\),v⊤​K​\(Z,⋅\)⟩ℋK−2​⟨v⊤​K​\(Z,⋅\),α⊤​K​\(X,⋅\)⟩ℋK\\displaystyle=\\left\\langle v^\{\\top\}K\(Z,\\cdot\),v^\{\\top\}K\(Z,\\cdot\)\\right\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\-2\\left\\langle v^\{\\top\}K\(Z,\\cdot\),\\alpha^\{\\top\}K\(X,\\cdot\)\\right\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\(77\)\+⟨α⊤​K​\(X,⋅\),α⊤​K​\(X,⋅\)⟩ℋK\\displaystyle\\qquad\+\\left\\langle\\alpha^\{\\top\}K\(X,\\cdot\),\\alpha^\{\\top\}K\(X,\\cdot\)\\right\\rangle\_\{\\mathcal\{H\}\_\{K\}\}\(78\)=v⊤​KZ​Z​v−2​v⊤​KZ​X​A−1​KX​Z​v\+v⊤​KZ​X​A−1​KX​X​A−1​KX​Z​v\\displaystyle=v^\{\\top\}K\_\{ZZ\}v\-2v^\{\\top\}K\_\{ZX\}A^\{\-1\}K\_\{XZ\}v\+v^\{\\top\}K\_\{ZX\}A^\{\-1\}K\_\{XX\}A^\{\-1\}K\_\{XZ\}v\(79\)=v⊤​KZ​Z​v−v⊤​KZ​X​A−1​KX​Z​v−σε2​v⊤​KZ​X​A−2​KX​Z​v\\displaystyle=v^\{\\top\}K\_\{ZZ\}v\-v^\{\\top\}K\_\{ZX\}A^\{\-1\}K\_\{XZ\}v\-\\sigma\_\{\\varepsilon\}^\{2\}v^\{\\top\}K\_\{ZX\}A^\{\-2\}K\_\{XZ\}v\(80\)=v⊤​ΣZ​v−σε2​v⊤​KZ​X​A−2​KX​Z​v\\displaystyle=v^\{\\top\}\\Sigma\_\{Z\}v\-\\sigma\_\{\\varepsilon\}^\{2\}v^\{\\top\}K\_\{ZX\}A^\{\-2\}K\_\{XZ\}v\(81\)=v⊤​ΣZ​v−σε2​‖A−1​KX​Z​v‖22≤v⊤​ΣZ​v\\displaystyle=v^\{\\top\}\\Sigma\_\{Z\}v\-\\sigma\_\{\\varepsilon\}^\{2\}\\left\\\|A^\{\-1\}K\_\{XZ\}v\\right\\\|\_\{2\}^\{2\}\\leq v^\{\\top\}\\Sigma\_\{Z\}vsince‖A−1​KX​Z​v‖22≥0\\left\\\|A^\{\-1\}K\_\{XZ\}v\\right\\\|\_\{2\}^\{2\}\\geq 0forKX​X≥0K\_\{XX\}\\geq 0and symmetric positive definiteAA\.

Hence with \([70](https://arxiv.org/html/2606.11650#S4.E70)\) and \([81](https://arxiv.org/html/2606.11650#S4.E81)\), we obtain:

\(82\)𝔼ε​\[e2\]≤‖f‖ℋK2​v⊤​ΣZ​v\+σε2​‖v⊤​KZ​X​A−1‖22\.\\mathbb\{E\}\_\{\\varepsilon\}\[e^\{2\}\]\\leq\\\|f\\\|\_\{\\mathcal\{H\}\_\{K\}\}^\{2\}v^\{\\top\}\\Sigma\_\{Z\}v\+\\sigma\_\{\\varepsilon\}^\{2\}\\bigl\\\|v^\{\\top\}K\_\{ZX\}A^\{\-1\}\\bigr\\\|\_\{2\}^\{2\}\.The decomposition of error into noise and deterministic components is readily interpretable\. The deterministic component captures how well the training data atXXcould reconstruct the functionalv⊤​f​\(Z\)v^\{\\top\}f\(Z\)if the observations were noiseless\. This error is largely controlled by the RKHS norm \(or complexity\) of the true functionffand by the posterior covariance at evaluation pointsZZ\. The noise component, controlled byσε2​‖v⊤​KZ​X​A−1‖22\\sigma\_\{\\varepsilon\}^\{2\}\\\|v^\{\\top\}K\_\{ZX\}A^\{\-1\}\\\|\_\{2\}^\{2\}, captures how much the observation noise is amplified through the kernel reconstruction and through the linear functionalv⊤v^\{\\top\}\.

Finally, we derive the diagonal posterior standard deviation from the full covariance form\. Recall thatΣZ\\Sigma\_\{Z\}is a positive semidefinite block matrix\. We therefore have

\(83\)\|\(ΣZ\)m​n\|≤\(ΣZ\)m​m​\(ΣZ\)n​n=σ​\(zm\)​σ​\(zn\)\.\|\(\\Sigma\_\{Z\}\)\_\{mn\}\|\\leq\\sqrt\{\(\\Sigma\_\{Z\}\)\_\{mm\}\}\\sqrt\{\(\\Sigma\_\{Z\}\)\_\{nn\}\}=\\sigma\(z\_\{m\}\)\\sigma\(z\_\{n\}\)\.Thus, we have:

\(84\)v⊤​ΣZ​v\\displaystyle v^\{\\top\}\\Sigma\_\{Z\}v=∑m,n=1Mvm​vn​\(ΣZ\)m​n\\displaystyle=\\sum\_\{m,n=1\}^\{M\}v\_\{m\}v\_\{n\}\(\\Sigma\_\{Z\}\)\_\{mn\}≤∑m,n=1M\|vm\|​\|vn\|​σ​\(zm\)​σ​\(zn\)\\displaystyle\\leq\\sum\_\{m,n=1\}^\{M\}\|v\_\{m\}\|\|v\_\{n\}\|\\sigma\(z\_\{m\}\)\\sigma\(z\_\{n\}\)=\(∑m=1M\|vm\|​σ​\(zm\)\)2\\displaystyle=\\left\(\\sum\_\{m=1\}^\{M\}\|v\_\{m\}\|\\sigma\(z\_\{m\}\)\\right\)^\{2\}Substituting \([84](https://arxiv.org/html/2606.11650#S4.E84)\) into \([62](https://arxiv.org/html/2606.11650#S4.E62)\) gives the desired error bound in the diagonal form[Equation63](https://arxiv.org/html/2606.11650#S4.E63)\. This bound can be interpreted as yielding large errors if the functional places a large weight on uncertain locations\.

## 5Results

In this section, we demonstrate the performance of our proposed method on various problems, including the 1D Poisson equation with analytic solution \([Section5\.1](https://arxiv.org/html/2606.11650#S5.SS1)\), the nonlinear diffusion equation with manufactured solution \([Section5\.2](https://arxiv.org/html/2606.11650#S5.SS2)\), and the 2D advection\-diffusion equation on a bell\-shaped complex geometry \([Section5\.3](https://arxiv.org/html/2606.11650#S5.SS3)\)\. We also present results for a semiconductorp−np\-ndiode problem \([Section5\.4](https://arxiv.org/html/2606.11650#S5.SS4)\) to illustrate the range over which the uncertainty estimates remain reliable for the proposed method\. For the training of data\-driven Whitney forms, we consider a lightweight cross\-attention transformer with dropout\.

### 5\.1Toy one\-dimensional linear example

Consider a toy example of the one\-dimensional Poisson equation onΩ=\(0,1\)\\Omega=\(0,1\):

\(85\)F=−∇u,∇⋅F=f,F=\-\\nabla u,\\qquad\\nabla\\cdot F=f,with Dirichlet boundary conditionsu​\(0\)=αu\(0\)=\\alpha,u​\(1\)=0u\(1\)=0, and constant source termf≡1f\\equiv 1\. The analytic solution is:

\(86\)u​\(x\)=−x22\+\(12−α\)​x\+α,F​\(x\)=x\+α−12\.u\(x\)=\-\\tfrac\{x^\{2\}\}\{2\}\+\\bigl\(\\tfrac\{1\}\{2\}\-\\alpha\\bigr\)x\+\\alpha,\\qquad F\(x\)=x\+\\alpha\-\\tfrac\{1\}\{2\}\.The goal is to learn the solutionuuand the associated fluxFFgiven the left boundary conditionα\\alpha\. In particular, we use this toy example because the flux depends linearly on the left boundary condition\.

In the implementation, to construct the fine\-scale discretization, we employ the classic Galerkin finite element withP0P\_\{0\}elements for solutionuuandP1P\_\{1\}Lagrange elements for the flux\. For a mesh withNNcells, we haveNc​e​l​l=NN\_\{cell\}=Ncell\-averaged values foruuandNf​a​c​e=N\+1N\_\{face\}=N\+1nodal flux values forFF\. Then the sparse basis evaluations in the fine scaleϕ0∈ℝNcell×Nq\\phi^\{0\}\\in\\mathbb\{R\}^\{N\_\{\\text\{cell\}\}\\times N\_\{q\}\}andϕ1∈ℝNface×Nq\\phi^\{1\}\\in\\mathbb\{R\}^\{N\_\{\\text\{face\}\}\\times N\_\{q\}\}are constructed with the quadrature weight exponent12\\frac\{1\}\{2\}\. In one dimension, this is analogous to the lowest\-order Raviart–Thomas spaceRT0\\mathrm\{RT\}\_\{0\}where each hat functionϕi1\\phi\_\{i\}^\{1\}is supported on the two cells adjacent to nodeii, and its derivatived​ϕi1/d​xd\\phi\_\{i\}^\{1\}/dxplays the role of∇⋅ϕiF\\nabla\\\!\\cdot\\\!\\phi\_\{i\}^\{F\}\. The resultingP0P\_\{0\}–P1P\_\{1\}pair yields the standard mixedP0P\_\{0\}–H​\(div\)H\(\\mathrm\{div\}\)formulation in 1D\. In addition, for a prescribed number of PoU, we haveN0=NPoUN\_\{0\}=N\_\{\\text\{PoU\}\}for the coarse 0\-forms andN1int=\(N02\)N\_\{1\}^\{\\text\{int\}\}=\\binom\{N\_\{0\}\}\{2\}for the interior coarse 1\-forms, respectively\. For the 1D problem, we explicitly set two boundary 1\-formsψ1left\\psi\_\{1\}^\{\\rm left\}andψ1right\\psi\_\{1\}^\{\\rm right\}, which gives theN1b​c=2N\_\{1\}^\{bc\}=2and total number of coarse 1\-formsN1=\(N02\)\+2N\_\{1\}=\\binom\{N\_\{0\}\}\{2\}\+2\.

Comparison of the ground truth with the prediction is shown in Fig\.[3](https://arxiv.org/html/2606.11650#S5.F3)\(left\)\. The reconstructed solution field from the coarse 0\-forms can approximate the ground\-truth solution with the MSE below10−510^\{\-5\}\. Fig\.[3](https://arxiv.org/html/2606.11650#S5.F3)\(right\) illustrates the GP inference of the left boundary flux with respect to different values ofα\\alpha\. For the 1D problem where the target boundaryΓ\\Gammais reduced to the left boundary pointx=0x=0, the coarse\-to\-fine reconstruction of flux in[Proposition4\.3](https://arxiv.org/html/2606.11650#S4.Thmtheorem3)is simplified to use the coarse boundary 1\-form for the left boundary\. Furthermore, as this problem is linear, the error bound is still tight near the training range even in the small extrapolation area, with only a slight increase in the most extreme cases,α=0\.5\\alpha=0\.5andα=3\.5\\alpha=3\.5\.

![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/1d_toy_solution.png)
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/1d_toy_posterior.png)

Figure 3:Solution reconstruction and quantified uncertainty for the toy example\.Left: Solution fielduuwith respect to different left boundary conditionsα\\alpha\. The black solid line is the analytic solution, and the red dashed line is the reconstructed solution\.Right: Inference result for the left boundary flux with respect to left boundary conditionsα\\alphawith shaded blue area as the error bound from Theorem[4\.4](https://arxiv.org/html/2606.11650#S4.Thmtheorem4)\.
### 5\.21D nonlinear example with manufactured solution

We next consider a nonlinear diffusion equation with piecewise nonlinear diffusivity and the same Dirichlet boundary conditions as in[Section5\.1](https://arxiv.org/html/2606.11650#S5.SS1):

\(87\)F=−k​\(u\)​∇u,∇⋅F=f,F=\-k\(u\)\\,\\nabla u,\\qquad\\nabla\\cdot F=f,where

\(88\)k​\(u\)=\{k0if​u≤u0,β​\(u−u0\)q\+k0if​u\>u0\.k\(u\)=\\begin\{cases\}k\_\{0\}&\\text\{if \}u\\leq u\_\{0\},\\\\ \\beta\\,\(u\-u\_\{0\}\)^\{q\}\+k\_\{0\}&\\text\{if \}u\>u\_\{0\}\.\\end\{cases\}The analytic solution isu​\(x\)=α​\(1−x\)u\(x\)=\\alpha\(1\-x\), which givesF​\(x\)=α​k​\(α​\(1−x\)\)F\(x\)=\\alpha\\,k\\bigl\(\\alpha\(1\-x\)\\bigr\)and

f​\(x\)=∇⋅F​\(x\)=\{0if​α​\(1−x\)≤u0,−β​q​α2​\(α​\(1−x\)−u0\)q−1if​α​\(1−x\)\>u0\.f\(x\)=\\nabla\\cdot F\(x\)=\\begin\{cases\}0&\\text\{if \}\\alpha\(1\-x\)\\leq u\_\{0\},\\\\ \-\\beta q\\alpha^\{2\}\\bigl\(\\alpha\(1\-x\)\-u\_\{0\}\\bigr\)^\{q\-1\}&\\text\{if \}\\alpha\(1\-x\)\>u\_\{0\}\.\\end\{cases\}This problem features a nonlinear constitutive lawk​\(u\)k\(u\), a spatially varying sourcef​\(x\)f\(x\)that depends on the parameterα\\alpha, and a kink inkkat the thresholdu=u0u=u\_\{0\}\. Therefore, we have a manufactured problem where the flux is a nonlinear function of the left boundary conditionα\\alphawith a transition point atx=1−u0αx=1\-\\frac\{u\_\{0\}\}\{\\alpha\}\.

![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/1d_nonlinear_pou.png)Figure 4:Conditional coarse bases for[Section5\.2](https://arxiv.org/html/2606.11650#S5.SS2)\. Coarse shape functions\{ψi0\}\\\{\\psi^\{0\}\_\{i\}\\\}evaluated on the fine\-scale nodes can adapt to different valuesα=\[1,2\.1,3\]\\alpha=\[1,2\.1,3\]\.![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/1d_nonlinear_posterior.png)Figure 5:Quantified uncertainty for the left\-boundary flux in the 1D nonlinear problem\.The learned surrogate predicts the boundary fluxFΓF\_\{\\Gamma\}as a function of the left Dirichlet boundary valueα\\alpha\. The growth of the uncertainty band outside the training range inα\\alphareflects reduced confidence in extrapolation\.The visualization of coarse shape functions\{ψj0\}\\\{\\psi\_\{j\}^\{0\}\\\}in Fig\.[4](https://arxiv.org/html/2606.11650#S5.F4)demonstrates that the conditional neural Whitney form allows an adaptive adjustment of the PoU\. Asα\\alphachanges, the coarse basis functions reallocate their support to better accommodate the condition\. One coarse basis function becomes much more dominant near the left side at larger values ofα\\alpha, suggesting that the solution structure there becomes more parameter\-sensitive\. For this example with both linear and nonlinear regions, the error bound grows more noticeably in the extrapolation range than in the previous linear example \(Fig\.[5](https://arxiv.org/html/2606.11650#S5.F5)\)\. Furthermore, in the left extrapolation region where the flux is still linear in the conditioning variable, the GP\-inferred flux deviates less from the ground\-truth than in the nonlinear extrapolation region on the right\. The error bound grows in both extrapolation directions because the distance in the GP feature space from the training samples increases regardless of whether the extrapolation is in the linear regime or the nonlinear regime\.

### 5\.3Conditioned 2D advection\-diffusion equation in complex geometry

We extend the problem to the steady advection\-diffusion equation in a complex geometry and an unstructured mesh:

\(89\)−∇⋅\(k​∇u\)\+𝜷⋅∇u=0in​Ω,\-\\nabla\\cdot\(k\\,\\nabla u\)\+\\boldsymbol\{\\beta\}\\cdot\\nabla u=0\\quad\\text\{in \}\\Omega,wherekkis the diffusion coefficient and𝜷=\(cos⁡θ,sin⁡θ\)\\boldsymbol\{\\beta\}=\(\\cos\\theta,\\sin\\theta\)is the advection direction withθ∈\[0,2​π\)\\theta\\in\[0,2\\pi\)\. The domain is a two\-dimensional “bell”\-shaped domainΩ\\Omegawith three distinct Dirichlet boundary conditions as shown in Fig\.[6](https://arxiv.org/html/2606.11650#S5.F6)\. We use an unstructured triangular mesh to discretize the complex bell geometry and parameterize the problem with bothθ\\thetaanduc​r​a​c​ku\_\{crack\}\. Specifically, we have the top handle of the bell \(Γ1\\Gamma\_\{1\}\) withu=1u=1, the crack boundary \(Γ2\\Gamma\_\{2\}\) withu=uc​r​a​c​ku=u\_\{crack\}, and the remaining bell\-body boundary \(Γ3\\Gamma\_\{3\}\) withu=0u=0\. The conditioning variableθ\\thetacontrols the advection direction and thus changes the transport pattern inside the domain\.

![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_mesh.png)Figure 6:Complex geometry demonstration\.Left figure shows the irregular triangular mesh and right figure shows the nonhomogeneous Dirichlet boundary conditions\.![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_solution.png)
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_loss.png)
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_post.png)

Figure 7:Quantified uncertainty on complex geometry\.Top: The solution to the advection\-diffusion equation, conditioned on the advection direction, is recovered\.Bottom:Dropout during training provides a high\-quality basis; turning off at 100k epochs allows refinement of the final basis \(bottom left\)\. The error in the flux through the crack at the bottom of the Liberty Bell is bounded by the estimator in[Equation63](https://arxiv.org/html/2606.11650#S4.E63)\(bottom right\)\.Results are shown in Fig\.[7](https://arxiv.org/html/2606.11650#S5.F7)\. In the implementation, we first train the transformer for the Whitney forms with a dropout rate of 0\.1 for 100,000 epochs and then turn off the dropout for another 20,000 epochs\. Empirically, this two\-stage training strategy alleviates overfitting and prevents the transformer from getting stuck in local minima during training\. Fig\.[8](https://arxiv.org/html/2606.11650#S5.F8)visualizes the Whitney 0\-forms learned by the neural PoU framework for three different values of the conditioning variables, showing that the learned partitions adapt naturally to the geometric structure\. Asθ\\thetaanduc​r​a​c​ku\_\{crack\}increase, some 0\-forms tend to localize near the crack interface to capture changes near the crack boundary efficiently, while others tend to align with the flow near other boundaries or within the domain\.

![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_pou_ex1.png)
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_pou_ex2.png)
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/2d_pou_ex3.png)

Figure 8:Conditioned basis adapts to advection\.The three basis functions are shown conditioned on anglesθ∈\[0\.45​π,0\.75​π,1\.05​π\]\\theta\\in\[0\.45\\pi,0\.75\\pi,1\.05\\pi\], respectively, highlighting the ability of the basis to adapt itself to the relevant conditioning\.
### 5\.4p−np\-ndiodes

Finally, we consider thep−np\-ndiode problem obtained from Charon TCAD simulations\[osti\_1575982\]on the rectangular device domainΩ=\(0,Lx\)×\(0,Ly\)\\Omega=\(0,L\_\{x\}\)\\times\(0,L\_\{y\}\)withLx=1L\_\{x\}=1andLy=0\.5L\_\{y\}=0\.5\. For each applied anode biasVaV\_\{a\}, the data consist of the electrostatic potential, electron and hole current densities, doping profiles, and the terminal currents at the cathode and anode\. The cathode bias is fixed and the anode bias is swept over\[−1,1\]\[\-1,1\]in thermal\-voltage units; multiplication byV0=2\.585×10−2​VV\_\{0\}=2\.585\\times 10^\{\-2\}\\,\\mathrm\{V\}gives the corresponding physical voltage scale\.

The two\-dimensional TCAD fields are reduced to one\-dimensional profiles in a way that preserves the terminal\-current constraint\. Letφ​\(x,y;Va\)\\varphi\(x,y;V\_\{a\}\)denote the electrostatic potential and letjx​\(x,y;Va\)j\_\{x\}\(x,y;V\_\{a\}\)denote the net conventional current density in the axial direction\. On the structured100×10100\\times 10cell grid, we define:

\(90\)u​\(x;Va\)=1Ly​∫0Lyφ​\(x,y;Va\)​𝑑y,J​\(x;Va\)=∫0Lyjx​\(x,y;Va\)​𝑑y\.u\(x;V\_\{a\}\)=\\frac\{1\}\{L\_\{y\}\}\\int\_\{0\}^\{L\_\{y\}\}\\varphi\(x,y;V\_\{a\}\)\\,dy,\\qquad J\(x;V\_\{a\}\)=\\int\_\{0\}^\{L\_\{y\}\}j\_\{x\}\(x,y;V\_\{a\}\)\\,dy\.Since the steady drift–diffusion solution satisfies charge conservation \(∇⋅j=0\\nabla\\cdot j=0\) and the top and bottom boundaries are insulating, integration over the cross section gives:

\(91\)dd​x​J​\(x;Va\)=0,J​\(x;Va\)≃Ia​\(Va\)\.\\frac\{d\}\{dx\}J\(x;V\_\{a\}\)=0,\\qquad J\(x;V\_\{a\}\)\\simeq I\_\{a\}\(V\_\{a\}\)\.
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/1d_diode_visualization.png)
![Refer to caption](https://arxiv.org/html/2606.11650v1/figures/1d_diode.png)

Figure 9:Quantified uncertainty for the I\-V response of a semiconductor component\.Top:schematic of the one\-dimensional p\-n diode and the prescribed anode–cathode orientation\.Bottom:Comparison of the ground\-truth current and predicted current in original scale \(bottom left\) and log scale \(bottom right\)\. Despite the appearance of a relatively good current reconstruction on a linear scale \(bottom left\), on a log scale the uncertainty estimator reveals the range of voltages over which the model may be reliably applied\.We use the same conditional Whitney–GP construction as in previous 1D examples, where the neural Whitney forms are generated by the transformer conditioned onVaV\_\{a\}\. Given that the current values vary over several orders of magnitude, we train the GP on the logarithm of the current magnitude and map the posterior mean back to the original scale by exponentiation, i\.e\., the posterior meanμ​\(Va\)\\mu\(V\_\{a\}\)is mapped back to the original current magnitude byexp⁡\(μ​\(Va\)\)\\exp\(\\mu\(V\_\{a\}\)\)\. The numerical results are shown in Fig\.[9](https://arxiv.org/html/2606.11650#S5.F9)with the inference in original scale on the left and in log scale on the right\. Despite the fairly good current reconstruction on a linear scale, the uncertainty estimator reveals the range of voltages over which the model may be reliably applied, and flags where the uncertainty becomes large relative to the scale of the quantity of interest\. This example should be considered an illustration and guideline for where the proposed model starts to become less trustworthy, i\.e\., the posterior uncertainty grows rapidly in the exponential regime\.

## 6Conclusion

In this work, we developed a structure\-preserving neural surrogate framework for PDE\-governed systems that incorporates a reduced finite element space and tractable uncertainty quantification with a closed\-form posterior error bound\. Our method separates the problem into two coupled parts: inP1, the trainable PoU Whitney forms generate the reduced\-orderH​\(div\)H\(\\mathrm\{div\}\)–L2L^\{2\}spaces that preserve the conservation law with the induced coarse divergence operator and support the graph interpretation; inP2, the GP learns the nonlinear state\-to\-flux map on the coarse graph by solving the optimal recovery problem with an equality constraint that enforces the conservation exactly\. In particular, the conditional Whitney forms provide an adaptive mechanism that enables generalization across different geometries, and the optimal recovery formulation yields a saddle\-point KKT condition with a fast Schur\-complement solve\. Numerical results demonstrate that the proposed method can reconstruct accurate solution fields and fluxes and provide meaningful uncertainty quantification for boundary flux functionals\. In future work, we will generalize the framework to more complex problems and real datasets for practical applications, such as large\-scale deployments for climate modeling\. In addition, the current structure\-preserving construction is based on the de Rham complex and is primarily designed for scalar\- and vector\-valued problems; extending the method to handle tensor\-valued variables is an important direction for future work\.

## Acknowledgments

H\. Zhang and Dr\. Kinch acknowledge support from the United States Department of Energy under the Advanced Scientific Computing Research \(ASCR\) program \(award number DE\-SC0024563\)\. Drs\. Kinch, Owhadi, Propp, and Trask acknowledge funding under the ASCR\-sponsored Mathematical Multifaceted Integrated Capability Centers program \(award number DE\-SC0023163\)\.

## References

## Appendix AProof for[Proposition3\.5](https://arxiv.org/html/2606.11650#S3.Thmtheorem5)

###### Proof A\.1\.

For a boundary edgee∈ℰγe\\in\\mathcal\{E\}\_\{\\gamma\}, the divergence∇⋅ϕe1\\nabla\\\!\\cdot\\boldsymbol\{\\phi\}\_\{e\}^\{1\}is supported entirely on its unique adjacent cellK​\(e\)K\(e\)and is zero everywhere else\. Letβeγ=\(∇⋅ϕe1,ϕK​\(e\)P0\)Ω\\beta\_\{e\}^\{\\gamma\}=\(\\nabla\\\!\\cdot\\boldsymbol\{\\phi\}\_\{e\}^\{1\},\\phi\_\{K\(e\)\}^\{P\_\{0\}\}\)\_\{\\Omega\}\. Since the coarse 0\-form and boundary 1\-form are defined as

ψi0=∑a∈𝒞hWi​a​ϕaP0,𝝍α,γ1,bc=∑e∈ℰγWα,K​\(e\)​ϕe1,\\psi\_\{i\}^\{0\}=\\sum\_\{a\\in\\mathcal\{C\}\_\{h\}\}W\_\{ia\}\\phi\_\{a\}^\{P\_\{0\}\},\\qquad\\boldsymbol\{\\psi\}\_\{\\alpha,\\gamma\}^\{1,\\mathrm\{bc\}\}=\\sum\_\{e\\in\\mathcal\{E\}\_\{\\gamma\}\}W\_\{\\alpha,K\(e\)\}\\boldsymbol\{\\phi\}\_\{e\}^\{1\},
we have

\(Dγ\)i​α\\displaystyle\(D\_\{\\gamma\}\)\_\{i\\alpha\}=\(∇⋅𝝍α,γ1,bc,ψi0\)Ω\\displaystyle=\(\\nabla\\\!\\cdot\\boldsymbol\{\\psi\}\_\{\\alpha,\\gamma\}^\{1,\\mathrm\{bc\}\},\\psi\_\{i\}^\{0\}\)\_\{\\Omega\}=∑e∈ℰγWα,K​\(e\)​\(∇⋅ϕe1,ψi0\)Ω\\displaystyle=\\sum\_\{e\\in\\mathcal\{E\}\_\{\\gamma\}\}W\_\{\\alpha,K\(e\)\}\(\\nabla\\\!\\cdot\\boldsymbol\{\\phi\}\_\{e\}^\{1\},\\psi\_\{i\}^\{0\}\)\_\{\\Omega\}=∑e∈ℰγWα,K​\(e\)​Wi,K​\(e\)​\(∇⋅ϕe1,ϕK​\(e\)P0\)Ω\\displaystyle=\\sum\_\{e\\in\\mathcal\{E\}\_\{\\gamma\}\}W\_\{\\alpha,K\(e\)\}W\_\{i,K\(e\)\}\(\\nabla\\\!\\cdot\\boldsymbol\{\\phi\}\_\{e\}^\{1\},\\phi\_\{K\(e\)\}^\{P\_\{0\}\}\)\_\{\\Omega\}=∑e∈ℰγβeγ​Wα,K​\(e\)​Wi,K​\(e\)\\displaystyle=\\sum\_\{e\\in\\mathcal\{E\}\_\{\\gamma\}\}\\beta\_\{e\}^\{\\gamma\}W\_\{\\alpha,K\(e\)\}W\_\{i,K\(e\)\}This showsDγ=Wγ​Bγ​Wγ⊤D\_\{\\gamma\}=W\_\{\\gamma\}B\_\{\\gamma\}W\_\{\\gamma\}^\{\\top\}, where\(Wγ\)i​e=Wi,K​\(e\),Bγ=diag⁡\{βeγ:e∈ℰγ\}\(W\_\{\\gamma\}\)\_\{ie\}=W\_\{i,K\(e\)\},\\ B\_\{\\gamma\}=\\operatorname\{diag\}\\\{\\beta\_\{e\}^\{\\gamma\}:e\\in\\mathcal\{E\}\_\{\\gamma\}\\\}\.

Then it remains to check the sufficient condition\. Let

D∂=\[D1​D2​⋯​Dr\]D\_\{\\partial\}=\[D\_\{1\}\\;D\_\{2\}\\;\\cdots\\;D\_\{r\}\]denote the boundary block of the divergence matrixDD, collecting the boundary groupsγ=1,…,r\\gamma=1,\\dots,r\. GivenD∈ℝN0×N1D\\in\\mathbb\{R\}^\{N\_\{0\}\\times N\_\{1\}\}, we have

rank⁡\(D\)≤N0\.\\operatorname\{rank\}\(D\)\\leq N\_\{0\}\.IfD∂D\_\{\\partial\}has full row rankN0N\_\{0\}, sinceD∂D\_\{\\partial\}is the column submatrix ofDDcorresponding to the boundary basis functions, we have

rank⁡\(D\)≥rank⁡\(D∂\)=N0,\\operatorname\{rank\}\(D\)\\geq\\operatorname\{rank\}\(D\_\{\\partial\}\)=N\_\{0\},and thereforerank⁡\(D\)=N0\\operatorname\{rank\}\(D\)=N\_\{0\}\.

LetG∂=∑γ=1rWγ​Bγ​Wγ⊤G\_\{\\partial\}=\\sum\_\{\\gamma=1\}^\{r\}W\_\{\\gamma\}B\_\{\\gamma\}W\_\{\\gamma\}^\{\\top\}be positive definite\. For anyx∈ℝN0x\\in\\mathbb\{R\}^\{N\_\{0\}\}such thatx⊤​D∂=0x^\{\\top\}D\_\{\\partial\}=0, we have

0=x⊤​\(∑γ=1rDγ\)​x=x⊤​G∂​x⟹x=00=x^\{\\top\}\\left\(\\sum\_\{\\gamma=1\}^\{r\}D\_\{\\gamma\}\\right\)x=x^\{\\top\}G\_\{\\partial\}x\\Longrightarrow x=0sinceG∂≻0G\_\{\\partial\}\\succ 0\. HenceNull⁡\(D∂⊤\)=\{0\}\\operatorname\{Null\}\(D\_\{\\partial\}^\{\\top\}\)=\\\{0\\\}, and by the rank–nullity theorem,D∂D\_\{\\partial\}has full row rankN0N\_\{0\}\.