Multi-output Gaussian process prediction of physical fields under linear equality constraints

Mahamat Hamdan Nassouradine; Clément Gauchy; Pierre-Emmanuel Angeli; Sébastien da Veiga

Source-linked AI summary

Multi-output Gaussian process prediction of physical fields under linear equality constraints

Mahamat Hamdan Nassouradine, Clément Gauchy, Pierre-Emmanuel Angeli, Sébastien da Veiga

arXiv:2608.25709v1Aug 26, 2026stat.MLcs.LG

TL;DR

Predicting multiple high-dimensional physical fields under linear equality constraints challenges Gaussian-process models because outputs must be handled jointly while preserving physics. The paper combines constraint-preserving row-wise PCA with analytically constrained multi-output GP regression, achieving competitive or superior performance, especially in data-scarce regimes.

Problem
Existing deductive GP strategies arbitrarily select one output to recover, making predictive accuracy and uncertainty quantification sensitive to that choice.
Method
The framework combines row-wise PCA, which preserves linear equality constraints in latent space, with a constrained multi-output GP whose kernel encodes them analytically.
Results
Across population dynamics and industrial CFD applications, the method achieved competitive or superior predictive performance, with its advantage most pronounced in data-scarce regimes.
Takeaways & Limitations
Joint modeling of all outputs with analytical constraint encoding provides a robust alternative to deductive approaches, particularly when training data are scarce.
Takeaways & Limitations
Optimal redistribution strategies for handling more general constraints remain an open direction for future work.

Abstract

from arXiv · show

We address the simultaneous prediction of multiple high-dimensional physical fields governed by linear equality constraints, a setting that arises in many real-world applications in physics machine learning. Gaussian process (GP) regression is a widely used surrogate modeling approach due to its effectiveness in small-sample regimes and its ability to provide uncertainty quantification. However, applying GP models in this setting raises two major challenges: the high dimensionality of the discretized output fields and the enforcement of the physical constraint in predictions. For the latter, a common strategy consists in deducing one output from the others via the constraint relation. Through a benchmark, we show that this deductive approach is sensitive to the arbitrary choice of which output to deduce, affecting both predictive accuracy and uncertainty quantification. Consequently, there is a need for an approach that treats all fields symmetrically while strictly respecting the underlying physics. Motivated by these limitations, we propose a robust framework for jointly modeling constrained multi-field data. Our approach first leverages a specific PCA procedure for multi-field data, coined row-wise PCA, which has the interesting property of preserving the constraint in the latent space. Since standard PCA strategies for multi-field data do not preserve such constraints, we investigate theoretically the optimality of the row-wise choice. In a second step, we consider a linearly-constrained multi-output GP approach based on a specific kernel parametrization which is trained on the latent space of row-wise PCA. The proposed framework is validated on a population dynamics problem and on an industrial CFD application, which involves the prediction of Reynolds stress tensor components under the incompressibility constraint.

1. Introdution

The paper introduces Row-CMO, a unified framework that reduces multi-field dimensionality while preserving linear equality constraints for constrained Gaussian-process prediction. It addresses deductive approaches’ sensitivity to the arbitrarily deduced field through row-wise PCA and constrained multi-output GP regression.

Motivation: Linear equality constraints matter because they ensure physically consistent predictions and are expected to improve metamodel accuracy in applications such as incompressible flows and conservation laws.Examples include velocity incompressibility and species mass fractions summing to unity.
Contribution: Row-CMO predicts multi-field outputs under linear equality constraints by preserving the constraint structure during dimensionality reduction and applying constrained GP regression in latent space.The framework uses PCA because nonlinear alternatives are less robust in small-sample regimes and incompatible with linear constraints.
Motivation: Benchmarking shows deductive approaches are heterogeneous, with performance potentially degrading significantly depending on which output is arbitrarily deduced.The deductive approach fits an unconstrained GP to the remaining Q−1 fields and derives the final output algebraically.
Contribution: The paper theoretically compares field-wise, column-wise, and row-wise multi-field PCA, attributing global methods’ excess reconstruction error solely to inter-field spectral heterogeneity.Row-wise PCA is selected for its constraint-preserving latent representation.
Validation: Row-CMO is validated on Lotka–Volterra population dynamics with a sum-to-zero constraint and CFD prediction of Reynolds stress components under incompressibility.The CFD case concerns the Buice–Eaton 2D diffuser.

2. Problem statement

The problem is to learn Q coupled, high-dimensional scalar fields from few noise-free simulator runs while enforcing known linear equality constraints exactly. Gaussian processes provide small-sample effectiveness and uncertainty quantification, while parametrization constrains every prior and posterior realization without privileging a deduced field.

Problem formulation: The simulator maps D-dimensional inputs to Q coupled scalar fields evaluated at S fixed coordinates, with only tens to hundreds of expensive, noise-free training runs.The objective is to approximate the resulting vector-valued mapping across the input domain.
Problem formulation: The fields satisfy a known linear equality constraint with nonzero coefficients, initially specialized to constant coefficients and a zero right-hand side.Input-dependent coefficients, nonzero constraint functions, and noisy observations are deferred to Section 4.6.
Gaussian-process modeling: Gaussian process regression is adopted because it performs effectively in small-sample regimes and provides principled uncertainty quantification.The framework uses both single-output and multi-output GP formulations, including matrix-valued kernels such as the LCM.
Constrained multi-output GPs: Constraint parametrization expresses the constrained function as a linear transformation of an unconstrained function, transferring the restriction to a matrix whose columns lie in the constraint null space.The transformation induces a valid GP because Gaussian processes are stable under linear maps.
Constrained multi-output GPs: Every prior or posterior realization satisfies the linear constraint exactly at every input, unlike an unconstrained ICM whose learned covariance cannot guarantee exact satisfaction.The non-unique parametrization is optimized jointly with kernel hyperparameters, making the construction stable to the choice of constrained entry while selecting the best-fitting representation.

3. Linear equality constraints and deductive approaches

The deductive strategy enforces linear equality constraints by selecting one output to compute from the others, but benchmark results show that this arbitrary choice affects predictive accuracy and favors jointly constrained modeling. In the three-output experiment, CMOGP provides the most stable and strongest overall benchmark performance, motivating constraint-preserving dimensionality reduction for high-dimensional fields.

Deductive approach: Deduction enforces the predictive-mean constraint by construction and can be combined with independent GPs or LCM-based MOGPs without changing their training procedures.Posterior samples can also be generated to comply with the constraint.
Experimental setup: The benchmark uses three scalar outputs satisfying 𝑓1+𝑓2+𝑓3=0, 200 replications, a fixed 200-point test set, and training sizes 𝑁∈{20, 50, 100}.The experiment deliberately combines Ishigami and Branin functions with different complexity and smoothness, while evaluating uncertainty through average 90% predictive-interval length.
RMSE across deduction scenarios: Independent GP accuracy worsens when the evaluated output is the deduced one, and this gap persists as the training set grows.Figure 1 compares RMSE distributions for 𝑓2 across 200 replications, training sizes 𝑁∈{20, 50, 100}, and all deduction choices.
Win rate analysis: CMOGP is the best-performing joint approach in the win-rate benchmark, while LCM is more stable than independent GP across deduction scenarios.CMOGP models all outputs jointly with a constrained kernel, whereas LCM still computes one output during post-processing.
Motivation: Because deductive approaches are unstable and constrained MOGP regression does not scale directly to large output dimension, the study motivates constraint-preserving dimensionality reduction.The stated large-output setting has dimension 𝑃=𝑄×𝑆.

4. Principal component analysis for multi-field data

The section compares three PCA strategies for reducing multi-field outputs and shows that row-wise PCA uniquely preserves linear constraints in latent, predicted, and reconstructed fields. Its reconstruction-error bounds quantify when multi-field methods avoid a penalty and provide a data-driven comparison of row-wise and column-wise PCA.

Motivation and scope: Theoretical guarantees motivate row-wise PCA for physical fields arising from the same PDE system, where similar spatial covariance structures are common.The approach is designed to reduce output dimensionality while preserving the constraint structure required by the constrained MOGP.
PCA strategies: Field-wise PCA minimizes each field’s reconstruction error independently but uses separate bases, whereas column-wise PCA captures inter-field correlations and yields scalar latent variables modeled by independent GPs.Row-wise PCA instead uses a shared spatial basis and produces vector-valued latent weights whose components can be modeled jointly by a MOGP.
Constraint preservation: Row-wise PCA is the only strategy that transfers the linear constraint exactly to latent weights and reconstructed fields, regardless of truncation rank.Consequently, it supports the constrained MOGP while column-wise and field-wise PCA do not generally preserve the constraint.
Row-wise versus column-wise PCA: Neither row-wise nor column-wise PCA dominates globally: Theorem 3 gives an exact criterion for comparing their global reconstruction errors on all fields.The bounds and criterion can be evaluated directly from data to identify similar-field, heterogeneous-field, or row-wise-versus-column-wise regimes.
Reconstruction-error bounds: Theorem 1 and Theorem 2 show that each multi-field method’s excess reconstruction error is governed by inter-field spectral heterogeneity relative to the field-wise optimum.When fields have similar spatial or observation structures, the excess error tends to zero; shared spatial bases make row-wise excess error vanish, while shared observation bases make column-wise excess error vanish.

5. Numerical experiments

Numerical experiments validate Row-CMO on Lotka–Volterra population dynamics and industrial CFD, showing predictive advantages over deductive alternatives and persistent sensitivity to deduction-output selection. The results also support the row-wise PCA diagnostic and demonstrate accurate constrained field prediction in a realistic high-dimensional application.

A priori diagnostic and predictive accuracy: Theorem 3’s dominance criterion provides a training-data diagnostic: column-wise PCA captures only marginally more variance initially, with the difference approaching zero beyond m≈4.In the CFD case, the criterion is positive but small, on the order of 10^-3, and monotonically decreasing.
Lotka–Volterra experiments: Row-CMO consistently achieves the lowest RRMSE at higher latent dimensions across all three Lotka–Volterra training-set sizes, outperforming deductive alternatives.The advantage appears despite near-equivalent training-data variance explained by row-wise and column-wise PCA for m≥5.
Lotka–Volterra experiments: The Lotka–Volterra results attribute Row-CMO’s generalization advantage to information pooling across fields and analytical constraint encoding in the constrained MOGP.The shared basis uses stacked observations, while the constrained model encodes the physical relation directly; the advantage is most pronounced in data-scarce regimes.
Lotka–Volterra experiments: Deductive methods show persistent performance heterogeneity across deduction scenarios for every training size and latent dimension, confirming sensitivity to the arbitrary deduced output.This variability is visible across Col-Indep, Fw-Indep, and Fw-LCM families.
Spatial distribution of predictions and errors: In industrial CFD, Row-CMO accurately reproduces the spatial structure, intensity, and extent of recirculation in the predicted Reynolds-stress-related fields, with highest uncertainty in the recirculation region.The representative visualization uses M=8 latent dimensions and compares predictions with simulations for τ11, τ22, and k.

6. Conclusion

The paper proposes a framework combining row-wise PCA with constrained multi-output GP regression for simultaneous prediction of high-dimensional physical fields under linear equality constraints. It achieves competitive or superior performance in population dynamics and industrial CFD applications, especially in data-scarce regimes, while motivating extensions to irregular observations and spatio-temporal data.

Conclusion: The proposed framework combines constraint-preserving row-wise PCA with constrained multi-output GP regression for simultaneous prediction of high-dimensional physical fields.It was validated on the Lotka–Volterra population dynamics system and an industrial CFD application predicting Reynolds stress tensor components under incompressibility.
Conclusion: The method achieved competitive or superior predictive performance against deductive alternatives, with its generalization advantage most pronounced in data-scarce regimes.Improved performance was attributed to information pooling through a shared spatial basis and analytical constraint encoding.
Conclusion: Future work includes handling fields observed at different input locations by exploiting constraint-based coupling and extending the methodology to spatio-temporal physical simulations.The spatio-temporal extension is suggested following Mak et al. (2018).

A. Proofs of the theoretical results

The proofs rely on the Davis–Kahan sin Θ theorem to control subspace perturbations and on the relation between projector distance and principal angles.

Theoretical tools: The proofs repeatedly invoke the Davis–Kahan sin Θ theorem.The theorem is stated for symmetric matrices using their leading eigenvectors and an eigenvalue-separation condition.
Theoretical tools: The argument also connects distances between rank-m projectors to the principal angles between their underlying column spaces.This relation translates subspace-angle bounds into projector-distance bounds.

A.1. Proof of Theorem 1 (Row-wise excess error bound)

The proof derives the row-wise excess error bound using trace formulations, projector identities, trace Cauchy–Schwarz, and a perturbation argument applying Theorem 4 to Σ_k and Σ_row.

The proof begins with a trace formulation and uses projector idempotence to establish the error expression.
Trace Cauchy–Schwarz and a perturbation bound from Theorem 4 are then applied with Σ_A = Σ_k and Σ_B = Σ_row.
Finally, the proof bounds projector distance using (39) and substitutes (42) into (41) to obtain the stated result.

A.2. Proof of Theorem 2 (Column-wise excess error bound)

The proof of Theorem 2 follows the argument for Theorem 1 in observation space. It uses SVD projection equivalence, applies Theorem 4 with the column-wise and target kernels, and substitutes the resulting perturbation expression into the excess-error bound.

Proof strategy: The proof operates in observation space and uses the SVD identity Y_kP_k,m = Π_k,mY_k to express equivalent left and right projections.This establishes the form of the column-wise excess error before applying the perturbation bound.
Kernel perturbation: Theorem 4 is applied with Σ_A = K_k and Σ_B = K_col, yielding the perturbation ΔK_k = K_col − K_k.
Bound derivation: Substituting the perturbation expression into the preceding excess-error relation produces the stated bound.

A.3. Proof of Theorem 3 (Row-wise vs. column-wise) · B. Uncertainty quantification of the standard 𝑘–ε turbulence model parameters

Theorem 3’s proof compares row-wise and column-wise PCA through their decomposed global errors and shared Frobenius norms. The standard k–ε model is then characterized by its turbulence assumptions and five calibrated parameters.

A.3. Proof of Theorem 3 (Row-wise vs. column-wise): A.3 shows that both row-wise and column-wise global errors decompose into total energy minus explained variance.
A.3. Proof of Theorem 3 (Row-wise vs. column-wise): Because Y_row and Y_col contain identical elements in different arrangements, their Frobenius norms coincide.
A.3. Proof of Theorem 3 (Row-wise vs. column-wise): Under the corollary’s shared-eigenvector condition, rank-m truncation yields identical row-wise and column-wise projection matrices.
B. Uncertainty quantification of the standard 𝑘–ε turbulence model parameters: The standard k–ε turbulence model is widely used in CFD because it provides reliable predictions across many flows at reasonable computational cost.
B. Uncertainty quantification of the standard 𝑘–ε turbulence model parameters: The model treats turbulence as an additional viscous effect through a spatially varying eddy viscosity coefficient ν_t.
B. Uncertainty quantification of the standard 𝑘–ε turbulence model parameters: Dimensional analysis expresses eddy viscosity using turbulent kinetic energy k, dissipation rate ε, and calibrated constant C_µ.
B. Uncertainty quantification of the standard 𝑘–ε turbulence model parameters: The turbulent kinetic-energy equation follows from Reynolds decomposition, whereas the dissipation-rate equation requires a simplified modeled form.
B. Uncertainty quantification of the standard 𝑘–ε turbulence model parameters: The model ultimately contains five experimentally or semi-analytically calibrated parameters, whose standard values are reported in Table 1.

B.1. Calibration of 𝐶µ · B.3. Prandtl numbers σ𝑘and σε · B.5. Sampling procedure for the standard 𝑘−ε turbulence model parameters

The standard k–ε parameters are calibrated and sampled from experimentally informed uncertainty distributions while enforcing physical and model-consistency constraints. The procedure combines Reynolds-stress measurements, turbulence-decay information, uncertain Prandtl numbers, and a constrained joint sampling algorithm.

B.1. Calibration of 𝐶µ: At y/H = 4, local equilibrium yields a nominal Cµ = 0.06 from Reynolds-stress measurements, close to the classical value.The calibration uses the parallel-shear-flow simplification and P_k = ε.
B.1. Calibration of 𝐶µ: Cµ uncertainty is estimated by drawing 1,000 samples from Gaussian Reynolds-stress distributions based on Buice and Eaton (2000b) measurements and uncertainties.The component uncertainties are 4% for τ11, 5% for τ22, and 10% for τ12; τ33 is assigned an ad-hoc uncertainty of 4%.
B.1. Calibration of 𝐶µ: Cε2 is calibrated from isotropic homogeneous decay, assuming k(t) ∝ t^−β and inferring β from 650 experimental and DNS values reported by Meldi and Sagaut (2012).A uniform β distribution is used to generate a 1,000-sample Cε2 distribution shown with a histogram and Gaussian KDE.
B.3. Prandtl numbers σ𝑘and σε: σk and σε receive 10% uncertainty based on Foken (2006), with uniform distributions truncated to enforce 2σk − σε ≤ 1.Only i.i.d. samples satisfying the interface-behavior constraint are retained.
B.3. Prandtl numbers σ𝑘and σε: The Von Kármán constant κ is modeled by fitting a Gaussian to n = 8 values from Foken (2006), while local equilibrium relates Cε1 and Cε2.The supplied passage introduces the Cε1–Cε2 relation under boundary-layer flow and ε = Pk.

C. Complete win rate analysis of the deductive approach benchmark · C.1. Predictive intervals

The complete benchmark shows that deductive-method win rates vary across deduction scenarios and training-set sizes, whereas CMOGP remains governed by its fixed RMSE. Predictive intervals have comparable coverage across methods, while CMOGP provides stable, shorter intervals across scenarios.

C. Complete win rate analysis of the deductive approach benchmark: Deductive-method win rates change across deduction scenarios at every training-set size, unlike CMOGP’s mechanically determined variation from its fixed RMSE.This scenario dependence is present for N=20, 50, and 100.
C. Complete win rate analysis of the deductive approach benchmark: Across 200 replications, Figure 11 compares each method’s RRMSE win rate for f1, f2, and f3 at N∈{20, 50, 100} across all deduction scenarios.The highlighted column identifies the scenario in which the evaluated output is deduced, while CMOGP is unchanged across scenarios.
C. Complete win rate analysis of the deductive approach benchmark: CMOGP is evaluated as the same model in every deduction scenario, providing a scenario-invariant reference for the benchmark comparison.The figure aggregates results over 200 independent replications.
C.1. Predictive intervals: Figure 12 reports the distribution of average 90% predictive-interval lengths for f2 across N∈{20, 50, 100} and the three deduction scenarios.The comparison covers all deduction scenarios and training-set sizes in the benchmark.
C.1. Predictive intervals: All methods achieve comparable empirical coverage rates, so shorter predictive intervals indicate more informative predictions.Interval length is therefore used to distinguish informativeness when coverage is comparable.
C.1. Predictive intervals: CMOGP maintains stable and shorter predictive intervals across all deduction scenarios.For deductive methods, intervals for the deduced output are obtained by propagating modeled-output posterior samples through (18).

D. Complete RRMSE results for the Lotka–Volterra experiment

Across all four Lotka–Volterra outputs and three training-set sizes, Row-CMO overtakes Col-Indep and Fw-Indep at moderate-to-high latent dimensions, while trailing at small dimensions. Its advantage widens as training data decrease.

Complete RRMSE results: Row-CMO surpasses Col-Indep and Fw-Indep at moderate-to-high latent dimensions across all four outputs and training-set sizes.At small latent dimensions, Row-CMO starts above the alternatives.
Complete RRMSE results: The performance gap between Row-CMO and the alternatives widens as the training-set size decreases from N=30 to N=10.
Complete RRMSE results: Figure 13 reports RRMSE versus latent dimension for outputs p, q, r, and s, comparing deduction scenarios k*∈{0, 1, 2, 3}; error bars show ±1 standard deviation over 10 seeds.