Numerical simulation of dynamical systems is usually organized as a causal march through time: each state is computed from the previous one. We explore a different formulation for coupled systems. For each subsystem type we train a neural surrogate mapping a full driving trajectory and initial condition directly to a full output trajectory; following classical waveform relaxation, coupled systems are assembled by enforcing self-consistency among these trajectories: simulation becomes a fixed-point problem over complete trajectories rather than a stepwise rollout. On coupled van der Pol oscillators and Hodgkin-Huxley neuron networks, sequential depth becomes the number of solver iterations: 4-10 Newton iterations where the reference integrator takes 1500 steps. The gradient likewise loses its time recursion: it becomes a linear system solved by GMRES at memory independent of solver depth. A single scalar measured from the learned operator, the spectral radius of its Jacobian, predicts in advance where the coupled solve will converge; past that boundary, unrolled backpropagation diverges and a Neumann adjoint fails, while the implicit gradient remains correct to 0.04%. We report where the approach succeeds and where surrogate error degrades it.
Zermelo's algorithm is a classical method for computing the maximum likelihood estimator in the Bradley--Terry (BT) model, but its convergence can be slow in practice. To accelerate computation, Newman introduced a family of Zermelo-type fixed-point iterations parameterized by $α$, with Zermelo's algorithm recovered at $α=1$. Empirical evidence suggests that the choice $α=0$ often converges substantially faster, making it a promising alternative, yet the mechanism underlying this acceleration remains elusive. This paper provides theoretical insight into this phenomenon through a systematic local convergence analysis. We derive closed-form expressions for local convergence factors under synchronous and asynchronous updates and analyze their dependence on $α$ via spectral analysis of the associated Jacobian matrices. For synchronous updates, we show that the algorithm may fail to converge when $α<1$, and its local convergence factor is quasi-convex in $α$ under the population BT model. In contrast, asynchronous updates are always locally convergent, and their local convergence factor is provably monotonically increasing in $α$ under the population BT model of consistently ordered bipartite comparison graphs, establishing the optimality of $α=0$ in this setting. We further establish asymptotic approximation results for the population convergence factors under the BT model, justifying their practical relevance. Numerical experiments on synthetic and real-world datasets confirm the theory. Our analysis complements existing convergence results and shows that the acceleration of $α=0$ arises not only from the parameter choice but, more importantly, from the use of asynchronous updates.
Non-expansive two-time-scale stochastic approximation is governed by a slow stochastic Krasnoselskii--Mann fixed-point iteration rather than by contraction to a unique equilibrium. We study this regime under a contractive fast map and a non-expansive reduced slow map. We first prove a finite-horizon lower bound showing that, for any prescribed slow stepsize schedule $(β_k)$, the classical KM residual scale $(\sum_{i<N}β_i(1-β_i))^{-1}$ is worst-case sharp for the corresponding unregularized KM update. Combined with the raw fast-tracking leakage scale, this explains the previously observed $k^{-1/4+o(1)}$ last-iterate mean-square residual exponent. We then introduce a residual-preconditioned slow oracle that cancels the first-order dependence on the fast tracking error. In a nested Tikhonov-KM algorithm, the uncorrected oracle yields total-sample rate $T^{-1/4+o(1)}$, while the corrected oracle yields $T^{-1/3+o(1)}$. This improvement comes from changing the slow-oracle bias from first order to second order in the fast error after all inner-loop samples are counted. Finally, we show that the repeated inner-loop cost of the nested method can be avoided in a smooth derivative-oracle model. A single-loop algorithm that tracks both the fast equilibrium and the leakage preconditioner online achieves $T^{-1/2+o(1)}$ with $O(1)$ primitive samples per iteration.
Multidisciplinary design analysis of coupled engineering systems requires the computation of equilibrium states in which all disciplinary coupling variables are mutually consistent. Conventional fixed-point iteration resolves this consistency problem separately at each design point, which can become expensive when disciplinary evaluations are costly and many analyses are required in outer-loop tasks such as multidisciplinary design optimization, uncertainty quantification, or digital twin updating. This paper introduces REMAL, a residual manifold surrogate modeling framework for coupled systems. Instead of approximating each discipline independently or directly learning converged coupling variables, the proposed method learns a surrogate model of the joint residual manifold via multitask Gaussian process models. An entropy-based active learning strategy selects additional residual evaluations near uncertain zero-contour regions, and equilibrium states for new design inputs are recovered by solving a nonlinear least squares optimization problem using only the trained surrogate. The method is evaluated on four engineering coupled system benchmarks: a satellite model, an aerostructural model, a finite-element gas-turbine heat-transfer and economics model, and a modified turbine model with added feedback coupling. Across these cases, REMAL consistently demonstrates the cost effectiveness when repeated evaluations of the fixed point across the design space are necessary. Theoretically, we show that, under mild assumptions, REMAL's predictive fixed point error is bounded.