Time Without Timesteps — Liyu Zerihun

Research paper with Mark Shinyoung Lee.

Read the paper

Overview

Almost every simulation of a differential equation is a march. You have the state now, an integrator estimates the state a moment later, and repeating that gives you a trajectory. It works, it is mature, and it is usually the right tool.

But the march is a choice about how to construct the answer, not a fact about the answer. A differential equation with its initial conditions picks out one trajectory, the one that satisfies the dynamics and every coupling constraint at once. The trajectory is a global object. The order a solver builds it in is incidental.

Time Without Timesteps, or TWT, takes that seriously for coupled systems. Instead of stepping every subsystem forward together, we train a neural surrogate per subsystem type that maps a whole driving input trajectory and an initial condition to a whole output trajectory. A coupled system is then a fixed point: every subsystem's trajectory has to equal what its surrogate predicts under the inputs its neighbors produce. Simulation becomes a self-consistency problem over complete trajectories.

This is not a new idea. It is waveform relaxation, from 1982 circuit simulation, with the inner numerical integration swapped for a learned operator. What we found is that the swap changes which solvers are practical, which in turn changes what the method can do.

The Circularity

Take K subsystems. Each one evolves under its own dynamics, driven by an input that neighboring subsystems produce:

x_i' = f_i(x_i, u_i(t), theta_i),    x_i(0) = x_i0
u_i(t) = g_i({x_j(t) : j in N(i)})

If every u_i were known, every subsystem could be solved on its own. The trouble is that the drives are functions of the trajectories you are trying to find.

Timestepping resolves this locally: advance everything one step, and the circularity never gets a chance to bite. Waveform relaxation resolves it globally: guess every trajectory, integrate each subsystem under the inputs those guesses induce, repeat until the guesses agree.

The classical method never became a general tool, and cost is a big reason. Every relaxation sweep is a full simulation of every subsystem, so each iteration costs as much as the whole problem, and the derivative of a sweep costs the same again.

Learned Module Operators

For each module type we train a surrogate S that approximates the isolated solution operator, whole trajectory in, whole trajectory out:

x_i_hat = S_phi(theta_i, u_i, x_i0)

Here u_i is a discretized input trajectory and x_i_hat is the full predicted output. The network is a residual MLP; one forward pass gives the entire trajectory. Time is still in the representation. For the oscillator system we work in a Chebyshev coefficient basis, 256 coefficients standing in for 1500 samples. The neuron system uses raw samples. The representation is an implementation detail; the formulation does not care.

Training data comes from simulating isolated modules with a conventional solver under known drives. Since the drives a module sees inside a coupled solve are not known in advance, we use a rolling buffer: early samples use simple drives, later samples draw their drives from weighted combinations of previously generated outputs. This turned out to matter more than any architectural decision, and it is also where most of our mistakes came from. More on that below.

The Solve

With trained surrogates, the coupled solution is a fixed point over the stack of all trajectories X:

X = F(X),    F(X)_i = S_phi(theta_i, g_i({X_j : j in N(i)}), x_i0)

At a zero-residual fixed point, every trajectory equals the one its own surrogate predicts under the inputs its neighbors induce.

The natural iteration is the classical one, Picard: X <- F(X). It converges only while the coupled operator is contractive,

rho(J_F) < 1,    J_F = (dS_phi/du) W

where W is the coupling operator and rho is the spectral radius. For the all-to-all diffusive coupling we use, ||W|| is exactly the coupling strength k_c. This is a property of the iteration, not of the problem.

So we also solve the fixed point with Newton's method on the residual G(X) = X - F(X):

(I - J_F) dX = -G(X),    X <- X + dX

The Jacobian is never formed. GMRES solves the linear system, and each GMRES iteration needs one product (I - J_F) v, which with a learned operator is a single forward-mode pass. Newton converges wherever I - J_F is nonsingular, no contraction needed.

This is the part that only works because the operator is learned. We ran the same Newton solve with an exact RK4 integrator in the module slot, and at moderate coupling it took 1548 seconds against 13.7 seconds for plain relaxation, because every Krylov iteration was a full numerical integration. The equilibrium view and the learned operator are not independent choices. Each one is what makes the other practical.

On 20 coupled van der Pol oscillators with all-to-all diffusive coupling, the coupled solve converges in 4 to 10 Newton iterations on trajectories that span 1500 timesteps. At weak to moderate coupling the trajectory error against ground truth is 1 to 6 percent. The same surrogate with coupling disabled is wrong by more than 100 percent everywhere.

Contraction Predicts Convergence

The result I find most useful is that you can predict where the coupled solve will fail before running it.

Power iteration on the isolated learned operator gives rho(dS_phi/du), with no coupled system involved. Since ||W|| = k_c, the contraction condition predicts a boundary:

k_c* = 1 / rho(dS_phi/du)

We measured rho across stiffness and overlaid the predicted boundary on the observed sweep.

Stiffness mu	rho (measured)	Predicted k_c*	Converged at k_c = 2
0.5	0.87	1.15	4/5
1.0	0.64	1.56	3/5
2.0	0.57	1.75	5/5
3.0	0.54	1.85	5/5
5.0	0.52	1.92	5/5

The ordering is exact. Iteration count, convergence rate, and trajectory error all follow rho.

Two things fall out of this. First, the true operator at mu = 1 has rho around 0.54, and the learned one has 0.64. Learned operators are more expansive than the dynamics they approximate, so the classical iteration gives up earlier than the physics would. Second, and more surprising, accuracy and composability are different properties that do not improve together. The surrogate is least accurate in isolation at mu = 5 and composes best there. Across three separately trained models, better fit came with larger rho every time.

The obvious fix, penalizing the Jacobian norm during training, made things worse. A random directional probe in n dimensions puts 1/n of its energy along any one direction, so what it estimates is the Frobenius norm, not the spectral radius. Constraining it crushed every direction except the one that mattered. Power iteration, not random probing. We report that as a negative result because it cost us time and would cost someone else time too.

Gradients Without a Backward March

Differentiable simulators exist, and gradients through coupled ODEs are routine. But both standard routes, unrolling the integrator or solving the continuous adjoint, inherit the sequential structure of the forward march. Unrolling stores every intermediate state. The adjoint sweeps backward across the full horizon.

At a fixed point the gradient is different in kind. For an objective J evaluated at the converged X*, the implicit function theorem gives

(I - J_F^T) v = dJ/dX*
dJ/dtheta = v^T (dF/dtheta)

There is no time recursion in this. It is one linear system, and it does not depend on how the forward solve was performed. We solve it with GMRES, at memory independent of how many Newton iterations the forward solve took.

The standard equilibrium-model backward pass expands (I - J_F^T)^-1 as a Neumann series. That reintroduces rho(J_F) < 1, so it fails exactly where Newton was brought in to succeed. GMRES does not have that problem.

	Coupling 0.1, N = 6	Coupling 0.1, N = 50	Coupling 0.5, N = 6
Cosine sim. vs. unrolled	1.0000	1.0000	-0.13
Rel. err. vs. finite diff.	0.06%		0.04%
Neumann adjoint	converges	converges	diverges
Unrolled memory, depth 2 to 160	102 to 146 MB	108 to 472 MB
Implicit memory	112 MB	169 MB	112 MB

Where the comparison is meaningful, implicit and unrolled gradients agree to five digits. Past the contraction boundary, unrolled backpropagation diverges by a factor of 83 and the Neumann adjoint returns non-finite values. The implicit gradient stays correct to 0.04 percent against finite differences. That is the regime the Newton solver was introduced to reach, and the implicit gradient is the only one that survives it.

To check the gradients do real work, we recovered the coupling weights of a six-oscillator ring from trajectories generated by the true physics. Parameter error fell from 83 percent to a median of 1.2 percent.

Neurons

The second system is a pair of Hodgkin-Huxley neurons with a recurrent synapse, where the surrogate predicts only the coupling variable s(t) and keeps voltage and gating latent. The coupling here is nonlinear, since synaptic current depends on both pre- and postsynaptic state.

The cleanest statement is in the activity numbers. A feedforward model predicts by construction that recurrence does nothing: its activity is pinned at 0.0551 no matter the recurrent weight. The true activity falls monotonically as recurrence strengthens, because the synaptic reversal potential makes strong input shunting during a spike. The fixed point tracks that dependence. The uncoupled view cannot.

The error numbers are less flattering, and we report them as they are. Self-consistency roughly halves feedforward error, but the exact-operator floor is under 2 percent, so most of what remains is the surrogate, amplified through a worse-conditioned fixed point. At the strongest recurrence the solve does not converge.

What the Method Costs

Throughout the paper we report three numbers rather than one: the self-consistency residual ||X - F(X)||, the trajectory error against ground truth, and the error of the same fixed point solved with an exact integrator as the module operator. The third is the floor attributable to the formulation. The gap between it and the surrogate result is what learning costs.

That decomposition paid for itself. Early on, the module simulator applied the drive at the end of each step and the reference applied it at the start. Both are valid O(dt) discretizations, but they define different fixed points, and the discrepancy grew from under 1 percent to over 10 percent with coupling. It looked exactly like surrogate error until the exact-operator check isolated it.

The residual also needs a caveat. It measures agreement under the learned surrogate, not satisfaction of the original equation. A residual of 1e-12 can sit next to 30 percent physical error, and at strong coupling it does.

Limitations

The systems are small: 20 oscillators, neuron pairs. Nothing in the formulation is tied to that scale, but the structural claims about sequential depth would be a lot more convincing at a thousand modules.

We report sequential depth, not wall-clock speedup. Neither our implementation nor the reference integrator is optimized, so a timing comparison would not mean anything yet. The claim is that iteration count is decoupled from trajectory length and that iterations are parallel across modules. The practical claim needs an implementation we have not built.

Surrogate error of 1 to 4 percent is nowhere near solver-grade, and the fixed point amplifies it. A converged fixed point is only as meaningful as the operators it is built from.

The dominant failure in practice was training coverage. Twice, an experiment probed a parameter range the data generator had never sampled, and twice we blamed the method before tracing it to the data. Composing learned operators makes this worse than ordinary surrogate use, because the fixed-point iteration decides what inputs the operator sees, and those inputs are not known ahead of time. Generating data that covers the distribution coupling induces, rather than one chosen beforehand, is the central unsolved problem for this class of method.

Direction

Windowing the time axis is the classical remedy for slow relaxation convergence, and it would reduce rho directly. We did not pursue it and it is the obvious next step.

Beyond that: learned preconditioners for the fixed-point solve, compressed trajectory representations for long horizons, and training objectives that control the spectral radius directly instead of through the proxy that failed.

The larger point is the one about composability. Standard surrogate evaluation measures how well an operator fits in isolation. For anything that composes learned operators, the norm of the Jacobian looks at least as predictive of end-to-end behavior, and it can be measured before any composition is attempted. That seems like it should generalize past this paper.