--- title: Adaptive Pseudo-Time Stepping url: https://www.emergentmind.com/topics/adaptive-pseudo-time-stepping type: topic --- # Adaptive Pseudo-Time Stepping Searching arXiv for recent and foundational papers on adaptive pseudo-time stepping and closely related adaptive time-stepping formulations. Adaptive pseudo-time stepping denotes a class of numerical and algorithmic strategies in which an evolution variable that is not necessarily physical time is advanced with nonuniform increments chosen to balance stability, robustness, accuracy, and computational cost. Across the literature, the “pseudo-time” interpretation appears in several distinct but structurally related settings: residual networks viewed as explicit ODE discretizations along depth, geometric time reparametrizations of Hamiltonian or particle systems, adaptive multistep integration for dissipative gradient flows, explicit stabilization strategies for stiff ODEs and SPDEs, and goal-oriented timestep control for weakly instationary hyperbolic PDEs. A unifying theme is that the step size is treated as a control variable for the evolution itself rather than as a passive discretization parameter. In residual networks, this control acts along network depth and is learned jointly with model parameters [1911.10305]. In geometric integration, it appears through a reparametrization \(dt/d\tau=g(z)\), so that a constant step in pseudo-time \(\tau\) induces a variable step in physical time \(t\) [1108.0322, 1211.3664]. In PDEs and rate-independent systems, it is tied either to error indicators, energy laws, or complementarity residuals [2008.00212, 2204.05860, 1405.3497]. ## 1. Depth as pseudo-time in residual and neural-ODE-like models A direct pseudo-time interpretation arises in residual architectures when a residual block is written as \[ \mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1, \] which is exactly the forward Euler discretization of \[ \dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T], \] with \(t_j=jh\), \(T=Dh\), and \(\mathbf{W}(t_j)\approx \mathbf{W}_j\) [1911.10305]. In this interpretation, \(t\) is not physical time but a pseudo-time along depth, and \(\mathbf{y}(t)\) is the feature representation as it evolves through the network. Adaptive pseudo-time stepping in this setting replaces the fixed step \(h\) by learned, layer-dependent increments \(\Delta t_j\): \[ \mathbf{y}_{j+1} = \mathbf{y}_j + \mathcal{F}(\mathbf{y}_j,\mathbf{W}_j)\,\Delta t_j, \qquad T=\sum_{j=0}^{D-1}\Delta t_j. \] The pseudo-time grid and the effective evolution time are therefore learned jointly with the network weights [1911.10305]. The motivation is explicitly compared with adaptive Runge–Kutta–Fehlberg (RKF), which updates the next step by \[ \Delta t_{j+1} = k \,\Delta t_j \left( \frac{Tol}{|y_{j+1}-\hat y_{j+1}|} \right)^{1/(p+1)}, \] with \(p=4\), but the network setting does not estimate local truncation error explicitly [1911.10305]. The central stability bound established for variable-step ResNets is \[ \|\mathbf{y}_D^\epsilon - \mathbf{y}_D\| \le \epsilon \cdot \prod_{j=0}^{D-1}\bigl(1+\|\mathbf{W}_j\|_2\,\Delta t_j\bigr), \] where \(\|\mathbf{W}_j\|_2\) is the spectral norm of the weight matrix in block \(j\) [1911.10305]. This identifies the product “step size \(\times\) spectral norm” as the relevant local amplification mechanism. The same work also shows that very small \(\Delta t_j\) can make layerwise gradients nearly identical, \[ \frac{\partial L}{\partial \mathbf{y}_n}\approx \frac{\partial L}{\partial \mathbf{y}_D}, \] thereby smoothing the feature evolution too strongly and reducing representational power [1911.10305]. This yields the characteristic trade-off: small pseudo-time steps improve stability and robustness, while larger steps improve discriminative transformation but can destabilize training. The implemented controller is a data-independent, parameter-driven LSTM. At each residual block, the convolution kernels \(\mathbf{w}_d\) are compressed, mapped through an input transform, processed recurrently, and converted into a channel-wise step-size vector \[ \Delta \mathbf{t}_d = g_{out}(W_{out}\mathbf{h}_d+\mathbf{b}_{out}), \] with sigmoid output so that each channel’s step lies in \((0,1)\) [1911.10305]. The residual update becomes \[ \mathbf{y}_{d+1} = \mathbf{y}_d + \mathcal{F}(\mathbf{y}_d,\mathbf{w}_d)\otimes \Delta \mathbf{t}_d. \] Because the controller is data-independent, it is discarded after training and only the learned step-size vectors are retained at inference, giving “no additional overhead in inference phase” [1911.10305]. On ImageNet, the method improves top-1 accuracy across ResNet, ResNeXt, and SENet variants while leaving inference GFLOPs essentially unchanged; for example, on ResNet-50 the baseline has 25.56M parameters and 3.86 GFLOPs, while TSC\(_{\text{LSTM}}\)-ResNet-50 has 27.83M parameters and 3.89 GFLOPs during training, but 25.57M parameters and 3.86 GFLOPs at inference [1911.10305]. The learned step-size profiles are small in earlier layers and larger in late layers, suggesting a staged pseudo-time allocation rather than a uniform depth discretization [1911.10305]. ## 2. Geometric reparametrization and structure-preserving pseudo-time A second major line of work interprets adaptive stepping as a time reparametrization. For Hamiltonian systems, naive variable-step symplectic integration typically fails for two distinct reasons. If the step depends on time alone, \(\Delta=\Delta(t)\), backward error analysis shows that the method can remain symplectic while developing parametric instabilities caused by resonance between oscillations of \(\Delta(t)\) and orbital motion [1108.0322]. If the step depends on phase-space position, \(\Delta=\Delta(q,p)\), the transformation \[ dt=\Delta(q,p)\,d\tau \] leads to equations in pseudo-time \(\tau\) that are generally no longer in canonical Hamiltonian form [1108.0322]. The proposed remedy is to integrate with constant step in \(\tau\) while preserving the transformed geometric structure either by an extended phase-space method or by a non-canonical mixed-variable generating function [1108.0322]. In the extended phase-space construction, time is promoted to a coordinate \(q_0=t\) with conjugate momentum \(p_0\), and an extended Hamiltonian \[ K(q,p,q_0,p_0)=\frac{1}{\rho(\mathbf z)}\bigl(H(q,p,q_0)+p_0\bigr) \] is introduced, with \(d\tau=\rho(\mathbf z)\,dt\) [1108.0322]. Constant pseudo-time steps \(h=\Delta\tau\) then induce adaptive physical steps \(dt=h/\rho\). The alternative non-canonical method directly preserves the two-form \(\rho(x,y)\,dx\wedge dy\) in one degree of freedom by constructing a Poisson integrator via generating functions [1108.0322]. In both cases, pseudo-time stepping is not an auxiliary implementation trick but the mathematically correct mechanism by which adaptivity is reconciled with geometric structure. The same reparametrization principle is used in accelerator simulation with the Boris–Buneman integrator [1211.3664]. There, adaptive pseudo-time stepping is realized through a Sundman transformation \[ \frac{dt}{d\tau}=g(z), \] with \(z=(\mathbf{x},\mathbf{p})\), so that a constant step \(\Delta\tau\) in pseudo-time gives a variable physical step \(\Delta t=g(z)\Delta\tau\) [1211.3664]. The choice of \(g\) is driven not by a local truncation error estimate but by a physically motivated state function proportional to a power of the maximum self-field acceleration: \[ g(\mathbf{p}_{1..N},\mathbf{f}_{1..N}^{\text{self}})\propto \left(\max_i\|\mathbf{a}_i\|\right)^{-\beta/2}, \] with \(\beta=1\), hence \(g\propto A_{\max}^{-1/2}\) [1211.3664]. Strong space charge therefore produces small physical steps, while weak space charge produces large steps. This is embedded into an adaptive multiple-time-stepping Boris–Buneman method in which the outer step \(h=\lambda\Delta\tau\) controls the expensive self-field updates, while an inner step \(h/m\) remains close to a prescribed \(\Delta t_{\text{inner}}\) for external fields [1211.3664]. In a photoinjector test, the adaptive method achieves a given emittance error with dramatically fewer self-field solves than a fixed-step MTS method; in a cyclotron test, where space charge varies only mildly, the gain disappears, indicating that the usefulness of adaptive pseudo-time depends strongly on scale separation and force variability [1211.3664]. These geometric and force-based formulations suggest a broad interpretation: pseudo-time stepping is a controlled reparametrization of evolution. The step size is encoded in \(dt/d\tau=g(z)\), and the numerical problem becomes one of integrating the \(\tau\)-system with a structure-preserving method. ## 3. Stability-driven adaptivity in stiff deterministic and stochastic systems A third strand uses adaptive pseudo-time steps to stabilize explicit or semi-explicit discretizations in stiff settings. In stiff ODEs, the basic difficulty is the classical explicit stability restriction. For \(\dot u+\lambda u=0\), explicit Euler is stable only if \[ 0\le k_n\lambda\le 2. \] The proposed explicit stabilization strategy is to take a large unstable step \(K\), followed by \(m\) small stable steps \(k\), so that the composite amplification \[ (1-k\lambda)^m(1-K\lambda) \] remains stable [1205.2751]. For \(K\lambda\gg1\), the number of stabilizing small steps scales only logarithmically, \[ m \approx \frac{\log(K\lambda)}{c}, \quad c:=k\lambda, \] which implies that a small amount of local pseudo-time refinement can stabilize a much larger effective step [1205.2751]. For systems with a spectral gap this yields large speedups; for parabolic spectra, Chebyshev and dyadic damping strategies extend the real-axis stability interval with cost reductions of order \((K\lambda_N)^{-1/2}\) and \((K\lambda_N)^{-1/3}\), respectively [1205.2751]. This suggests a pseudo-time interpretation in which large outer steps advance the slow dynamics and short bursts of small steps quench stiff modes. In SDEs with multiplicative, non-globally Lipschitz coefficients, \[ dX(t)=X(t)f(X(t))\,dt + X(t)g(X(t))\,dW(t), \] the adaptive rule \[ h_n = \frac{\bar h}{1+\lfloor f(X_n)\rfloor + \lfloor g^2(X_n)\rfloor} \] is designed so that \[ h_n f(X_n)\le \bar h,\qquad \sqrt{h_n}g(X_n)\le \sqrt{\bar h}, \] thereby uniformly bounding the drift and diffusion responses per step [1706.03098]. The analysis is delicate because the Brownian increments over random steps are no longer independent, so the step times are constructed as stopping times and the proof uses a discrete Itô formula with conditional moments together with martingale convergence lemmas [1706.03098]. Under \(\sup_{u\neq0}2f(u)/g^2(u)=\beta<1\), the adaptive Euler–Maruyama scheme satisfies \[ \lim_{n\to\infty}X_n=0 \quad\text{a.s.} \] for sufficiently small \(\bar h\), while if \(\liminf_{u\to0}2f(u)/g^2(u)>1\), then \(\mathbb P[\lim X_n=0]=0\) [1706.03098]. The same rule also yields positivity with arbitrarily high probability on any finite horizon [1706.03098]. Here pseudo-time stepping is not guided by a classical error estimate but by a Lyapunov-driven control of qualitative dynamical properties. The corresponding SPDE literature uses state-dependent, admissible time-step rules to control explicit exponential integrators in the presence of one-sided Lipschitz but non-globally Lipschitz drift [1812.09036]. For semilinear SPDEs \[ dX = [-AX+F(X)]\,dt + B(X)\,dW, \] the adaptive explicit exponential schemes \[ X_h^{n+1} = e^{-\Delta t_{n+1}A_h}\bigl(X_h^n + P_hB(X_h^n)P_J\Delta W_{n+1}\bigr) + \Delta t_{n+1}\varphi_1(\Delta t_{n+1}A_h)P_hF(X_h^n) \] and \[ X_h^{n+1} = e^{-\Delta t_{n+1}A_h}\bigl(X_h^n + P_hF(X_h^n)\Delta t_{n+1} + P_hB(X_h^n)P_J\Delta W_{n+1}\bigr) \] use adaptive \(\Delta t_{n+1}\) chosen so that along the discrete trajectory \[ \|F(X_h^n)\|^2 \le R_1 + R_2\|X_h^n\|^2 \] holds [1812.09036]. Admissible step controls include \[ \Delta t_{n+1}\le \frac{\delta}{\|F(X_h^n)\|^{1/\theta}}, \qquad \Delta t_{n+1}\le \delta\,\frac{\|X_h^n\|}{\|F(X_h^n)\|}, \qquad \Delta t_{n+1}\le \frac{\delta}{\|X_h^n\|}, \] among others [1812.09036]. Under these rules, the full discretization achieves strong convergence \[ \|X(T)-X_h^N\|_{L^2(\Omega;H)} \le C_0\Big(h^{1+r-\epsilon} + \lambda_{J+1}^{-(1+r)/2+\epsilon} + \Delta t_{\max}^{1/2-\epsilon}\Big), \] which is essentially the same temporal order as in the globally Lipschitz case [1812.09036]. Numerically, these adaptive schemes outperform fixed-step tamed or stopped alternatives at equal computational effort in Allen–Cahn and Swift–Hohenberg SPDEs [1812.09036]. A closely related but distinct phenomenon appears in the Poisson–Nernst–Planck equations with generalized Frumkin–Butler–Volmer boundary conditions [1703.10297]. The paper develops adaptive VSBDF2 and VSSBDF2 schemes using a coarse–fine local truncation error estimator. For the semi-implicit scheme, when the underlying dynamics approach a steady state, the adaptive time step does not grow indefinitely but stabilizes at a limiting size \(dt_\infty\), and the numerical solution only “nearly” converges to the steady state [1703.10297]. The paper states that in the companion article the linearized SBDF2 scheme about the steady-state solution is shown to be conditionally stable, and that this conditional stability is the cause of the adaptive time-stepper’s behaviour [1703.10297]. By contrast, the fully implicit BDF2 adaptive scheme is not subject to such time-step constraints, but the required nonlinear solve makes it significantly slower [1703.10297]. This suggests that in pseudo-time marching to steady states, adaptivity alone does not remove the fundamental stability ceiling of a conditionally stable base integrator. ## 4. Energy-dissipative pseudo-time stepping in gradient flows and rate-independent systems In dissipative PDEs, adaptive pseudo-time stepping is often governed by discrete energy laws rather than by classical stability regions. For the phase field crystal (PFC) equation, \[ \partial_t\Phi = \Delta\bigl(\Phi^3-\epsilon\Phi+(1+\Delta)^2\Phi\bigr), \] the adaptive variable-step BDF2 discretization is \[ D_2\phi^n = \Delta_h\mu^n,\qquad \mu^n = (1+\Delta_h)^2\phi^n + (\phi^n)^3 - \epsilon\phi^n, \] with variable-step BDF2 operator \[ D_2 v^n = \sum_{k=1}^n b_{n-k}^{(n)}\delta v^k \] and step ratios \(r_k=\tau_k/\tau_{k-1}\) [2008.00212]. The key structural condition is the zero-stability bound \[ 0t_{k-1}\) in sticking or rate-independent slip [2204.05860]. Adaptivity is driven by two residuals: one for local stability/complementarity, \[ I_1^k = \int_{s_{k-1}}^{s_k} \hat t'(s)\,{}_{\mathbb V^*}\{-D_z\mathcal I(\hat t(s),\hat z(s)),\partial\mathcal R(0)\}\,ds, \] and one for the energy-dissipation balance, \(I_2^k\), built from the discrete energy identity [2204.05860]. Steps are accepted if both residuals are below a tolerance, doubled if both are below half the tolerance, and halved otherwise [2204.05860]. The resulting scheme increases the pseudo-time step during sticking and in viscous jumps, and numerical experiments show a significant reduction in the number of steps relative to a uniform pseudo-time discretization [2204.05860]. This is one of the clearest cases where “pseudo-time stepping” is literal: the step size is attached to an auxiliary curve parameter rather than to physical time. A recent filtered implicit framework also fits this energy-stability-oriented perspective. A variable-step pre-filter/post-filter method based on Implicit Euler constructs a second/third-order embedded pair, \[ \tilde y_n = y_n - \frac{\alpha_n}{2}\kappa_{n-1},\qquad \frac{y_{n+1}-\tilde y_n}{k_n}=f(t_{n+1},y_{n+1}),\qquad y_{n+1}^{3\text{rd}} = y_{n+1}-\beta_n(\kappa_n-\kappa_{n-1}), \] with \[ \alpha_n = \frac{k_n^2}{k_{n-1}k_{n-2}}, \] and an explicit variable-step formula for \(\beta_n\) [2504.00405]. The difference between the second- and third-order outputs gives the estimator \[ EST = |y_{n+1}^{3\text{rd}}-y_{n+1}^{2\text{nd}}|, \] and a simple halving/doubling controller adjusts \(k_n\) [2504.00405]. The method is proposed in the standard ODE initial-value setting but is described as directly relevant to pseudo-time stepping toward steady states because it retains the Implicit Euler solve while allowing larger, adaptively selected steps once the solution changes slowly [2504.00405]. ## 5. Goal-oriented and locality-aware adaptivity in hyperbolic and particle systems For weakly instationary hyperbolic conservation laws, pseudo-time adaptivity is often best understood as goal-oriented global time-step control. In scalar conservation laws and the 2D Euler equations, the error in a target functional is represented as \[ J(u)-J(u_h) = (R(u_h),\varphi-\varphi_h)_\Omega, \] where \(R(u_h)\) is the residual and \(\varphi\) solves the adjoint problem [1405.3497, 1405.6510]. Using space–time \(L^2\) projections, the error splits into a temporal part \(\eta_k\) and a spatial part \(\eta_h\): \[ J(u)-J(u_h)=\eta_k+\eta_h \] up to boundary and linearization terms [1405.3497, 1405.6510]. The main innovation is that, because only \(\varphi-\varphi_h\) matters by Galerkin orthogonality, it is sufficient to compute the dual gradient \(w=\nabla\varphi\), which satisfies a conservation law instead of a transport equation: \[ \partial_t w + \partial_x(a\,w)=\partial_x\psi \] in 1D, or more generally a conservative dual PDE in multiple dimensions [1405.6510, 1404.4503, 1405.3497]. This allows the same finite-volume or DG machinery as in the primal solver. The time indicator \(\eta_k^n\) is then used to design new step sizes by an equidistribution principle. In the scalar case, the update is written as \[ k_m := k_n\, \frac{T_k(L+1)/T}{\bar\eta_k^n}, \] so that stationary regions, where \(\bar\eta_k^n\) is tiny, receive very large timesteps, while transient regions receive small ones [1405.6510]. In the Euler bump problem, the combined space–time adaptive method takes timesteps with CFL numbers of order \(10^2\) in stationary intervals and small timesteps when a perturbation enters the flow field [1404.4503, 1405.3497]. The method becomes even more efficient when coupled with a hybrid explicit–implicit strategy: use explicit steps when the equidistributed timestep would imply a small CFL and implicit steps when the timestep is large [1404.4503]. This is directly analogous to adaptive pseudo-time continuation, except that the controller is driven by an adjoint-weighted functional error rather than by a residual norm. A more localized multirate interpretation appears in CFD–DEM particle integration. The method partitions particles into boxes using orthogonal recursive bisection and assigns each box a local timestep based on whether particles are close to collision partners [1802.09579]. For particle \(i\), the contact time is \[ \tau_p = \sqrt{\frac{m_i}{k}}, \] and if a potential collision partner lies within \(3r_p\), the local timestep is set to \[ \Delta t_{p,i} = \frac{\tau_{p,i}}{50}. \] Otherwise the particle uses a much larger sub-iteration step \(\Delta t_{sp}=\Delta t_f/\text{nsubit}\) [1802.09579]. Each box takes its own local timestep \(\text{ltstep}[i]\), which is the minimum over the particles in that box, and advances independently within each fluid time slab [1802.09579]. This is a clear example of local pseudo-time allocation in a multiphysics setting: the global physical time is fixed by the fluid step \(\Delta t_f\), while the DEM subsystem evolves in locally adaptive pseudo-time substeps whose density reflects local collisional stiffness. A related space–time adaptive strategy appears in convection–diffusion finite differences [2101.00659]. There, the scheme starts with the highest-order centered spatial discretization and RK\(_4\) in time, using a global time step given by the largest CFL allowed by linear stability, then detects local instabilities a posteriori and selectively switches troubled nodes to an upwind spatial scheme and a more dissipative RK\(_D\) time scheme while keeping the same global \(\Delta t\) [2101.00659]. This suggests a different kind of adaptive pseudo-time stepping: instead of shrinking the global step, one modifies the local scheme so that the large step remains admissible. ## 6. Cross-cutting principles and recurrent trade-offs Across these disparate formulations, several recurrent principles define adaptive pseudo-time stepping. The first is that step-size control is tied to a structural quantity rather than to time alone. In residual networks, the decisive combination is \(\|\mathbf W_j\|_2\Delta t_j\) [1911.10305]. In Boris–Buneman particle acceleration, it is the maximum self-field acceleration \(A_{\max}\) through \(g\propto A_{\max}^{-1/2}\) [1211.3664]. In non-Lipschitz SDEs and SPDEs, it is the size of the drift and diffusion responses through bounds such as \[ h_n = \frac{\bar h}{1+\lfloor f(X_n)\rfloor + \lfloor g^2(X_n)\rfloor} \] or \[ \Delta t_{n+1}\le \frac{\delta}{\|F(X_h^n)\|^{1/\theta}} \] [1706.03098, 1812.09036]. In rate-independent systems, it is the complementarity residual and the energy-dissipation residual [2204.05860]. In hyperbolic PDEs, it is the adjoint-weighted time contribution to a goal functional [1405.6510, 1404.4503, 1405.3497]. The second principle is that pseudo-time is often more fundamental than physical time. In Hamiltonian and accelerator problems, the evolution is explicitly rewritten in terms of \(\tau\) with \(dt/d\tau=g(z)\) [1108.0322, 1211.3664]. In rate-independent systems, \(s\) is the primary evolution variable and physical time is recovered as part of the solution [2204.05860]. In residual networks, depth itself plays the role of pseudo-time [1911.10305]. Even when the papers work in physical time, as in PFC or weakly instationary Euler flows, the adaptive mechanism behaves as if it were allocating pseudo-time density to where the solution or the target functional changes most [2008.00212, 1405.3497]. The third principle is that adaptivity can preserve or destroy structure depending on how it is introduced. Variable steps inserted naively into symplectic schemes can cause resonant instabilities or destroy canonical Hamiltonian form [1108.0322]. Semi-implicit BDF2 can exhibit a limiting step and only “nearly” converge to steady state because of conditional stability [1703.10297]. By contrast, adaptive variable-step BDF2 for the PFC model preserves a modified energy law under the sharp ratio bound \(r_k<3.561\) [2008.00212], and the rate-independent pseudo-time algorithm converges to parametrized balanced viscosity solutions as the tolerance tends to zero [2204.05860]. The fourth principle is that training-time or offline adaptation can be separated from runtime cost. The ResNet controller is used only during training; after training, the learned pseudo-time schedule is frozen and incurs negligible inference overhead [1911.10305]. In accelerator simulation, the adaptive rule is algebraic in the current state and cheaper than repeated local error estimates or rejected steps [1211.3664]. In adjoint-based Euler timestep control, the dual solve is performed on a coarse grid, and the expensive fine-grid run uses the resulting step sequence [1404.4503, 1405.3497]. This suggests a broad practical pattern: adaptive pseudo-time stepping is often most useful when the adaptation mechanism shapes the evolution law while leaving the dominant runtime kernel unchanged. A plausible implication is that “adaptive pseudo-time stepping” is best regarded not as a single algorithm but as a design paradigm. Its core ingredients are: a surrogate evolution variable; a structurally meaningful step-selection signal; and a discretization or controller that allows nonuniform stepping without destroying the qualitative behavior one wishes to preserve. The literature surveyed here shows that this paradigm can be instantiated by learned channelwise depth increments [1911.10305], state-dependent geometric time rescaling [1108.0322, 1211.3664], positivity- and stability-preserving explicit step control [1706.03098, 1812.09036, 1205.2751], energy-law-aware variable-step multistep schemes [2008.00212, 2504.00405], and adjoint-driven goal-oriented timestep allocation [1405.6510, 1404.4503, 1405.3497]. Source: https://www.emergentmind.com/topics/adaptive-pseudo-time-stepping