chc.deep_galerkin

Stability tier: experimental.

Deep Galerkin Method -- a neural PDE solver, meeting classical Galerkin/Marchuk FEM on one PDE, and coupled forward-backward into a mean-field game.

The Deep Galerkin Method (Sirignano-Spiliopoulos) trains a network to satisfy a PDE by minimising its residual at random points -- a mesh-free Galerkin scheme. solve_poisson_dgm applies it to the same 1-D Poisson BVP -V''(x) = f(x), V(0)=V(1)=0 that chc.galerkin solves with a variational-difference FEM (progonka), so the neural Galerkin can be checked against the analytic and the classical one. The bridge from classical Marchuk/Galerkin FEM to learning-based PDE solvers.

The second half is the mean-field game that Poisson solve was only a stepping stone to. A backward HJB for the value V and a forward Fokker-Planck for the density rho, joined by the optimal feedback alpha* = -(b/r) V_x and by the population mean m(t) = int x rho(x,t) dx, are solved as one coupled system by DGM (cf. Carmona & Lauriere, SIAM J. Numer. Anal. 59(3), 2021) -- solve_mfg_dgm. Both boundary conditions are hard constraints rather than penalties: rho(0,.) is the prescribed Gaussian by construction and V(T,.) is the prescribed terminal cost evaluated at the network's own terminal mean.

The falsifiable gate is that the linear-quadratic case has a CLOSED FORM: LQMeanFieldGame.solve returns it exactly, from a stationary Riccati root plus a 2x2 trace-free two-point boundary value problem. That reduction also exposes an obstruction -- the fixed point degenerates where one denominator vanishes, always on the anti-monotone branch c > 1 + r a^2/(q b^2) -- which is what lq_mean_field_certificate measures the DGM against. Derived in validation/lq_mean_field.mac, proved in proofs/lq_mean_field.v.

PopulationNoise = Literal['independent', 'common'] module-attribute

Whether each agent gets its own Brownian path or the whole population shares one.

MeanFieldOptimizer = Literal['adam', 'lbfgs'] module-attribute

Adam on freshly drawn collocation points every step, or L-BFGS on one fixed draw.

ScalarMLP

Bases: Module

A small tanh MLP x -> V(x) (scalar in, scalar out) for a 1-D field.

LQMeanFieldGame dataclass

A scalar linear-quadratic mean-field game with a closed-form solution.

dX = (a X + b alpha) dt + sigma dW under E[int (q/2)(X - c m)^2 + (r/2) alpha^2 dt + (q_T/2)(X_T - c_T m_T)^2], with m(t) the population mean. The terminal weight q_T is pinned to the stationary Riccati root so the Riccati coefficient is constant and the remaining two-point boundary value problem has an exact matrix-exponential solution -- that is what makes this a gate rather than itself a numerical solve. The general case has the same structure with a time-varying transition.

riccati_root property

The stationary root of (b^2/r) P^2 - 2 a P - q = 0; also the terminal weight.

closed_loop_rate property

A = a - b^2 P/r = -sqrt(a^2 + q b^2/r) -- negative whatever the sign of a.

branch_threshold property

The coupling c above which the transition matrix turns oscillatory.

lambda_squared property

lam^2 = A^2 - q c b^2/r = a^2 + (q b^2/r)(1-c); negative on the oscillatory arm.

obstruction_gain property

k = A + q_T c_T b^2/r. The fixed point survives every horizon when k < lam.

obstruction_horizon()

First horizon at which the mean-field fixed point degenerates, or inf if none.

On the oscillatory branch this is finite for every terminal weight; on the real branch it is finite exactly when k > lam. Proved in proofs/lq_mean_field.v.

fixed_point_denominator(horizon=None)

cosh(lam T) - k sinh(lam T)/lam -- zero exactly at the obstruction horizon.

solve(n_time=401)

The exact solution on a uniform time grid.

Z is the one component without a closed form (it is a quadrature of S^2 and m^2); it is integrated by the trapezoid rule on the same grid and affects only the level of V, never the control.

MeanFieldSolution dataclass

The exact LQ mean-field equilibrium: V = P x^2/2 + S x + Z against a Gaussian rho.

value_gradient(t, x)

V_x = P x + S -- exact, no quadrature, and the only part the control depends on.

finite_population_rms(n_agents)

sqrt(v(t)/N) on times: what N agents cost against the density.

Exact at every t and every N, not asymptotic. A mean-field feedback is a function of an agent's own state and of deterministic coefficients, so the closed-loop agents are independent Ornstein-Uhlenbeck processes and m_N is a sample mean of N i.i.d. draws; the 1/sqrt(N) is then an identity and the constant is the closed-loop variance. Derived in validation/finite_population_gap.mac, measured by finite_population_gap_certificate.

simulate_population(n_agents, *, replicates=1, seed=0, n_step=400, mean_feedback=0.0, noise='independent')

Run n_agents under this feedback; return (times, empirical mean paths).

The plant stepped here is the primitive one, dX = (a X + b alpha) dt + sigma dW by Euler-Maruyama. Only the control is taken from the solution, so the paths are an independent check on finite_population_rms rather than a restatement of it; the Euler bias on the variance is |A| dt relative, two orders below the Monte-Carlo error at the default step.

mean_feedback adds (kappa/b) m_N to every agent's control. That is the difference between an N-player equilibrium and a mean-field control applied to N players: it couples the agents, moves the closed-loop rate to A + kappa and biases the mean away from m. noise="common" drives every agent with one Brownian path instead of its own, which is the independence assumption the identity rests on -- it is exposed so it can be violated on purpose.

hjb_residual(t, x)

Residual of the backward HJB at the closed form. Zero up to the Z quadrature.

fokker_planck_residual(t, x)

Residual of the forward Fokker-Planck at the Gaussian, in log form (divided by rho).

FieldMLP

Bases: Module

A small tanh MLP (t, x) -> R for a space-time field.

MeanFieldDGM

Bases: Module

A Deep Galerkin solution of the coupled HJB / Fokker-Planck system.

Both boundary conditions are structural rather than penalised: the density carries the prescribed initial Gaussian as an exact factor, and the value function carries the prescribed terminal cost evaluated at the network's own terminal mean -- which is where the mean-field coupling enters the value side. What is left to minimise is only the two interior residuals.

quadrature property

The integration nodes, derived rather than stored.

An array field of an eqx.Module is an inexact-array leaf, so eqx.filter(model, eqx.is_inexact_array) hands it to the optimiser along with the weights -- and an optimiser that moves the quadrature grid corrupts every mean and mass computed on it, silently and without touching the density.

mean(t)

int x rho / int rho on the quadrature grid -- normalised, so mass drift cannot contaminate the coupling.

MeanFieldCurve dataclass

Evidence for Result 49 -- the coupled solve, and what its own residual fails to see.

BlindnessSweep dataclass

Result 55 (d): one solve per horizon on the anti-monotone instance, scored three ways.

Every score is read off the trained model; exact_slopes, errors and control_errors need the closed form, and are what the scores are judged by. The ranks are Spearman's, each score against errors.

FinitePopulationGap dataclass

Evidence that the mean-field density is priced at finite N, not merely approached.

CongestedMeanFieldGame dataclass

LQMeanFieldGame with a cubic congestion coupling: the reduction survives, the linearity does not.

The running cost pays for distance from a congestion-shifted target, (q/2)(x - c m - gamma m^3)^2, so the value function stays quadratic in x -- the Riccati coefficient never sees gamma -- while the reduced two-point boundary value problem for (m, S) becomes nonlinear:

m' = A m - (b^2/r) S,   S' = q (c m + gamma m^3) - A S,   S(T) + P c_T m(T) = 0

That is the smallest change that breaks LQMeanFieldGame.solve's closed form without breaking the reduction it rests on, which is what makes it a usable gate for a generic dual-weighted estimator. At congestion = 0 it is the LQ game exactly, and adjoint_weighted_error is then exact rather than first-order -- Result 55, recovered as a special case. Derived in validation/nonlinear_dwr.mac.

terminal_row property

B_T = (P c_T, 1), the row the terminal condition sends to zero. gamma-free.

reduced_field(state)

F(m, S), in jax so the adjoint is jax.vjp of this and of nothing else.

trajectory(slope_initial, times)

RK4 the reduced system from (m_0, slope_initial); shape (len(times), 2).

solve(n_time=401)

Shoot on S(0) until the terminal row vanishes; return (times, trajectory).

The LQ closed form supplies the first iterate, so the secant converges in a handful of steps even at a congestion that moves the answer by tens of percent. There is no closed form to fall back on here -- that is the point of the class, and it is why this returns the two arrays rather than a MeanFieldSolution: three of that type's five fields do not exist without one.

NonlinearDwrCurve dataclass

Evidence for A16: the adjoint built by autodiff, and what the linearisation costs.

solve_poisson_dgm(source, width=32, steps=4000, n_collocation=128, seed=0)

Deep Galerkin solve of -V''(x) = source(x) on [0,1], V(0)=V(1)=0.

Minimises the mean-squared PDE residual at random collocation points plus the boundary term.

solve_mfg_dgm(game, width=32, steps=3000, n_time=16, n_space=16, n_quadrature=96, learning_rate=0.002, fp_weight=1.0, mass_weight=1.0, half_width=None, seed=0, optimizer='adam')

Solve the coupled forward-backward mean-field system by Deep Galerkin.

By default the spatial box comes from problem data alone -- the initial mean, the initial variance and the stationary variance sigma^2/(2|A|) -- never from the closed-form solution, so the comparison in lq_mean_field_certificate stays honest. Pass half_width to widen it when the population is known to travel further than that; the certificate uses it near the obstruction horizon deliberately, to hand the solver a domain it could not have guessed and show that this still does not save it.

optimizer="lbfgs" runs steps L-BFGS iterations with a zoom line search on ONE draw of n_time x n_space collocation points -- the draw Adam would have taken first -- because a quasi-Newton curvature estimate assumes the objective does not change between iterations, and resampling every step would change it. learning_rate is Adam's alone; the line search picks L-BFGS's step.

dual_weighted_error_estimate(game, model, n_time=401)

An a-posteriori estimate of |S_hat(0) - S(0)| from the model alone -- Result 55.

Result 49 measured a residual that FALLS as the answer degrades near the mean-field obstruction, so a residual-based stopping rule reports its cleanest convergence exactly where the solution is worst (under Adam; the fall is the optimiser's, see below). The remedy is not a bigger residual budget but the right functional. Because the reduced fixed point is AFFINE, the error is an exact quotient (validation/mean_field_dwr.mac STEPs 2-3, proofs/mean_field_dwr.v)::

S_hat(0) - S(0) = eps/den(T) - int_0^T z(s).g(s) ds,   z(s) = Phi(T-s)^T v / den(T)

with g = y' - M y the reduced defect of the model's own y = (m, S), eps = v.y(T) the terminal-consistency defect (structurally zero here), and z the exact ADJOINT solution. Everything is read from the trained network and the cost parameters: the transition matrix is integrated from the model's OWN closed-loop rate A = a - b^2 V_xx/r, never from the reference solution, so this is an estimator and not an oracle comparison.

Two facts make it work where the raw residual does not. The dual weight carries 1/den(T) and nothing else that can blow up, so it diverges exactly at the obstruction; and the reduced residual is homogeneous of degree 1 in (m, S), so a bounded approximator facing a diverging solution keeps a bounded residual by construction (proofs/mean_field_dwr.v, bounded_approximator_is_blind). Conditioning the residual by 1/|den| alone is NOT enough -- measured rank correlation with the error 0.41 against -0.67 for the raw residual and 1.00 for this estimator; what matters is the projection onto the adjoint mode, not the scalar rescaling.

Those three numbers are one configuration. Over five seeds, widths 32, 64 and 128, and Adam against L-BFGS (residual_blindness_sweep; just paper-4 in causaldyn-bench), this estimate ranks the error perfectly in 28 of 30 configurations and at 0.976 in the other two; conditioning is positive in all 30 and perfect in none; and the raw residual ranks the error backwards in 15 of 15 Adam configurations and FORWARDS in 13 of 15 L-BFGS ones -- its size is bounded by the homogeneity, its order is decided by the optimiser. As a value rather than a rank the estimate is within 12% of the true error under Adam and within 80% under L-BFGS: rank with it anywhere, quote it as a number only where it was measured.

Scope: the estimate is exact for this affine family up to quadrature and autodiff error (measured within 6% at a horizon where the error is 72 and the raw residual is at its smallest). For a non-quadratic problem the same construction is the standard dual-weighted residual and is first-order, not exact.

lq_mean_field_certificate(steps=2500, seed=0)

Measure the coupled DGM against the exact LQ mean-field equilibrium.

Three arms. The first checks that the gate itself is exact -- the closed form annihilates both PDE residuals, and the obstruction horizon bisected from the denominator matches the formula. The second checks that the neural solve reproduces it on a monotone instance. The third is the one that can fail loudly: on the anti-monotone branch, as the horizon approaches the obstruction, the DGM's error grows while its own residual shrinks, so a residual-based stopping rule reports success exactly where the answer is worst.

residual_blindness_sweep(horizons=(0.3, 0.4, 0.5, 0.6, 0.68, 0.72, 0.76, 0.79), *, width=32, optimizer='adam', steps=2500, seed=0)

Rank the error of one solve per horizon by three scores as the horizon nears T*.

lq_mean_field_certificate asserts the blindness on two horizons; this is the curve behind it, on the same anti-monotone instance and in the same box, so the horizon is the only thing that changes (T* = 0.8036). The three scores are the raw residual, the residual conditioned by 1/|den(T)|, and dual_weighted_error_estimate -- which reads everything off the model and never calls game.solve().

What is and is not proved about the ranks. proofs/mean_field_dwr.v proves that a model whose reduced state stays within B keeps its reduced residual under a bound that does not depend on the horizon, while its error exceeds any level close enough to T* (bounded_approximator_is_blind): the raw residual cannot rank the error there. That it ranks it BACKWARDS -- rank_residual < 0 -- is not implied, and is what this sweep measures, per seed, width and optimizer. Measured over seeds 0-4 and widths 32, 64, 128: it does under Adam in 15 of 15 configurations and does not under L-BFGS in 13 of 15.

finite_population_gap_certificate(replicates=400, seed=0, n_step=400)

Price Result 49's modelling gap: the Fokker-Planck density is the N -> infinity limit.

Four arms, in increasing order of what they can tell you.

The first measures the empirical-mean gap of a simulated population against the closed form sqrt(v(t)/N), and fits the exponent. The second is the falsification: drive every agent with one shared Brownian path and the exponent must collapse to zero, because independence -- not large N -- is what makes the identity exact.

The third and fourth ask what an N-player equilibrium would change. A mean-feedback gain kappa moves the closed-loop rate to A + kappa: arm three checks the variance against the moved constant and against the unmoved one, so the second number is there to be wrong. Arm four takes kappa = c/N, which is the scaling a real Nash equilibrium has, and measures the resulting bias under common random numbers -- it is O(1/N), so it never catches up with the O(1/sqrt N) fluctuation, which is why the exponent is the part that transfers.

Derived in validation/finite_population_gap.mac; the exponent's Monte-Carlo standard error is derived there too, and reported here rather than assumed.

adjoint_weighted_error(field, times, trajectory, rate, terminal_row, *, free_component=1)

Estimate J(yhat) - J(y) for J(y) = y(0)[free_component], from yhat alone.

The generic form of Result 55. Pairing the linearised primal with the transposed adjoint gives d/dt (z . delta) = z . g for any field (validation/nonlinear_dwr.mac STEP 2, residual 0), and integrating it across the interval leaves

J(yhat) - J(y) = (eps - int_0^T z . g dt) / z(0)[free_component]

with g = yhat' - F(yhat) the defect the caller supplies as rate, eps = B_T yhat(T) the terminal defect, and z the solution of -z' = F'(yhat)^T z from z(T) = B_T. The only thing needed of F is a vector-Jacobian product, which is jax.vjp -- so the hand-written affine transition matrix of dual_weighted_error_estimate is replaced by a linearisation of whatever field is passed in.

free_component names the component of y(0) the boundary value problem does not pin; the other one is an initial condition and cannot be in error. The returned value is SIGNED, because for a nonlinear field the sign of the remainder is part of the measurement.

Exactness is a property of the field, not of this function. For an affine F the second variation vanishes and the formula has no remainder at any perturbation size; otherwise the remainder is O(eta^2) in the size of the error, so the estimate is first-order. nonlinear_dwr_certificate measures both exponents rather than asserting them.

nonlinear_dwr_certificate(congestion=0.6, n_time=801, seed=0)

Measure the generic dual-weighted estimator, and price the linearisation it rests on.

Result 55 shipped an exact error quotient for the LQ mean-field game and recorded that the exactness belonged to the affine reduced problem. This replaces the hand-written affine adjoint with jax.vjp of whatever field is handed over, and measures the four statements that makes, on a manufactured error yhat = y + eta p whose true value is known exactly.

The perturbation is manufactured rather than taken from a trained network deliberately: a network's error is not known, so it could not separate "the estimate is first-order" from "the network is bad". Here the true error is eta p_S(0) by construction and the only thing being measured is the estimator.

seed is accepted for signature parity with the other certificates in this module and is unused: nothing here is random.