Making the non-differentiable differentiable

01: Motivation

Why differentiability is the organizing goal

Modern machine learning is, at bottom, gradient descent on a computational graph. A discrete choice, hard threshold, collision test, or external solver can block that graph and separate what can be learned end-to-end from what must be hand-tuned, sampled, or bolted on afterward. The whole enterprise reviewed here is the project of dissolving those walls.

A differentiable simulator computes gradients of physical states with respect to inputs such as control signals, poses, material properties, and boundary conditions. That enables end-to-end optimization and tight coupling with learning and control. The payoff is threefold: inverse problems (recover parameters from observations) become gradient descent problems, reducing expensive black-box search; control and design can be optimized directly through the physics; and hybrid solvers become possible, where classical numerical methods run alongside trainable neural components and the gradient flows through the whole pipeline. 1 5

Two senses of "differentiable" are worth separating up front, because most confusion in this field comes from conflating them.

Differentiable in principle means a gradient exists almost everywhere and autodiff will hand you a number. That number can still be useless. Three everyday examples, all of which technically "have a gradient": 8 26 36

• A hard threshold or mask, m = (x > 0). Its derivative is exactly zero everywhere it is defined, so autodiff faithfully returns 0. That valid number carries no direction for changing x. Step along it and the optimizer stays put.
• A long chaotic rollout, such as a many-step pendulum or turbulent flow. A number comes back, often exploded or vanished: tiny perturbations compound, so the direction becomes wildly amplified and unhelpful.
• A single-sample score-function estimate. It is unbiased on average. On any given step, its variance is so large that the returned direction is nearly random.

Differentiable in a useful way means the gradient is informative: nonzero where you need to move, finite, low-variance, and well-conditioned enough that stepping along it actually lowers the loss. Replace that hard threshold with a smooth sigmoid and the gradient now points toward the decision boundary; an MSE loss on a smooth model gives a gradient that reliably reduces error. Almost every tool in Section 3 is a device for turning a case from the first list into a case from the second. 12 13

Much of the sophistication in the field targets the second condition. A function can be technically differentiable and still yield gradients that are worthless for optimization.

02: Lineage

The historical and architectural arc

Foundations: autodiff and differentiable programming

The substrate for everything else is automatic differentiation (AD) and the broader "differentiable programming" paradigm: entire pipelines beyond neural nets are implemented so reverse-mode AD can backpropagate through them. This fused backpropagation (from deep learning) with adjoint methods (from the numerical solution of differential equations). Once solvers themselves are written in AD-capable frameworks (JAX, PyTorch, Julia), the boundary between "simulator" and "model" begins to dissolve. Traditional engines resisted this because they often estimate gradients by finite differences (incompatible with backprop, poorly scaling) and use explicit integration with penalty contact (non-smooth). Accurate, cheap gradients arrived when solvers exposed differentiable internals and avoided finite-difference probes from outside. 1 2

The generative era (VAE / GAN): the discrete-sampling problem

The first wall hit at scale was stochastic nodes with discrete latent variables. Categorical variables are the natural way to represent discrete structure, yet a categorical sample blocks backpropagation. This motivated most of the gradient-estimation toolkit in Section 3: score-function/REINFORCE (unbiased, high variance), reparameterization (low variance, originally continuous-only), and continuous relaxations such as the Gumbel-Softmax / Concrete distribution (Jang et al.; Maddison et al., 2016), which propagated into VAEs, discrete-sequence GANs, neural architecture search, and combinatorial optimization. The era's lasting contribution was a durable set of estimators. 12 13

The transformer / operator-learning era

In parallel, the transformer became the dominant inductive-bias-light backbone for neural surrogates. More consequentially for physics, neural operators arose, learning function-to-function mappings (input field → solution field): DeepONet (Lu et al.), a branch/trunk decomposition grounded in an operator-approximation theorem; and the Fourier Neural Operator (Li et al.), a spectral method, generally the strongest baseline across PDE benchmarks that degrades where sharp discontinuities create Gibbs artifacts. These are natively differentiable, so the question shifts from "does the gradient exist" to "does the learned operator respect physics and stay stable under rollout." 19 20

The diffusion / score-based / flow-matching era

The current frontier for generative physics modeling is diffusion and its relatives (VAE → GAN → normalizing flows → score-based diffusion → flow matching / stochastic interpolants). In scientific ML these serve as fast surrogates (calorimeter-shower and particle-physics fast simulation), probabilistic field generation and reconstruction (turbulence, PDE emulation, diffusion-based ensemble weather forecasting, score-based data assimilation), and physics-informed / hard-constrained sampling. The conceptual shift is to reframe simulation as sampling from a learned distribution, a natural fit for physics that is genuinely stochastic or under-determined by sparse data. 21

The differentiable-simulator arc (running in parallel)

Independently, a simulator-centric line matured: analytical LCP-contact gradients for simple rigid bodies (de Avila Belbute-Peres et al., 2018; Degrave et al., 2019); high-DOF and deformable systems including fluids, elastic bodies, cloth, and soft multi-body dynamics (Hu et al.; Liang et al.; Qiao et al.); end-to-end perception-to-physics via differentiable rendering + simulation (Jatavallabhula et al., 2021); and platform-level differentiability (MuJoCo XLA, Tiny Differentiable Simulator, Warp, DiffPD, Nimble, Genesis), plus domain codes for accelerators (Cheetah), subsurface flow (DPFEHM), neuroscience (Jaxley), and detectors. The two arcs meet wherever a learned component and a physics solver share one gradient path. 3 4 5 6 7

03: Toolkit

The hard components and the tools for each, illustrated

Real pipelines rest on algorithmic components implementing discrete decisions or relying on discrete intermediate structures; these break gradient flow, and the fix is a differentiable proxy via smoothing or relaxation. Every tool below is really a different answer to one question: what do I send backward when the forward op has no usable slope? That shared problem looks like this:

Tool 1: Continuous relaxation (Gumbel-Softmax / Concrete)

Send a relaxed vector forward. First reparameterize the randomness out of the way with the Gumbel-Max trick (a categorical draw equals argmaxᵢ(log πᵢ + gᵢ) with fixed Gumbel noise gᵢ = −log(−log uᵢ)), then relax the hard argmax into a temperature-scaled softmax:

The temperature τ is the dial between a smooth blur and a near-discrete sample. It is biased (temperature-dependent), low-variance, and standard for categorical latents you control. 12 13 14

In the lineageThe generative (VAE / GAN) era invented this for discrete latents, where categorical draws block ordinary backprop. The same trick carried into hard attention and neural-architecture search in the transformer era.

Tool 2: Straight-through estimator (STE)

The cheapest trick, behind quantized networks and VQ-VAE, keeps the hard op in the forward pass so the network really sees a discrete z. On the way back, it substitutes a convenient derivative:

Forward and backward literally disagree about what the op is. The estimator is biased, cheap, and often good enough. 15 27

In the lineageStraight-through lets VQ-VAE and VQ-GAN codebooks, the discrete tokenizers under modern transformers and latent diffusion, send a gradient to their encoder while emitting genuinely discrete tokens forward.

Tools 3 & 4: Score-function vs. reparameterization

Both differentiate the same expectation, ∇θ E[f(z)]. The fork determines what kind of downstream function you can use:

Score-function treats f as a black box, so it works for non-differentiable rewards and pays for that flexibility with high variance. Reparameterization follows a smooth sample path, giving low variance when f is smooth and pθ is reparameterizable (Gaussian VAEs; Gumbel-Softmax is the discrete cousin). 12 26

In the lineageReparameterization made the Gaussian VAE trainable and forms the mathematical core of diffusion, where every denoising step reparameterizes Gaussian noise. Score-function is the generative era's fallback for discrete or black-box outputs (early discrete-latent inference; sequence models trained on non-differentiable metrics), now reused to fine-tune diffusion models on non-differentiable rewards.

Tool 5: Randomized / stochastic smoothing

The workhorse for contact and friction in physics, hard renderers, and black-box solvers averages the function over random input perturbations, convolving sharp edges into a smooth surface:

Stein's identity gives the gradient of the smoothed function using only evaluations of f. No derivative of f is needed, so it works even where f is a flat staircase or a black box. Radius σ is the dial: σ → 0 recovers the sharp function; larger σ is smoother and more biased. 10 11

In the lineageThis is the differentiable-simulator arc’s answer to contact and friction, and the route by which diffusion guidance reaches through a non-differentiable classifier or reward.

Tool 6: Implicit differentiation / adjoint

The one that matters most for physics, because simulators run many steps. Whenever the output is the solution of an equation, such as a fixed point, an optimization optimum, or an integrated ODE/PDE, you can either unroll (backprop through every iteration, storing every state, memory growing with steps) or differentiate the solution implicitly. If z* satisfies F(z*, θ) = 0:

This powers neural ODEs, differentiable optimization layers (OptNet), and adjoint-method PDE solvers. 17 18

In the lineageIt powers the transformer era's implicit-depth models (deep equilibrium networks, neural ODEs), the probability-flow ODE behind diffusion sampling, and adjoint PDE solvers in differentiable simulators. Use it for equation-defined outputs.

The taxonomy behind the pictures

The hard components fall into eight categories: (1) discrete/categorical sampling; (2) hard selection: argmax, top-k, quantization, rounding, integer casts; (3) combinatorial operators: sorting, ranking, matching, shortest paths, dynamic programming; (4) control flow and thresholds: branches, masks, cuts, Heaviside steps; (5) non-smooth physics: contact, collision, friction, fracture; (6) rendering/rasterization; (7) black-box/external components; and (8) long-horizon/chaotic composition, a pathology of the full chain. The tool families map onto them as follows.

Which era reaches for which tool, and why

Because the tools were invented to solve concrete problems in the architectures of Section 2, each era leans on a characteristic subset. Several tools recur across eras, which is why the same six ideas keep appearing.

Routing a hard component to a tool

The top split is the kind of hard component: discrete choice, non-smooth function, equation-defined output, or continuous random node. For a discrete choice, use score-function with a variance-taming baseline when the downstream value is a black box or non-differentiable reward; use straight-through when the forward pass must stay genuinely discrete; use Gumbel-Softmax when a soft relaxation can be annealed toward the discrete limit. For a non-smooth or black-box function, exploit an analytic/structured gradient where one exists, fall back to randomized smoothing when no internals are usable, or replace the component with a learned surrogate when it is swappable. An equation-defined output goes to implicit / adjoint. A continuous stochastic node goes to reparameterization.

Long-horizon or chaotic composition belongs to the whole chain. It rides on top of whatever tool you choose and is handled with truncated backprop, checkpointing, gradient clipping, and the audit checks below. 36

04: Gradient audit

How to know whether the differentiable tool is telling the truth

Before asking whether a differentiable pipeline is accurate, ask whether its gradient is worth trusting. A forward simulator can match observations and still produce a gradient that points the optimizer the wrong way. Treat differentiability as a claim that has to be audited.

4.1  Start with the baseline and cost audit

Begin with stronger questions than "Can I differentiate this?" What is the classical baseline? For PDE-constrained problems, that is often a classical adjoint; finite differences or random search are weaker baselines. Is this a repeated-query problem where gradients will be reused thousands of times? And what is the total cost to a validated result once simulations, solver rewrites, training, hyperparameter search, final high-fidelity checks, and human tuning are counted? A paper that reports only neural inference time, while hiding the expensive steps needed to reach a validated answer, has not shown a real speedup. 17 18 37

4.2  Classical adjoint first

For smooth PDE inverse problems and many engineering design loops, the right baseline is adjoint PDE-constrained optimization, which often computes sensitivities at a cost nearly independent of the number of design variables. Differentiable ML can still help as an implementation route for adjoint-style gradients inside AD frameworks, as an amortization layer when many related solves are needed, or as a coupling layer when a neural component must share one gradient path with the solver. 17 18

4.3  Gradient checks that should be routine

At minimum, audit the gradient before trusting any inverse, control, or design result: random-direction finite differences, adjoint consistency on small problems, gradient-to-loss correlation under small steps, norm and signal-to-noise tracking, sensitivity to timestep, mesh, solver tolerance, relaxation temperature τ, smoothing radius σ, and rollout horizon, plus direct probes near contacts, collisions, thresholds, and topology changes. The final check is whether the candidate that optimizes the relaxed or learned surrogate still works in the high-fidelity solver that actually matters. 6 8 17 18 37

4.4  Identifiability comes before optimization

For sparse and noisy experimental inverse problems, a gradient can optimize an under-identified parameter. Physics regularization only narrows information already present in the experiment. Before reporting a point estimate, check sensitivity rank, parameter correlations, posterior width, and whether several parameter settings explain the data equally well. If the inverse problem is ill-posed or multi-modal, report uncertainty or posterior samples alongside the optimizer endpoint. 24 25

4.5  Constraints are four different things

Keep these separate: hard constraints, soft penalties, projection or correction into the feasible set, and post-hoc repair after the model has already produced an invalid result. They have different gradients, different failure modes, and different claims. A soft penalty gives a soft objective. A repaired design only shows post-hoc feasibility, leaving the differentiable path itself to be checked.

05: Evaluation

Benchmarks and how to judge success

Evaluation here is unusually treacherous: low pointwise error can coexist with physically nonsensical output, and vice versa. A few benchmark suites and a richer-than-usual metric set have emerged.

Forward surrogate / operator suites. PDEBench is the most-cited broad suite (Burgers, diffusion-sorption, diffusion-reaction, incompressible and compressible Navier–Stokes, Darcy, shallow water) with ready datasets and baselines (FNO, U-Net, PINN, gradient-based inverse); FNO is often a strong baseline in these suites and degrades where sharp discontinuities create spectral artifacts. APEBench targets autoregressive rollout stability; CFDBench and Multiphysics Bench add fluid-scale data and strongly-coupled multiphysics (where all learned solvers degrade); PDEInvBench is a dedicated 2026 TMLR benchmark for PDE inverse problems, with splits for in-distribution and out-of-distribution parameter recovery. 20 22 23

Metrics beyond RMSE. Normalized RMSE, maximum (worst-case) error, conserved-quantity error (does it respect conservation laws?), boundary error, and Fourier-space error (which frequency bands fail: high-frequency failure is the classic surrogate weakness). For diffusion/probabilistic models: calibration and proper scoring rules (CRPS), plus whether the model captures the invariant statistics of the attractor beyond pointwise trajectories. 20 21 23

06: Data & the estimators

What data size and quality mean for learning the non-differentiable parts

For the tools in Section 3, "data" behaves differently from ordinary supervised scaling, because the object being learned through is a gradient estimator sitting on a discrete or non-smooth node. Two quantities matter: (i) the number of Monte-Carlo samples used to estimate the gradient through that node at each step, and (ii) the noise in the signal feeding the node. The right frame throughout is bias vs. variance of the estimator, and data touches only one of the two.

6.1  Data size here means samples-per-gradient, and it only buys down variance

For a stochastic node, "more data" usually means more samples per estimate; training-set size is usually secondary. The canonical Monte-Carlo gradient survey (Mohamed, Rosca, Figurnov & Mnih, JMLR 2020) makes the ranking concrete: the score-function / REINFORCE estimator is unbiased and high-variance, so its error falls only slowly with the sample count N and needs baselines/control variates to be usable; the reparameterization / pathwise estimator reaches far lower variance at the same N (its variance is bounded by the squared Lipschitz constant of the cost). The data lever is the pairing of sample count and estimator: a small N with reparameterization can beat a large N with score-function. 26

6.2  The bias floor: data reduces variance; scarcity hides failure

Relaxation (Gumbel-Softmax), straight-through, and any biased smoothing are systematically biased. Adding samples shrinks the scatter around the wrong answer while the offset remains (the amber floor in Fig 15). You reduce that bias by annealing the temperature τ or shrinking the smoothing radius σ, which re-inflates variance. Bias–variance is a knob you set; extra data alone leaves the tradeoff in place. At low sample counts, the failure can become invisible: Suh et al. (2022) show that a first-order gradient through a discontinuity can have near-zero empirical variance while its true bias is large. Scarce data can conceal a broken estimator. 8 12 26

6.3  Data quality means noise amplification

Noise enters the non-smooth nodes in three ways. It inflates score-function variance (noisy rewards or labels widen an already high-variance estimator). It acts as an extra perturbation on thresholds and contact, raising the effective smoothing radius and therefore the bias. Any loss that differentiates the data, including PDE residuals in PINNs, derivative terms in equation discovery, and contact impulses, amplifies noise because a numerical derivative scales the error like ε/Δx. The robust recipe is denoise, then differentiate: reconstruct a clean, physically-consistent trajectory before taking gradients. Noise also enters implicit/adjoint methods through the solution z*, where it can degrade the conditioning of the linear solve. 24 28

6.4  Discrete capacity is limited by usable gradient signal

The "model size / complexity" question for these tools is really: how big can the discrete component be before it stops learning? For discrete bottlenecks trained with straight-through, including VQ-VAE codebooks, hard attention, and discrete tokens, only the codes that actually receive gradients survive; the rest go dead. This codebook collapse (van den Oord et al., 2017; Zheng & Vedaldi, 2023) means naively scaling the codebook or category count fails: large codebooks collapse to a used fraction unless propped up by reinitialization, EMA updates, or clustered assignment. The usable capacity of a non-differentiable component is set by the portion the estimator can reach. 27

6.5  Synthesis

For the differentiability tools, the outcome is governed by the estimator's bias–variance profile, and the three inputs act on different parts of it: sample count buys down variance (fastest for reparameterization, slowest for score-function); noise raises variance and, through differentiation, bias, so it must be denoised; and discrete capacity is bounded by the gradient signal each part receives. Systematic bias from relaxation, straight-through, and smoothing is controlled through τ, σ, and the annealing schedule. Those knobs deliberately trade back against the variance that sample count was reducing.

07: Engineering the tools

Making the tools actually train: the tricks that make "differentiable" work

Every tool in Section 3 manufactures a gradient where none was usable. In doing so, it perturbs the thing you are optimizing. The perturbation lands in one of two places: the estimator's statistics (it becomes biased or high-variance) or the loss landscape (it becomes ill-conditioned, discontinuous, or imbalanced across terms). The techniques below are targeted controls; each one counters a specific distortion introduced by a specific tool. Read the section as distortion → remedy.

7.1  Relaxation (Gumbel-Softmax): schedule the temperature

Distortion: statistics (bias). A relaxed categorical sample is exact only in the limit τ → 0; at any finite temperature it is a biased stand-in for the discrete draw, and the softmax leaks probability mass to the wrong categories.

Remedy: temperature scheduling. Start hot (τ ≈ 1–5) for smooth, well-behaved gradients and exploration, then anneal toward a small floor (τ ≈ 0.1–0.5); Jang et al. use the schedule τ = max(0.5, exp(−rt)) over the training step, with exponential, linear, or multiplicative decay as alternatives. The schedule is itself a bias–variance dial: anneal too fast and the model collapses prematurely onto a suboptimal discrete choice; anneal too slowly and it stays soft and indecisive. The straight-through Gumbel variant takes a hard argmax on the forward pass, so training and evaluation see the same discrete object, while backpropagating the soft gradient; τ can even be made a learned parameter. 30

7.2  Straight-through & vector quantization: prop the codebook up with auxiliary losses

Distortion: statistics and capacity. STE's forward/backward mismatch is biased, and the discrete bottleneck it feeds tends to collapse: because only selected codes receive gradients, most of a large codebook goes dead (§5.4).

Remedy: commitment loss, EMA, resets. VQ-VAE adds a commitment loss β‖z_e − sg(e)‖² (β ≈ 0.25, where sg is stop-gradient) that pulls encoder outputs toward their assigned codes; strip it out and codebook utilization collapses toward 0% as the codebook grows. The codebook itself uses an exponential-moving-average online k-means update (decay γ), which is empirically faster and more stable than gradient descent; a stop-gradient keeps the encoder and codebook updates decoupled; and a dead-code reset respawns any code unused for many steps as a noisy copy of a live encoder output. k-means initialization and starting with a small codebook help too. 27 31

7.3  Score-function & smoothing: reduce the variance

Distortion: statistics (variance). These estimators are unbiased and noisy; a single sample points in a nearly random direction (§1, §5.1).

Remedy: baselines, more samples, antithetics. Subtract a baseline / control variate b: it leaves the estimate unbiased yet can cut its variance sharply (REINFORCE-with-baseline; REBAR and RELAX build differentiable control variates for discrete variables). Average over more Monte-Carlo samples per step (variance falls roughly like 1/N). For randomized smoothing specifically, antithetic or coupled perturbations and randomized quasi-Monte-Carlo lower the variance at a fixed sample budget. 10 11 32

7.4  PINN soft-constraints: repair the loss landscape

Distortion: landscape. Turning a PDE into a soft penalty builds a multi-term objective with badly imbalanced gradients: the residual term dominates and starves the data and boundary terms, the problem becomes ill-conditioned, and time-dependent PDEs have their causality violated (Krishnapriyan et al. 2021; Wang, Teng & Perdikaris 2021).

Remedy: adaptive weighting, causal scheduling, better optimizers. Adaptive loss weighting reads the back-propagated gradient statistics each step and rescales the term weights λ_i so no term's gradient swamps the others, with an EMA running-average to tame the estimate's variance (learning-rate annealing, Wang et al. 2021, with η = 10⁻³, α = 0.1; relatives include GradNorm, inverse-Dirichlet, SoftAdapt, ReLoBRaLo, and self-adaptive weights); an NTK-based variant balances the per-term convergence rates (Wang, Yu & Perdikaris 2020). Causal / curriculum scheduling forces earlier times to be learned first and ramps difficulty (Wang, Sankaran & Perdikaris 2024; Krishnapriyan's curriculum and sequence-to-sequence reformulations recover 1–2 orders of magnitude in accuracy). Residual-based adaptive sampling concentrates collocation points where the residual is large; Fourier features restore high frequencies; and training typically runs Adam for the bulk, then switches to L-BFGS / second-order for the final, ill-conditioned high-accuracy push. 33 34 35

7.5  Differentiable simulators: tame the Jacobian spectrum

Distortion: statistics, via composition. Chaining differentiable steps multiplies their Jacobians; if the system's Jacobian spectrum contains large eigenvalues, the gradient explodes or vanishes and its variance blows up. Metz et al. (2021, Gradients are Not All You Need) trace this chaos-based failure precisely to that spectrum. Contact adds hard discontinuities on top.

Remedy: truncate, clip, checkpoint, soften. Truncated backprop shortens the unroll. A 400-step Brax rollout can fail outright, while a short window trains; the tradeoff is bias from dropping the long Jacobian-product terms. Gradient clipping / normalization rescales g whenever ‖g‖ exceeds a threshold, capping the explosion. Gradient checkpointing recomputes activations so a long rollout fits in memory. For contact, softening (blurring contact forces over a lengthscale) and chunking the trajectory at contact events restore usable gradients (§8 covers when these remedies fail and a zeroth-order method wins). 8 36

7.6  The five levers, and why Adam is the default

Across all six tools, the same handful of levers recur. Adam is the near-universal default because its per-parameter adaptive rates already absorb some of the gradient imbalance these tools create. It is a cushion, with the targeted fixes above still required.

08: The target regime

Physics-informed modeling from noisy, sparse experimental data

This is the regime the review is ultimately aimed at: integrate real experimental data that is noisy, sparse, expensive, and partially observed, and still recover physically meaningful models and parameters. The first question is whether the observations identify the parameters at all. Sparse sensors, noisy measurements, and hidden states can make many parameter settings equally plausible even when the optimizer fits the observations. Before reporting a point estimate, run sensitivity and rank checks, inspect parameter correlations, compare independent restarts, and ask whether adding a sensor or perturbation would collapse the posterior. Physics regularization can shrink the solution space around information already present in the experiment. 24 25

Why physics-informed methods fit. By embedding governing equations into training, PINNs and relatives are mesh-free, natively handle inverse problems (infer parameters, boundary conditions, closure coefficients from sparse and noisy data that would underdetermine a purely supervised model), and perform data assimilation (each source improves the other: simulation raises the effective resolution of the data; data raises the credibility of the simulation). 24

Handling noise and uncertainty explicitly. The strongest versions are Bayesian or ensemble-based. B-PINNs and related physics-informed Bayesian models carry uncertainty from both the measurement and the physical model, yielding posteriors over weights and over unknown equation parameters. The reporting rule is simple: if the inverse problem is ill-posed, multi-modal, or sensor-limited, report a posterior, interval, ensemble, or calibration diagnostic alongside the best-fit parameter vector. 24 28

Squeezing signal from extreme scarcity. Recent frameworks recover constitutive laws (hyperelastic strain-energy densities) from a single material test using only sparse displacement and reaction-force data, combining PINNs with finite-element discretization, a two-stage Adam→L-BFGS optimization, and sparse (ℓ_p) regression to keep the discovered model parsimonious. Multi-fidelity PINNs fuse cheap low-fidelity physics with a few high-fidelity measurements; transfer learning ports models across regimes when long-term high-fidelity data is unattainable. 29

Differentiable simulators as the assimilation engine. Where a fully differentiable simulator exists, it can be embedded directly in the training loop so gradient-based optimization fits physical parameters to data. Examples include differentiable multiphase reservoir flow, Cheetah for accelerator system identification, and Jaxley for fitting detailed biophysical neuron models to recordings and tasks. The published Jaxley paper reports large GPU-accelerated biophysical networks, multilevel checkpointing, gradients that cost several forward-pass equivalents, and a 2,000-neuron, 1-million-synapse benchmark with 8,000 time steps. Differentiability turns black-box assimilation into gradient-based fitting. The gradient cost and memory strategy still have to be reported. 3 4 5

09: The honest ledger

Where simple ML or engineering intuition wins

Many papers claim that differentiable physics or physics-informed ML "outperforms" a conventional approach. The more useful question is where the advantage is real, where it is fragile, and why. Three lines of evidence draw the boundary sharply. 37

9.1  The weak-baseline problem: much of the reported advantage is an artifact

A systematic review of ML-for-fluid-PDE research (McGreivy & Hakim, Nature Machine Intelligence, 2024) found that among articles claiming to beat a standard numerical method, 79% (60 of 76) compared against a weak baseline, usually an under-tuned or low-resolution classical solver that any competent numericist would improve. Combined with reporting biases (negative results under-reported), the authors conclude the field is overoptimistic: the literature makes ML look more successful at solving fluid PDEs than it actually is. The practical reading: where a mature, fast classical solver already exists for a forward PDE, differentiable or neural surrogates frequently show little speed-accuracy advantage once the baseline is tuned properly. The "win" was often measured against a straw man. 37

9.2  When the gradient itself becomes the problem

The central promise of a differentiable simulator is that its first-order gradient beats derivative-free (zeroth-order) estimates. Suh et al. (ICML 2022, oral) show this promise depends intricately on the physics. For systems that are stiff, discontinuous (contact-rich), or chaotic, the simulator's first-order gradient can be biased and actually hurt optimization, converging to worse or unsafe solutions than a plain zeroth-order or REINFORCE-style estimator. Empirical variance alone misses this bias. Their fix (an α-order estimator interpolating toward zeroth-order near non-smoothness) is essentially a concession: rich contact, chaos, and stiffness are exactly where the nominal first-order gradient is weakest, and simpler sampling-based methods are competitive or better. 8 37

9.3  PINNs stumble on problems classical methods find trivial

Krishnapriyan et al. (NeurIPS 2021) demonstrate that standard PINNs, which impose the PDE as a soft loss penalty, fail to learn even simple convection or reaction dynamics. A textbook finite-difference or finite-element scheme solves those problems in milliseconds. The failure comes from the soft PDE regularizer, which makes the loss landscape ill-conditioned and hard to optimize. Curriculum training and sequence-to-sequence reformulations recover 1–2 orders of magnitude in accuracy. The lesson stands: for a well-posed forward PDE with a good classical solver, a PINN is typically slower to train, less accurate, and less reliable. Its value is in inverse, sparse, or constrained settings. 37

9.4  Cases where simple ML or engineering intuition simply wins

• Abundant, clean data for a forward map. If you have plentiful high-quality input→output pairs and only need a fast approximation, a plain regression surrogate (or even the tuned numerical solver) is usually simpler and just as good. The physics machinery earns nothing extra.
• One-off forward solves. Differentiable pipelines pay off when a gradient is reused many times (inverse design, control, assimilation). For a single forward simulation, building and validating a differentiable stack costs more than it returns. Run the classical solver.
• Well-characterized engineering regimes. Where reduced-order models, dimensionless correlations, or lookup tables already capture the behavior (heat-transfer correlations, drag curves, standard control laws), engineering intuition plus a simple fit matches or beats a heavy differentiable model at a fraction of the effort and with better interpretability.
• Smooth, low-dimensional relationships. If the map is smooth and low-dimensional, linear models, splines, or gradient-boosted trees with a few physically-motivated features are strong, robust baselines. Overoptimistic studies often omit exactly these comparisons.

9.5  When it earns its place

The advantage is real and conditional. It appears when two things hold at once: (i) the landscape is smooth enough that gradients are trustworthy (ruling out the stiff, contact-rich, and chaotic cases above), and (ii) there is either no good classical alternative or an amortization or inverse structure that reuses the gradient many times. Concretely, it wins for repeated inverse-design and control loops, for sparse and noisy inverse problems where physics regularizes an otherwise ill-posed fit and no classical inverse solver exists, and for end-to-end pipelines (perception → simulation → control) where the gradient must cross module boundaries that no classical method spans. A non-smooth landscape, a fast classical solver, or a one-off forward solve shifts the advantage back to the simpler tool. That is the whole map: 8 37

10: Open problems

Challenges facing the field today

The field now has a credible toolkit. Its hardest regimes remain exactly the ones that matter most in practice: contact-rich mechanics, long chaotic horizons, scarce and noisy inverse data, coupled multiphysics, and generative samplers whose deployment cost is still high. The open problems below are the current edge of what gradients can be trusted to do. 5 8 21 23 24 35 36 37

• Useful gradients over merely existing gradients. Null-almost-everywhere, infinite, or high-variance gradients still require per-problem relaxation/smoothing choices, and the smoothing radius or temperature remains problem-specific.
• Correctness and validation of gradients. Contact-rich simulators can be forward-correct yet gradient-wrong; standard gradient-validation benchmarks are only now appearing.
• Long-horizon and chaotic optimization. Exploding/vanishing gradients, decorrelation, and a physics-imposed predictability ceiling. Fixes (checkpointing, gradient decay/clipping, multi-step penalty losses, truncated backprop, intermediate goals) are heuristic.
• Contact, friction, fracture. Non-smooth impulsive physics remains the hardest case; smoothing trades exactness for informativeness and the "right" differentiable contact model is unsettled.
• Data scarcity, noise, and ill-posedness. Inverse problems are frequently ill-posed; physics regularization depends on identifiability already present in the experiment, and calibrated UQ at scale is still expensive.
• Multiphysics coupling. Strong inter-field coupling degrades every current learned solver.
• Scale helps only in limited ways. Geometry, operating conditions, and accuracy demands cap corpus size; combining foundation-model scale with physical structure and stability guarantees is open.
• Cost of differentiable sampling. Diffusion surrogates' many function evaluations remain a deployment bottleneck; consistency and flow-matching variants help and are still becoming standard.

11: Practitioner's guide

A compact decision guide

The shortest honest summary of the field is that the right tool depends less on the buzzword than on the failure mode: whether the hard part is discrete, non-smooth, equation-defined, noisy, or structurally unidentifiable from the data. The table is meant to shorten that routing step before you spend a week building the wrong differentiable stack. 8 12 17 18 24 25 37

• Before differentiating: identify the strong classical baseline, estimate total validated cost, and decide whether this is a repeated-query problem where gradients will be reused.
• Before optimizing: run the gradient audit: directional finite differences, adjoint consistency where possible, gradient-to-loss correlation, resolution and tolerance sensitivity, and final high-fidelity validation.
• Discrete latent or categorical choice? Start with Gumbel-Softmax (+ STE for a hard forward); move to REINFORCE-with-baseline if you need unbiasedness or a black box; use generalized or variance-reduced Gumbel for exotic distributions.
• Hard selection, quantization, or integer physics? Straight-through, or relax during optimization and re-discretize afterward.
• Combinatorial operator? Purpose-built differentiable sorting, ranking, or matching, or stochastic smoothing of the black box.
• Non-smooth physics, rendering, or black box? Randomized or stochastic smoothing, or a learned differentiable surrogate.
• Equation- or optimization-defined layer? Prefer implicit differentiation or adjoint for a converged equation solve; short unrolls can still be useful for short horizons, nonconverged iterative processes, or when solver internals are part of the model.
• Long or chaotic rollout? Truncated backprop + checkpointing + gradient decay; tune the unroll horizon empirically; consider multi-step penalty losses.
• Noisy, sparse experimental data? Check identifiability first; then use physics-informed loss for regularization + Bayesian or ensemble treatment for calibrated uncertainty + multi-fidelity or transfer; denoise before differentiating.
• Multi-modal inverse problem? Use posterior sampling, diffusion or flow priors, or ensembles; report more than a single deterministic design unless the application genuinely needs one.
• Geometry or topology inverse design? Treat the representation as part of the model: constraints, manufacturability, and final-solver validation often matter more than the optimizer.
• Before claiming a win: benchmark against a tuned classical solver and a simple ML baseline; most reported advantages evaporate otherwise (§9.1).