Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Structured Simulators and Neural Operators › Graph-Based Physical Models

Michael BrenndoerferJune 22, 202610 min read

World Models HandbookAll chapters

Part 4 of 12Space, Physics, and Causality

Chapter 22Structured Simulators and Neural Operators › Graph-Based Physical Models

Reading section 4 of 11 · Graph-Based Physical Models · 10 minSection 4 of 11 · 10 min

Section 4 of 11 · 10 min read

Structured Simulators and Neural Operators › Graph-Based Physical ModelsOne section · Sections in this chapter4 / 11 ↓

Sections in this chapter11

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Differentiable simulators assume you know the equations. Graph-based physical models assume something much weaker: you know that the system is a network of local interactions. You do not know the force law, the material model, or even in detail what the entities are; you only know that influence is local, that it happens between pairs, and that the same rule applies everywhere. This is the dynamics-side twin of the object-centric representations from Part III, Chapter 5. There we argued that state should be decomposed into objects and relations; here we build the transition function out of the same decomposition. The payoff of matching the two is that the representation and the dynamics now agree about what the world is made of, and learning can happen in the natural coordinates of the problem.

The canonical particle example is the Graph Network Simulator (GNS) of Sanchez-Gonzalez and colleagues (2020). MeshGraphNets, introduced separately by Pfaff and colleagues, applies related graph machinery to mesh-based simulation. In these systems the "objects" are particles or mesh nodes, the edges connect spatial neighbors, and learned updates stand in for parts of a hand-designed simulator. The inductive bias has three parts:

State is a set of nodes, each carrying local information (position, velocity history, material properties).
Interactions are edges, and the dynamics of a node depend on its neighbors through learned message-passing functions, not on a global flattened vector.
The update predicts a derivative (typically acceleration), which is integrated with a fixed physical scheme such as semi-implicit Euler.

The third part deserves emphasis because it imports a useful kinematic structure. Predicting acceleration and integrating it enforces a second-order update from acceleration to velocity to position. It does not by itself enforce Newton's second law, momentum conservation, or force symmetry; those require additional parameterization or constraints. A model that predicts positions directly has no comparable integration scaffold.

A message-passing round looks like

m_{ij} = \phi\big(h_i, h_j, e_{ij}\big), \qquad h_i' = \psi\left(h_i, \sum_j m_{ij}\right),

where:

$h_i$ : the feature vector (for example, position and velocity history) of node $i$ before the update
$h_j$ : the feature vector of neighboring node $j$
$e_{ij}$ : the edge features carrying relative geometry between the pair, for example $x_j - x_i$
$m_{ij}$ : the message sent from node $j$ to node $i$ , computed by the learned interaction network $\phi$
$\sum_j m_{ij}$ : the sum of all incoming messages to node $i$ , aggregating its neighborhood
$\psi$ : the learned update network that combines a node's own features with its aggregated messages
$h_i'$ : the updated node features after one round of message passing

Here $\phi$ and $\psi$ are small neural networks, $e_{ij}$ carries relative geometry (for example $x_j - x_i$ ), and the sum aggregates all incoming messages. Two design details carry most of the physical weight. First, if all positional inputs enter through relative displacement rather than absolute position, translating the whole system leaves the interaction input unchanged. Other features, such as masses, material types, or relative velocities, may still affect the interaction. Second, summation matches additive force composition when messages are parameterized as pairwise force contributions; it does not alone enforce antisymmetry or conservation. Stacking $M$ rounds lets information travel $M$ hops per step for a fixed graph.

Relational inductive bias

An architectural constraint that encodes the hypothesis "the dynamics decompose into interactions between entities, with the same interaction law shared across entities." Graph networks are the standard carrier of this bias for physical worlds.

The payoff of this bias shows up in two places. First, sample efficiency: the model must learn one interaction law rather than a separate rule for every configuration of the system. A pair of neighboring particles observed in one corner of the domain constrains the same function that governs particles elsewhere, so every pair observation in the training data informs the single shared law. Second, size and shape transfer: a model trained on a chain of 8 masses can drive a chain of 12, because nothing in its parameters is tied to the number 8. The learned function $\phi$ takes a pair of features and returns a message; it neither knows nor cares how many such pairs exist. A fixed-width dense model over a flattened state does not natively accept another entity count; supporting size variation requires padding or masking, a maximum-size parameterization, shared structure, or architectural redesign.

Let us demonstrate the core mechanism in its simplest possible form. We will simulate a chain of masses connected by springs with an unknown rest length, assume we know only that (a) forces are pairwise between neighbors, (b) the same law $\varphi$ applies to every pair, and (c) damping is linear and known, and then recover $\varphi$ from observed motion by linear regression. This strips the message-passing machinery down to its identifiable core: what survives when you remove the neural networks is the hypothesis "one shared pairwise law," and that hypothesis alone is enough to make the problem solvable from modest data.

In[11]:

Code

import numpy as np

rng = np.random.default_rng(1)

k_spring, L0, b_damp, m = 2.0, 1.0, 0.3, 1.0
dt_chain, steps_chain = 0.01, 400


def simulate_chain(x0, v0, steps=steps_chain):
    """Semi-implicit Euler for a free chain with neighbor springs and damping."""
    n = len(x0)
    x, v = x0.copy(), v0.copy()
    xs = np.zeros((steps + 1, n))
    vs = np.zeros_like(xs)
    acc = np.zeros_like(xs)
    xs[0], vs[0] = x, v
    for t in range(steps):
        pair_force = k_spring * (x[1:] - x[:-1] - L0)  # on left mass from right
        a = np.zeros(n)
        a[:-1] += pair_force
        a[1:] -= pair_force
        a -= b_damp * v
        v = v + dt_chain * a
        x = x + dt_chain * v
        xs[t + 1], vs[t + 1], acc[t] = x, v, a
    return xs, vs, acc

We generate training trajectories with random initial conditions. The regression below sees only the two end masses, where the net spring force equals one pairwise force. This restriction exposes an identifiability property rather than taking a shortcut.

In[12]:

Code

n_train_traj, n_masses = 40, 8
train_data = []
for _ in range(n_train_traj):
    x0 = np.linspace(0, (n_masses - 1) * L0, n_masses) + rng.normal(
        0, 0.15, n_masses
    )
    v0 = rng.normal(0, 0.3, n_masses)
    train_data.append(simulate_chain(x0, v0))

Here is the subtle point. At a single instant, the acceleration equation for an interior mass $i$ determines only the difference $\varphi(d_{\text{right}}) - \varphi(d_{\text{left}})$ of the two pairwise forces. Any component of $\varphi$ that is constant in $d$ cancels from those instantaneous interior equations. If observations contain only such equations, any law of the form $\varphi(d) + c$ fits them equally well: the constant is a gauge freedom. Boundary equations break this symmetry because an end mass has only one neighboring force; a periodic chain with only pairwise differences retains it. Longer trajectories do not remove the offset by themselves when every observed acceleration still depends only on force differences, although additional boundary, constitutive, or force measurements can. Learning interaction laws from motion capture often has flavors of this problem: collective motion underdetermines individual forces unless boundary observations, contact events, direct force measurements, or a prior pins them down. Full GNS models bypass explicit force-law identification by supervising accelerations from relative inputs, but that does not make latent messages or physical forces uniquely identifiable. The underdetermination belongs to the inverse problem and can persist inside the architecture.

In[13]:

Code

D_list, F_list = [], []
for xs, vs, acc in train_data:
    n_rows = acc.shape[0] - 1  # terminal acceleration was not evaluated
    d0 = xs[:-1, 1] - xs[:-1, 0] + rng.normal(0, 0.01, n_rows)
    f0 = m * (acc[:-1, 0] + b_damp * vs[:-1, 0])  # phi(d0)
    dL = xs[:-1, -1] - xs[:-1, -2] + rng.normal(0, 0.01, n_rows)
    fL = m * (acc[:-1, -1] + b_damp * vs[:-1, -1])  # -phi(dL)
    D_list += [d0, dL]
    F_list += [f0, -fL]

D = np.concatenate(D_list)[:, None]
F = np.concatenate(F_list)
coef = np.linalg.lstsq(np.column_stack([D, np.ones_like(D)]), F, rcond=None)[0]
w1, w0 = float(coef[0]), float(coef[1])

Out[14]:

Console

Learned interaction law: phi(d) = 1.992 * d + -1.992
True interaction law:    phi(d) = 2.000 * d + -2.000

The regression recovers both the slope (stiffness) and the intercept (which encodes the rest length) from end-mass observations, up to the injected separation noise. Here $d$ is the mass separation, so $\varphi(d)=k(d-L_0)$ vanishes at $d=L_0$ and has intercept $-kL_0$ at $d=0$ , a fully compressed ideal spring. The learned law absorbs that offset automatically. Now the transfer test. We roll the learned pairwise law forward with the same semi-implicit integrator on a held-out chain of 8 masses and on a chain of 12 masses, using no size-specific learned parameters.

In[15]:

Code

def rollout_learned(x0, v0, w1, w0, steps=steps_chain):
    """Integrate the chain using the learned pairwise law phi(d) = w1*d + w0."""
    n = len(x0)
    x, v = x0.copy(), v0.copy()
    xs = np.zeros((steps + 1, n))
    xs[0] = x
    for t in range(steps):
        d = x[1:] - x[:-1]
        f_pair = w1 * d + w0
        a = np.zeros(n)
        a[:-1] += f_pair
        a[1:] -= f_pair
        a -= b_damp * v
        v = v + dt_chain * a
        x = x + dt_chain * v
        xs[t + 1] = x
    return xs


x0a = np.linspace(0, 7 * L0, 8) + rng.normal(0, 0.15, 8)
v0a = rng.normal(0, 0.3, 8)
xs_true8, _, _ = simulate_chain(x0a, v0a)
xs_pred8 = rollout_learned(x0a, v0a, w1, w0)
rms8 = np.sqrt(np.mean((xs_pred8 - xs_true8) ** 2, axis=1))

x0b = np.linspace(0, 11 * L0, 12) + rng.normal(0, 0.15, 12)
v0b = rng.normal(0, 0.3, 12)
xs_true12, _, _ = simulate_chain(x0b, v0b)
xs_pred12 = rollout_learned(x0b, v0b, w1, w0)
rms12 = np.sqrt(np.mean((xs_pred12 - xs_true12) ** 2, axis=1))

snap_t = 200  # snapshot time index for the trajectory comparison

Out[16]:

Visualization

Line plot of mass positions with truth and learned prediction overlapping closely. — Mass positions of a held-out 8-mass chain at the snapshot time: ground truth versus rollout under the learned pairwise interaction law. The learned law tracks the true positions closely across all eight masses.

Two error-versus-time curves for 8 and 12 mass chains rising slowly from near zero. — Rollout error versus horizon for held-out chains of 8 and 12 masses. The model trained only on 8-mass chains transfers to 12 masses with errors of the same order.

The learned law rolls out with tiny error on both tested chain sizes. Because the same local law applies at both sizes, data from 8-mass chains transfers to this matched 12-mass case. Transfer to other sizes still depends on coverage of the relevant local configurations and dynamics. The error curves also illustrate compounding rollout error: small per-step mismatches accumulate over the horizon. Full GNS models retain shared local updates while replacing this linear regression with neural message passing. Their learned interactions remain trustworthy only where training data or structural constraints support them.