Physics-Guided Mixture Models

• Physics-guided mixture models are probabilistic methods that incorporate physical constraints into mixture structures to maintain physically plausible multimodal representations.
• They leverage physics-based priors—via residual penalties, gating signals, and tailored component families—to filter out nonphysical solutions and enhance regime separation.
• These models have demonstrated success in diverse areas such as PDE inversion, signal processing, astrophysical clustering, and geophysical inversion.

Physics-guided mixture models are mixture-based probabilistic or conditional-computation models in which the mixture structure is constrained, regularized, or routed by physically meaningful information rather than being learned solely from statistical fit. Across recent work, this idea appears in several forms: conditional mixture density networks with distribution-level physics priors for intrinsically multimodal scientific responses, Gaussian-mixture weighting of PDE residuals in curriculum-guided PINNs, spectrum-aware mixture-of-experts for GNSS interference recognition, evidential mixture-of-experts for FLiDAR in scattering media, multivariate GMM priors in multi-physics inversion, and Gibbs-repulsive Bayesian mixtures derived from random matrix theory (Han et al., 11 Feb 2026, Yang et al., 19 May 2026, Zeng et al., 19 Jan 2026, Erbas et al., 23 May 2025, Astic et al., 2020, Cremaschi et al., 2023). The common objective is to preserve multimodality, heterogeneity, or regime structure while excluding spurious solutions that violate governing equations, physical supports, geometry, or interpretable domain structure.

1. Conceptual basis

A central motivation for physics-guided mixtures is that multimodality in scientific systems is often intrinsic rather than observational. Conditional responses can be multimodal because of latent regime switching, bifurcations, heterogeneity, or non-unique mappings; examples include cusp or pitchfork bifurcations in nonlinear ODEs, slow–fast stochastic dynamics with parameter-dependent bistability, shock-wave physics with elastic, plastic, and phase-transformation branches, and reaction–diffusion PDEs with multiple stable steady states (Han et al., 11 Feb 2026). In such settings, a single conditional mean is structurally inadequate, while an unconstrained flexible density can assign probability mass to nonphysical modes.

Mixtures are attractive because they give an explicit and likelihood-tractable representation of conditional or latent heterogeneity. In the scientific MDN setting, the components can align with physical regimes and the gating probabilities become regime probabilities (Han et al., 11 Feb 2026). In inverse problems, a multivariate GMM can encode petrophysical facies over vectors of physical properties, thereby coupling surveys that otherwise depend on different forward physics (Astic et al., 2020). In astronomy, physically derived component families such as the isothermal sphere, King profile, and Plummer sphere replace generic Gaussians when the target structure is governed by dynamical arguments rather than convenience (Kuhn et al., 2017).

Physics guidance does not have a single formal meaning. In the current literature it includes at least four distinct mechanisms. First, physics can enter as a soft penalty or residual in the training loss, as in MDNs or PINNs (Han et al., 11 Feb 2026, Yang et al., 19 May 2026). Second, it can define the routing signal in a mixture-of-experts architecture, for example through power spectral density features in GNSS or temporal segmentation of photon histograms in FLiDAR (Zeng et al., 19 Jan 2026, Erbas et al., 23 May 2025). Third, it can determine the component family itself, as in astrophysical cluster profiles or Coulomb-repulsive priors derived from log-gas ensembles (Kuhn et al., 2017, Cremaschi et al., 2023). Fourth, it can be injected through feature engineering and feasible-set projections, such as normalized energy surrogates in additive manufacturing or sector constraints and MPC smoothing in trajectory prediction (Basterrech et al., 30 Oct 2025, Li et al., 25 Jul 2025).

2. Main model families

The present landscape includes several distinct but related architectures.

Family	Physics guidance	Representative papers
Conditional MDNs	Distribution-level residuals on component means, weighted by gating probabilities	(Han et al., 11 Feb 2026)
Residual-space GMMs for PINNs	GMM fit to PDE residual distribution; curriculum weights from responsibilities and precisions	(Yang et al., 19 May 2026)
Physics-guided MoE	Routing from PSD, temporal segments, or physics-derived complexity cues	(Zeng et al., 19 Jan 2026, Erbas et al., 23 May 2025, Srivastava et al., 2022)
Bayesian repulsive mixtures	Gibbs measures with logarithmic repulsion and closed-form partition functions	(Cremaschi et al., 2023)
GMM priors in inversion	Multivariate class priors over physical-property vectors	(Astic et al., 2020)
Regime-unmixing mixtures	Pixel-wise nonlinear activation from observable physical features	(Pacheco et al., 5 May 2026)

In conditional density modeling, the canonical representation is

$$p(y|x)=\sum_{k=1}^K \pi_k(x)\, p_k(y|x;\theta_k(x)),$$

with Gaussian components used in the main scientific MDN instantiation for analytical tractability and moment formulas (Han et al., 11 Feb 2026). The same paper defines the total objective as

$$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$

so the mixture retains exact likelihood training while incorporating physical structure at the distribution level.

In mixture-of-experts form, the generic aggregation is

$y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$

with nonnegative gating weights summing to one (Zeng et al., 19 Jan 2026). In PhyG-MoE, the gate is driven by PSD features extracted from multitaper spectra, while the experts operate on STFT spectrograms (Zeng et al., 19 Jan 2026). In the quantum PG-MoE for Schrödinger’s equation, the gate is non-parametric and hard-routes basis configurations to experts according to the magnetization quantum number $S_z$, so the decomposition is dictated by the Hamiltonian’s physical structure rather than by learned soft assignments (Srivastava et al., 2022).

In Bayesian clustering and random-component mixtures, the physics guidance can reside in the prior over component locations. Coulomb-repulsive mixtures use Gibbs measures of the form

$$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$

with

$$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$

so the atoms behave like interacting particles in a log-gas (Cremaschi et al., 2023). The key computational feature is that the Hermite, Laguerre, and Jacobi ensembles have closed-form partition functions, which makes variable-$K$ Bayesian inference feasible without doubly intractable normalization.

Other variants retain the GMM formalism but redefine what the mixture acts on. In CGMPINN, a 1D GMM is fit to the scalar PDE residuals at collocation points, and the resulting responsibilities induce adaptive sample weights in a dual curriculum (Yang et al., 19 May 2026). In PGRU, the mixture is not over endmember abundances alone but over nonlinear residual mechanisms—GBM, PPNM, and Hapke—combined through learned attention and modulated by a pixel-wise regime scalar $\xi_i\in[0,1]$ predicted from physically observable features (Pacheco et al., 5 May 2026).

3. Mechanisms of physics injection

A distinctive feature of the scientific MDN formulation is that physics acts at the distribution level rather than on a single sampled prediction. Let $R(y,x)$ denote a physics residual. The general penalty is

$$L_{\mathrm{phys}} = \mathbb{E}_x\!\left[\sum_{k=1}^K \pi_k(x)\int R(y,x)\,p_k(y|x)\,dy\right],$$

while the implementation used in the paper evaluates the residual at the component means,

$$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$0

This $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$1-weighting preserves the semantics that dominant modes are enforced more strongly and inactive modes are not unduly penalized (Han et al., 11 Feb 2026). The residuals used include monotonicity for Hugoniot curves and a steady-state PDE residual for the Chafee–Infante equation (Han et al., 11 Feb 2026).

In CGMPINN, the physics signal is the empirical distribution of PDE residuals rather than the predicted state itself. A GMM is periodically fit to the residuals $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$2, responsibilities $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$3 are computed, difficulty scores are defined by responsibility-weighted residual magnitudes, and a shared curriculum parameter drives both easy-to-hard weighting and precision-based variance modulation (Yang et al., 19 May 2026). The weighted PDE loss is proved uniformly equivalent to the standard PDE loss under bounded GMM variances, and the time-varying objective admits a sublinear convergence guarantee for the gradient norm (Yang et al., 19 May 2026).

In MoE systems, physics frequently enters through the gate. PhyG-MoE uses a PSD encoder with two $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$4 Conv–Pooling blocks to route GNSS samples into three computational regimes: sparse peaks to MobileNetV4, wider structured occupancy to SK-GhostNet, and chaotic or saturated spectra to TransNeXt (Zeng et al., 19 Jan 2026). EvidenceMoE encodes the temporal physics of photon transport by assigning an Early Expert to the rising edge and peak, a Late Expert to the tail, and a Global Expert to the whole sequence; paired Evidence-Based Dirichlet Critics output Beta parameters, quality scores, and corrective residuals, and a Decider Network fuses the expert outputs through sigmoid gating (Erbas et al., 23 May 2025). PhysVarMix adds hard sector-specific boundary conditions and an MPC smoother after mixture sampling, so that diverse candidate trajectories remain within a feasible geometric sector and satisfy an Ackermann-style dynamics model during refinement (Li et al., 25 Jul 2025).

Physics can also guide the component family or support. In astronomy, an isothermal ellipsoid matches the central cusp and extended wings of a young star cluster more naturally than a multivariate normal, which may need multiple concentric Gaussian components to approximate a single physical cluster (Kuhn et al., 2017). In hyperspectral unmixing, PGRU predicts

$$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$5

where the regime residuals $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$6 come from GBM, PPNM, and Hapke, and $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$7 is obtained from spectral curvature, NDVI, NDVI gradient, EMP/DMP, LBP, and related observable cues (Pacheco et al., 5 May 2026). In additive manufacturing, the physics enters through the engineered surrogate $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$8, intended to reflect absorption or diffusion balance and to regularize multimodal marginals before GMM classification (Basterrech et al., 30 Oct 2025).

4. Representative domains and empirical behavior

In conditional multimodal scientific modeling, the MDN framework with distribution-level physics priors is evaluated on bifurcation phenomena, stochastic differential equations, atomistic shock dynamics, and reaction–diffusion steady states (Han et al., 11 Feb 2026). In the nonlinear ODE bifurcation example, with comparable capacity and equal training budget, the MDN with 457 parameters captured well-separated modes at $$L_{\mathrm{total}} = L_{\mathrm{NLL}} + \lambda_{\mathrm{phys}} L_{\mathrm{phys}} \; (+ \lambda_{\mathrm{class}} L_{\mathrm{class}}),$$9, while a conditional flow matching model with 521 parameters placed noticeable mass in intermediate regions between modes; for $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$0, the MDN better captured the broadened shape, and for unimodal $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$1, both were good (Han et al., 11 Feb 2026). In the Chafee–Infante case, adding the steady-state residual pulls component means toward the true equilibria while variances increase to cover remaining transient data mass, preserving likelihood fit while enforcing physical plausibility (Han et al., 11 Feb 2026).

In PDE solving, CGMPINN uses a residual-space GMM to guide a dual curriculum and reports the lowest relative $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$2 and maximum absolute errors among all compared methods across six benchmarks at comparable cost (Yang et al., 19 May 2026). On 1D Poisson, CGMPINN achieves $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$3, relative $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$4, and $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$5, reducing $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$6 by $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$7 relative to the standard PINN; on the heat equation, the reduction versus standard PINN is $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$8; on the damped wave equation, the runtime is the lowest among all methods (Yang et al., 19 May 2026). The ablation study shows that neither GMM weighting alone nor curriculum alone matches the full combination of GMM, curriculum, and precision modulation (Yang et al., 19 May 2026).

In signal processing and sensing, physics-guided MoE designs explicitly align model capacity with physical complexity. PhyG-MoE evaluates 21 GNSS jamming categories and reports an overall accuracy of $y(x)=\sum_{k=1}^K g_k(x)\,f_k(x;\theta_k),$9, with 11.19 M active parameters, 32.27 M total capacity, 1.71 G FLOPs, and 4.08 ms inference time; it improves overall accuracy by $S_z$0 over TFPENet while reducing FLOPs by $S_z$1 (Zeng et al., 19 Jan 2026). EvidenceMoE, trained on realistically simulated FLiDAR data generated from photon transport models in tissue, achieves $S_z$2 for depth and $S_z$3 for fluorescence lifetime, while heteroscedastic experts alone require about 500 epochs to converge versus about 70 for EvidenceMoE (Erbas et al., 23 May 2025). In hyperspectral unmixing, PGRU improves over LMM, GBM, and PPNM on Samson, Jasper Ridge, and Urban, with physical coherence $S_z$4, $S_z$5, and $S_z$6, respectively (Pacheco et al., 5 May 2026).

In autonomous systems and scientific computing, the same design principle appears in different guises. PhysVarMix combines a variational mixture head with sector projection and MPC smoothing; on Lyft it reports Collision $S_z$7, Off-road $S_z$8, Discomfort $S_z$9, and $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$0 error $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$1 m, while removing MPC smoothing drives Discomfort to $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$2 (Li et al., 25 Jul 2025). The quantum PG-MoE for high-dimensional eigen-solvers is 150x smaller than the baseline network and remains competitive in generalization by decomposing the wavefunction coefficients into experts indexed by $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$3 sectors and adding a variational-energy loss based on the Hamiltonian (Srivastava et al., 2022). In geophysical inversion, a dynamic multivariate GMM prior over density contrast and magnetic susceptibility differentiates two kimberlite facies in a joint gravity–magnetics inversion that single-physics inversions cannot separate (Astic et al., 2020). In model discrepancy analysis for a concrete bridge, a parameter-cluster mixture identifies systematic top-sensor deviations when irradiation is omitted, and the ensuing humidity-extended boundary model reduces the reported discrepancy from $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$4 to $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$5 (Villani et al., 16 Mar 2026).

5. Inference, optimization, and evaluation

The estimation procedures span both classical latent-variable inference and end-to-end gradient optimization. GMM-based models use the standard EM pattern: responsibilities

$$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$6

are computed in the E-step, followed by updates of weights, means, and covariances in the M-step (Kuhn et al., 2017, Basterrech et al., 30 Oct 2025). In the civil-engineering mixture for model improvement, EM operates on a Bayesian marginal posterior over assignment variables, parameter vectors, and Dirichlet-distributed mixture weights, with soft assignment probabilities

$$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$7

and a numerical M-step because the forward model is a simulator rather than a Gaussian kernel (Villani et al., 16 Mar 2026). In Coulomb-repulsive Bayesian mixtures, MCMC is feasible because the partition functions for the Hermite, Laguerre, and Jacobi ensembles are available in closed form, which stabilizes birth–death moves and updates of the repulsion parameter $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$8 (Cremaschi et al., 2023).

Neural conditional mixtures and MoE systems are usually optimized by stochastic gradient methods. The scientific MDN case studies use small fully connected MLPs with ELU activations, softmax outputs for $$p(\theta_1,\dots,\theta_K)\propto \exp(-\beta H(\theta)),$$9, positive parameterizations for $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$0, and ADAM for 20,000–50,000 iterations at learning rate $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$1 (Han et al., 11 Feb 2026). CGMPINN uses Adam for coarse exploration followed by L-BFGS for fast local convergence, with the GMM refit every $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$2 iterations and optional ReLoBRaLo balancing across PDE, boundary, and initial-condition terms (Yang et al., 19 May 2026). PhyG-MoE uses AdamW with learning rate $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$3, weight decay $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$4, a 10-epoch warm-up, cosine annealing to epoch 50, and a load-balancing auxiliary loss to avoid expert collapse (Zeng et al., 19 Jan 2026).

Evaluation is likewise domain-specific but reveals recurring themes. Conditional-density models emphasize validation NLL, calibration curves, CRPS, PIT histograms, coverage, and physics violation scores (Han et al., 11 Feb 2026). PINN variants report relative $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$5, maximum absolute errors, and training cost, together with theoretical guarantees such as uniform equivalence and explicit weighting-induced bias terms (Yang et al., 19 May 2026). Sensing applications use OA, F1, NRMSE, and runtime or FLOPs (Zeng et al., 19 Jan 2026, Erbas et al., 23 May 2025). In trajectory prediction, collision rate, off-road rate, discomfort, and $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$6 error play the role of physically meaningful downstream metrics (Li et al., 25 Jul 2025). In discrepancy analysis, KS-based deficiency and model discrepancy quantify whether the simulator remains statistically compatible with the measurements under physically meaningful parameter clusters (Villani et al., 16 Mar 2026). A plausible implication is that physics-guided mixtures are best assessed by paired statistical and physical criteria rather than by likelihood alone.

6. Limitations and directions

Several limitations recur across the literature. If the governing equations or constraints are misspecified, the physics prior can bias the model away from the data (Han et al., 11 Feb 2026). Weighting parameters such as $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$7, curriculum saturation, or evidential regularizers can be highly consequential: over-regularization distorts the density or the PDE objective, under-regularization weakens the physical benefits, and inappropriate evidence sensitivity can degrade performance, as seen when $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$8 worsens EvidenceMoE depth estimation (Han et al., 11 Feb 2026, Yang et al., 19 May 2026, Erbas et al., 23 May 2025). Mixture identifiability also remains problematic: label switching, mode collapse, and sensitivity to $$H(\theta)=\sum_i U(\theta_i)-\sum_{i<j}\log|\theta_i-\theta_j|,$$9 appear in MDNs, Bayesian mixtures, and residual-space GMMs alike (Han et al., 11 Feb 2026, Cremaschi et al., 2023, Yang et al., 19 May 2026).

Model-class restrictions are another persistent concern. Many implementations still use Gaussian components for tractability, even when supports are bounded, tails are asymmetric, or within-component structure is itself multimodal (Han et al., 11 Feb 2026, Basterrech et al., 30 Oct 2025). Router interpretability can also be partial rather than explicit: PhyG-MoE’s gating learns thresholds implicitly from PSD topology rather than computing an interpretable scalar entropy, and EvidenceMoE’s sigmoid gating does not enforce sum-to-one normalization (Zeng et al., 19 Jan 2026, Erbas et al., 23 May 2025). Some methods remain computationally heavier than simpler baselines because they evaluate multiple experts, multiple physical regimes, or repeated simulator calls in EM (Erbas et al., 23 May 2025, Pacheco et al., 5 May 2026, Villani et al., 16 Mar 2026). Epistemic uncertainty is often absent: the physics-guided MDN for shock dynamics captures aleatoric and multimodal variability but does not explicitly capture epistemic uncertainty in data-sparse regions (Han et al., 11 Feb 2026).

The proposed extensions are correspondingly diverse. The MDN framework explicitly suggests richer component families, per-component normalizing flows or copulas, skewed or truncated distributions, Dirichlet-process mixtures, variational Bayesian MDNs, hierarchical regime models, adaptive $K$0, sparsity penalties, and weak-form or uncertain physics operators learned jointly with the mixture (Han et al., 11 Feb 2026). The additive-manufacturing study points toward multi-resolution hierarchical models, frequency-domain parameter optimization, and multimodal observations such as microtomography and acoustic emission (Basterrech et al., 30 Oct 2025). The geophysical inversion framework suggests spatially aware mixture models, nonlinear rock-physics maps, and deep generative priors over facies-consistent property fields (Astic et al., 2020). PGRU states that formal ablations and blind-unmixing extensions with endmember variability remain future work (Pacheco et al., 5 May 2026). This suggests that the field is converging on a broader principle: mixtures are most useful when physical structure is not merely appended as a penalty, but is used to define what the components mean, where they apply, and how they are evaluated.