PaperlayerView exact revision on arXiv

Source-linked AI summary

Energetic Variational Deep Learning for Solving High-Dimensional Stationary Fokker-Planck Equations

Fanze Kong, Zhouyuan Qin, Shuangquan Xie

arXiv:2610.03325v1math.NAcs.NA

TL;DR

Stationary Fokker–Planck equations are challenging when densities lack closed forms, particularly in high dimensions. The paper combines Monte Carlo density data with variational physics losses and staged optimization. The combined method gives the lowest MAPE in all six reported examples and improves on the strong-form formulation in reported noise tests.

  • Problem

    Stationary distributions for non-gradient drift or state-dependent diffusion typically lack closed forms, and classical PDE methods face the curse of dimensionality in high-dimensional problems.

    Source: §1 · p. 2

  • Method

    The framework combines Monte Carlo density data with weak variational physics losses and uses Monte Carlo pretraining, joint loss optimization, and L-BFGS fine-tuning.

    Source: §1 · p. 2

  • Results

    Weak/Integral+MC+L-BFGS gives the lowest MAPE in all six examples, while the weak/variational hybrid consistently improves upon the strong-form formulation in reported noise tests.

    Source: §4.7 · p. 26, §5 · p. 28

  • Takeaways & Limitations

    Using a family of test functions, the method also treats stationary problems with non-gradient drift that does not satisfy pointwise orthogonality.

    Source: §5 · p. 28

    Extension of the method to time-dependent problems remains for future study.

    Source: §5 · p. 29

Abstract

from arXiv · show

This paper develops a deep learning framework for solving high-dimensional stationary Fokker-Planck equations by combining Monte Carlo weak supervision, energy-variational physics constraints, and quasi-Newton refinement. Monte Carlo statistics anchor the scale of the stationary density, while a free-energy variational loss enforces equilibrium physics; the formulation also applies to stationary systems with pointwise orthogonal rotational components and, using a family of test functions, to systems with more general non-gradient drift. Numerical results demonstrate accurate and robust density approximations across a range of high-dimensional stationary Fokker-Planck problems, including rough-potential and noisy-data settings.

1 Introduction

High-dimensional stationary Fokker–Planck equations are difficult to solve numerically, especially when invariant densities lack closed forms. The paper proposes a Monte Carlo and variational deep-learning framework, with a scale-invariant loss and extensions beyond equilibrium systems.

Motivation

  • For general stochastic dynamics, invariant densities often lack closed forms, while classical PDE methods face the curse of dimensionality in high-dimensional problems.Strong-form residuals can also bring optimization difficulties from higher-order derivatives and sensitivity to noisy data, motivating weak and variational formulations.Source: §1 · pp. 1–2

Proposed framework

  • A progressive training strategy uses Monte Carlo pretraining, joint optimization of Monte Carlo and variational losses, and L-BFGS fine-tuning.The variational loss is reformulated in degree-zero homogeneous form to address scale imbalance during joint optimization.Source: §1 · p. 2

Variational loss

  • Dividing the variational loss by total mass makes it scale-invariant and degree-zero homogeneous, with the paper presenting this loss as more robust to rough potentials and noisy data.The formulation is motivated by the free-energy dissipation structure of equilibrium steady states.Source: §1 · pp. 2–3

Proposed framework

  • The framework combines Monte Carlo data with an integral variational physics loss to learn high-dimensional steady-state densities with fewer training points and lower computational cost.The weak-form loss is derived from free-energy dissipation for equilibrium steady states.Source: §1 · p. 3

Scope

  • The framework covers equilibrium steady states and nonequilibrium steady states with pointwise-orthogonal potential and rotational components.Source: §1 · p. 3

2 Mathematical Formulation

The formulation derives physics losses from stationary Fokker–Planck structure for equilibrium and two nonequilibrium cases, then discretizes the losses and Monte Carlo data for training. Scale-invariant physics losses complement Monte Carlo data, which anchor the density scale.

2 Mathematical Formulation

  • The formulation distinguishes equilibrium zero-current states, pointwise-orthogonal nonequilibrium states, and general nonequilibrium states treated with test functions.For the pointwise-orthogonal case, the drift is decomposed into a potential gradient and a rotational component orthogonal to the potential gradient.Source: §2 · pp. 3–4
  • For equilibrium, zero flux and vanishing dissipation yield a variational physics loss for a neural-network density approximation on a truncated domain.The loss uses density normalization and a regularized density near zero to avoid numerical instability.Source: §2 · p. 4
  • For pointwise-orthogonal nonequilibrium states, the dissipation identity supplies a variational physics loss using the orthogonality and divergence-free conditions on the rotational component.Source: §2 · pp. 4–5
  • For general nonequilibrium states, a family of test functions measures weak stationary-equation residuals while allowing nonzero current and requiring only the original drift, density, and density gradient.The loss does not require the pointwise orthogonality condition used in the preceding case.Source: §2 · p. 5

2.1 Discretization

  • The Monte Carlo reference data come from Euler–Maruyama trajectories after burn-in, sampled at fixed intervals and projected onto histogram bins.The paper defines a pointwise MSE for the Monte Carlo data and uses bin centers and sizes to describe the reference bins.Source: §2.1 · pp. 5–6
  • Physics-loss integrals are discretized with Sobol quasi-Monte Carlo points in low dimensions and importance sampling in the high-dimensional example.The high-dimensional physics points are sampled from a proposal density with corresponding importance weights.Source: §2.1 · pp. 6–7
  • The physics losses are invariant to positive rescaling of the density, so their value does not depend on its global amplitude.This invariance does not enforce probability normalization; Monte Carlo pretraining determines the density scale. The paper also reports error estimates and convergence of quadrature approximations for the minimization problems in Cases 1 and 2.Source: §2.1 · p. 7

3 Learning Methods

The section presents a three-stage procedure: Monte Carlo data pretraining, joint data-and-physics optimization, and optional L-BFGS refinement. It also describes Monte Carlo reference construction and high-dimensional physics sampling.

3.1 Neural Network Architecture and Monte Carlo Setup

  • Monte Carlo reference data pair histogram densities with bins selected by density-weighted and uniform sampling; high-dimensional models use clipped log-density outputs mapped back to positive densities.Source: §3.1 · pp. 7–8

3 Learning Methods

3.2 Stage A: Monte Carlo Data Pretraining

  • Stage A fits scaled Monte Carlo data to initialize density scale and dominant spatial structure before physics constraints, retaining the lowest data-fitting-loss checkpoint.Source: §3.2 · pp. 8–9

3.3 Stage B: Joint Optimization: Variational Physics and MC Data Losses

  • Stage B alternates physics and Monte Carlo-data updates with separate Adam optimizers, fixed initial loss scales, and an independent monitoring set for checkpoint selection.Source: §3.3 · p. 9

3.4 Stage C: L-BFGS Fine-Tuning

  • Stage C fine-tunes L-BFGS variants from Stage B using the same physics sets and a high-dimensional data-fit safeguard; Gaussian-uniform importance sampling corrects spatial integrals.Source: §3.4 · pp. 10–11

4 Numerical Experiments

The numerical experiments compare seven methods using effective-region MAPE, with analytical stationary densities reserved for final evaluation, and assess weak-form training with Monte Carlo data and L-BFGS refinement. They also examine noisy-data conditions and the sensitivity of weak and strong physics losses to density perturbations.

4.1 Double-Well Potential in Two Dimensions

  • 0.3452% is the lowest MAPE for Weak+MC+L-BFGS, versus 1.8557% for Strong+MC+L-BFGS; weak-form error falls from 7.2102% (Data-only) to 1.2187% after Stage B.The final error is also below the original MC histogram’s 9.8451%.Source: §4.1 · p. 13, Figure 5 · p. 14

4.2 Rough Double-Well Potential with Orthogonal Rotation in Two Dimensions

  • 0.4813% is the lowest MAPE for Weak+MC+L-BFGS, versus 3.1981% for Strong+MC+L-BFGS; weak-form error falls from 11.3059% (Data-only) to 2.4990% after Stage B.The final value is below the MC histogram’s 9.6156% MAPE; this non-gradient test has nonzero flux and assesses the pointwise-orthogonal extension.Source: §4.2 · pp. 15–16, Figure 9 · p. 17

4.3 Four-Dimensional Nonequilibrium Steady State with a Test-Function Family

  • 4.2620% is the MAPE for WeakFamily+MC+L-BFGS, versus 7.8237% for Strong+MC+L-BFGS and 10.3307% for Data-only NN.The fixed test-function family imposes stationary-current balance without requiring R · ∇V = 0, while retaining the MC anchor and three-stage training.Source: §4.3 · pp. 17–18, Figure 13 · p. 19

4.4 Rough Quadruple-Well Potential in Four Dimensions

  • 1.7430% is the lowest MAPE for Weak+MC+L-BFGS, versus 6.3159% for Strong+MC+L-BFGS; weak-form error falls from 8.2217% (Data-only) to 3.5404% after Stage B.The final value is also below the MC histogram’s 9.6676% MAPE.Source: §4.4 · p. 20, Figure 17 · p. 21

4.5 Rough Single-Well Potential in Eight Dimensions

  • 1.2406% is the lowest MAPE for Weak+MC+L-BFGS, versus 4.2581% for Strong+MC+L-BFGS; weak-form error falls from 19.6102% (Data-only) to 10.7225% after Stage B.The eight-dimensional evaluation is reported on the (x1, x2) slice with remaining coordinates fixed at zero.Source: §4.5 · p. 22, Figure 21 · p. 23

4.6 Two-Dimensional Equilibrium Positive-Feedback System

  • 0.7518% is the lowest MAPE for Weak+MC+L-BFGS, versus 3.2612% for Strong+MC+L-BFGS; weak-form error falls from 5.3825% (Data-only) to 3.4041% after Stage B.The result supports the framework’s effectiveness for the equilibrium gene-regulatory model.Source: §4.6 · p. 25, Figure 25 · p. 26

4.7 Cross-Dimensional Summary and Ablation Studies

  • Weak/Integral+MC+L-BFGS gives the lowest MAPE in all six benchmarks, and weak-form training improves on Data-only before L-BFGS further reduces error.The stage-ablation pattern holds across all six examples.Source: Table 1 · p. 26, §4.7 · p. 26

4.8 Robustness to Noisy MC Data

  • Weak+MC+L-BFGS has lower MAPE than Strong+MC+L-BFGS in all nine two-dimensional noisy-data tests; at α = 0.2, it gives 2.5418% and 3.2789% in four and eight dimensions.The corresponding strong-form errors at α = 0.2 are 12.1804% and 35.9194%; WeakFamily+MC+L-BFGS also has lower MAPE at all four tested roughness scales.Source: Table 2 · p. 27, §4.8 · p. 27, Table 3 · p. 27
  • The experiment perturbs uθ only for physics-loss evaluation, holding coordinates, drift, MC labels, and evaluation points fixed; weak loss uses first derivatives while strong loss uses second derivatives.The same normalization is retained across noise levels to compare physics losses on a common scale.Source: §4.8 · pp. 27–28
  • At α = 0.025, weak MAPE is 4.155% and strong MAPE is 27.973%; weak changes by less than 0.30% across noise levels, versus 16.793% for strong.Without perturbations, weak MAPE is 4.201%.Source: §4.8 · p. 28
  • Weak+MC+L-BFGS achieves lower MAPE in all 15 tested noise conditions, suggesting greater robustness to noise in the MC data.Source: §4.8 · p. 28

5 Conclusion

The framework combines Monte Carlo density data with a variational physics constraint, extending from equilibrium to several classes of non-equilibrium stationary systems. Numerical experiments report effectiveness across dimensions, rough potentials, and noisy observations, while time-dependent problems remain future work.

5 Conclusion

  • The three-stage framework uses Monte Carlo data for weak supervision and a variational loss to enforce physical consistency.For equilibrium steady states, the loss follows from free-energy dissipation; it also extends to non-equilibrium states with pointwise orthogonal drift decomposition.Source: §5 · p. 28
  • Numerical experiments report effectiveness across low- and high-dimensional problems, rough potentials, and noisy Monte Carlo observations.The reported noise tests consistently favor the weak/variational hybrid over the corresponding strong-form formulation.Source: §5 · p. 28
  • A family of test functions also treats non-gradient stationary problems without pointwise orthogonality; the four-dimensional example has nonzero stationary current.In that example, the drift is not pointwise orthogonal to the potential gradient.Source: §5 · p. 28
  • Extending the method to time-dependent problems remains a question for future study.Source: §5 · p. 29

A.1 Problem formulation

The formulation seeks an equilibrium stationary density for a Fokker–Planck equation on a bounded domain. Its variational loss characterizes the equilibrium up to scale, while Monte Carlo data determine that scale.

A.1 Problem formulation

  • The stationary problem uses flux J = −D∇u −u∇V and seeks an equilibrium with J = 0.The setting is a bounded Lipschitz domain with diffusion coefficient D > 0 and potential V ∈ C1(Ω).Source: §A.1 · p. 29
  • The equilibrium density is u∗ = Z−1e−V/D and is normalized to have unit mass.Source: §A.1 · p. 29
  • The weak loss vanishes for every positive multiple of the equilibrium density, so the Monte Carlo term constrains its scale.The integrals use equally weighted Sobol quadrature in 2D and 4D and importance-sampling weights in 8D; the analysis concerns local minimizers and does not establish convergence of Adam or L-BFGS.Source: §A.1 · p. 29

A.2 Discrete loss and local error estimates

The discrete weak loss combines a penalty for deviations from equilibrium with Monte Carlo data fitting. A local quadratic analysis shows that it cannot increase mean-square error, while errors invisible to the weak operator and the density’s scale require particular care.

A.2 Discrete loss and local error estimates

  • The discrete loss pairs a flux-based weak penalty with Monte Carlo data fitting.The flux is approximated by Rhuh, and the data and quadrature weights enter the discrete formulation.Source: §A.2 · pp. 29–30
  • The weak loss reduces mean-square error components associated with positive eigenvalues of the discrete operator, while zero-eigenvalue components remain unchanged.For the quadratic model, the weak loss is combined with the Monte Carlo data-fitting term; the minimizer is unique under the stated positive-weight and finite-covariance assumptions.Source: §A.2 · p. 30
  • In the quadratic model, the combined-loss minimizer has mean-square error no greater than the original Monte Carlo error.The operator contracts error components according to its eigenvalues; errors in its nullspace are unchanged.Source: §A.2 · p. 30
  • The weak loss cannot determine the multiplicative constant because it vanishes for every c u∗; Monte Carlo fitting determines c.If Rheh = 0, the quadratic minimizer retains eh, so adding the weak loss does not reduce that error.Source: §A.2 · pp. 30–31
  • For the quadratic model, the mean-square error ratio tends to zero only under a condition on the discrete operator’s nullspace dimension.The necessary condition given is dim(ker Rh)/mh → 0.Source: §A.2 · p. 31

A.3 Mean-square error estimate for the nonlinear loss

The nonlinear-loss analysis extends the error bound to positive local minimizers. Interpolation toward the discrete equilibrium provides the key comparison with the original Monte Carlo error.

A.3 Mean-square error estimate for the nonlinear loss

  • The estimate applies to positive local minimizers of the original discrete loss.Source: §A.3 · p. 31
  • The proof compares the local minimizer with points interpolated toward the positive discrete equilibrium, along which the weak loss does not increase.The resulting inequality bounds the minimizer’s error by the original Monte Carlo error.Source: §A.3 · p. 31
  • At every positive local minimizer of the original discrete loss, incorporating the weak loss does not increase mean-square error.This extends the mean-square error bound established for the quadratic model.Source: §A.3 · p. 31

A.4 Local error analysis of the neural-network approximation

The analysis derives a mean-square error bound for a local quadratic model of the neural-network loss and extends it to linearized network outputs at observation points. The bound does not require the output Jacobian to have full column rank.

A.4 Local error analysis of the neural-network approximation

  • Near θ*, the network outputs and weak loss are assumed twice continuously differentiable, and the weak-loss Hessian is positive semidefinite.The Hessian is positive semidefinite because the weak loss is nonnegative and vanishes at θ*.Source: §A.4 · pp. 31–32
  • The analysis linearizes network outputs at θ* and retains the quadratic weak-loss term to estimate mean-square error.The Monte Carlo error is defined from the observed values and the target density at the observation points.Source: §A.4 · pp. 31–32
  • Every minimizer of the local quadratic model satisfies a mean-square error bound without requiring the output Jacobian to have full column rank.The bound extends to the linearized network error JDδθ at observation points.Source: §A.4 · p. 32

A.5 Uniform convergence of the quadrature approximation

The quadrature approximations of the weak loss converge uniformly over compact parameter sets when the network density and its gradient are jointly continuous and positive. Importance sampling has almost-sure convergence under the stated sampling conditions.

A.5 Uniform convergence of the quadrature approximation

  • Sobol points generated in a containing box and restricted to Ω also support the stated conclusion for bounded Lipschitz domains.Source: §A.5 · p. 33
  • Sobol quadrature uses affinely rescaled Sobol points on rectangular domains, while importance sampling uses iid draws from q > 0 almost everywhere with weights 1/q.For importance sampling, the strong law gives convergence of the weighted sums; normalization factors cancel in the ratio.Source: §A.5 · p. 33
  • Uniformity follows by approximating the compact families of numerator and denominator integrands with finite collections, then passing to their quotient.The normalized denominator is uniformly bounded away from zero for sufficiently large N.Source: §A.5 · p. 33
  • Sobol and importance-sampling approximations converge uniformly over θ ∈ Θ to the weak loss, with almost-sure convergence for importance sampling.The result assumes a compact parameter set and positive, jointly continuous uθ and ∇uθ on Ω × Θ.Source: §A.5 · pp. 32–33

A.6 Conclusion

The weak loss identifies the stationary density only up to positive scale, while Monte Carlo data determine that scale; the combined discrete loss does not increase mean-square error. The analysis also establishes quadrature convergence, but its neural-network error guarantee is limited to a local quadratic approximation.

A.6 Conclusion

  • The combined discrete loss has mean-square error no greater than the Monte Carlo data, and the quadratic model gives the exact ratio of the two errors.Source: §A.6 · p. 33
  • The weak loss is zero exactly for positive rescalings of the stationary density, so Monte Carlo data determine the proportionality constant.Source: §A.6 · p. 33
  • For neural networks, the error bound is proved only for the local quadratic approximation, not for the algorithm's computed network or L-BFGS convergence.Source: §A.6 · p. 33
  • The weighted sums used in computation converge uniformly in θ ∈ Θ to the integrals defining Φ as quadrature points increase.Source: §A.6 · p. 34
Loading 2610.03325v1…