Large Distant Gradients Need Not Be Reliable: reliability-weighted credit assignment for long-horizon autoregressive forecasting

Introduction

In autoregressive forecasting, training over long closed-loop rollouts provides supervision at distant forecast horizons, while backpropagation through time (BPTT) carries gradients from those distant losses backward through many earlier autoregressive steps. Classical analyses show that these distant gradients can grow rapidly in recurrent models [8, 40], and related sensitivity amplification has been observed in neural PDE solvers and weather models [31, 41]. But a large distant gradient need not be a reliable one ([39]). We ask when distant BPTT gradients acquire high leverage without acquiring correspondingly reliable learning signal, and how long-horizon supervision can be retained without trusting every backward contribution equally.

To see why gradient magnitude can be misleading, first consider the history-dominated case: an autonomous system without external forcing or control inputs. The loss at each future step produces a gradient that is propagated backward through earlier applications of the same model, and the contributions from all forecast steps are added to form the parameter update. The same Jacobian products that amplify a contribution’s predictable learning signal also amplify its unpredictable innovation. Consequently, a distant contribution can become large enough to dominate the summed update even when most of its magnitude reflects unreliable variation. Gradient magnitude therefore measures backward sensitivity, not necessarily useful credit.

Under full BPTT, all delayed contributions remain open and are summed into the update (Figure 1a). Existing controls address neighboring problems but intervene at different stages of the training process. Gradient clipping ([40]) and adaptive optimizers such as Adam ([22]) rescale gradients after forecast-step contributions have already been combined and therefore cannot identify which delayed contribution is unreliable (Figure 1b). Forward Jacobian regularization instead limits amplification by constraining the learned dynamics ([41]; [36]) (Figure 1c), whereas truncation or shorter rollouts limit delayed backpropagation but also remove distant supervision ([1]) (Figure 1d). None of these approaches directly asks which delayed contributions remain reliable before they are aggregated. This motivates a backward-only intervention that preserves the complete rollout, all horizon losses, and the learned forward dynamics while selectively reducing unreliable delayed credit (Figure 1e).

Internal-DW implements this reliability-weighting idea inside residual blocks. At a block $y = x + F_{\theta}(x)$, ordinary BPTT maps an incoming gradient through $I + J_{F}^{\mathsf{T}}$; Internal-DW replaces this backward operator with $\alpha I + mJ_{F}^{\mathsf{T}}$, assigning separate reliability weights to the identity and nonlinear routes while leaving the block’s forward value unchanged. The two gains are chosen jointly by a box-constrained $2 \times2$ Wiener solve that trades distortion of the conditional gradient signal against unpredictable variation. Because signal and noise are not separately observed, the solve requires an explicit noise estimate ([32]; [6]; [37]). We therefore make the identifying assumption explicit through two noise samplers: DW-Generic uses centered training residuals as a proxy for unpredictable error, whereas DW-Prior incorporates domain-supported structure when available.

We first test the mechanism in a controlled system where gradient signal-to-noise ratio (SNR) is known. Delayed gradients grow even as their SNR falls, the estimated gains reduce held-out route-gradient risk, and Internal-DW achieves the lowest forecast error among the tested methods. On application data, where the true signal–innovation split is unavailable, held-out gradient tests ask whether large per-horizon gradients transfer across examples, while controlled-noise tests ask whether the estimated gains decrease as noise increases. Across four history-dominated, weak-drive testbeds, Internal-DW improves forecasting and extends or preserves the useful training-horizon range. Its benefit weakens when usable history is limited or when important driven information is not represented by the sampler. These cases define boundaries of the current estimator, not of the broader delayed-credit reliability problem.

We therefore view unreliable delayed credit as the broader training problem, and Internal-DW as one reliability-weighted correction for it. The current estimator is best suited when useful credit is primarily propagated through state history; separating that credit from newly injected drive remains an important extension.

Related Work

Long autoregressive rollouts provide distant supervision but can make BPTT gradients large, noisy, and expensive to compute. Existing approaches intervene at different stages of long-horizon training. Scheduled sampling, DAgger, Professor Forcing, and pushforward training change the distribution of states on which the model is trained ([7, 45, 23, 10]). Jacobian penalties, stable recurrent operators, and recurrent training instead constrain or stabilize the forward dynamics ([31, 41, 36, 57]). Studies of training horizons examine trade-offs among forecast accuracy and learnability ([1]). TBPTT truncates backward propagation, yielding a cheaper but biased gradient ([48]); stop-gradient unrolling retains long-rollout states without propagating through time ([25]), while damped BPTT prescribes decay during gradient accumulation ([19]).

Long-term credit assignment also arises in reinforcement learning and agent training. RUDDER redistributes delayed rewards through return decomposition, whereas TRACE derives turn-level rewards along long agent trajectories ([3, 50]). GAE uses a single global temporal parameter $\lambda$ to balance shorter, lower-variance and longer, lower-bias estimates of policy advantage ([47]). Differentiable world models optimize policies by backpropagating through imagined trajectories, and Ma et al. redesign the model to provide more direct long-range gradient paths ([17, 29]). Internal-DW instead keeps the supervised autoregressive rollout, horizon losses, and learned forward dynamics fixed, and estimates the reliability of transported internal VJPs before their route contributions are aggregated. Appendix A discusses additional connections to finite-time sensitivity, recurrent-gradient estimators, covariance identification, and backward-message methods.

Reliable Gradient Transport in Long-Horizon Rollout

Long-horizon BPTT can give delayed gradient contributions high leverage without making them reliable. We formalize this mismatch, derive Internal-DW as a backward-only reliability weighting of residual routes that leaves the forward rollout unchanged, and show why estimating its gains requires an explicit noise model whose assumptions determine the scope of the current estimator.

Large Delayed Gradients Need Not Be Reliable

Jacobian transport controls a delayed gradient’s magnitude, not its reliability, and one amplified contribution can dominate the summed update. For the $K$-step loss $\mathcal{L}=K^{-1}\sum_{k=1}^{K}\lVert r_{t+k}\rVert^{2}$, where $r_{t+k}=\hat{x}_{t+k}-x^{\star}_{t+k}$, let $\mathrm{D}_{t+k}=\mathrm{d}\hat{x}_{t+k}/\mathrm{d}\theta$. Define $\mathbf{J}_{t}:=\partial\hat{x}_{t}/\partial\theta$, and $\Phi_{a,q}=\mathbf{J}_{a+q-1}\cdots\mathbf{J}_{a}$, with $\Phi_{a,0}=\mathbf{I}$. The exact gradient and its chain-rule expansion are

Because the same parameters are reused at every rollout step, each step injects a new local contribution while transporting earlier ones through Jacobian products. Appendix B formalizes this forced sensitivity recurrence and shows that, under a stationary local linearization, sensitivity can plateau, accumulate, or grow exponentially. (1) sums contributions from all such paths, so a dominant transported term can control the update’s magnitude and direction ([33]).

Write $r_{t+k}=b_{t+k}+\eta_{t+k}$, where $b_{t+k}=\mathbb{E}[r_{t+k}\mid\hat{x}_{t+k-1}]$ is the learnable conditional-mean component and $\eta_{t+k}$ is unpredictable innovation. The same Jacobian product transports both components, so amplification increases the leverage of innovation without making it informative. Even without exponential growth, innovation comparable to the predictable component can perturb the update. We therefore seek a backward-only linear route weighting that minimizes mean-squared error relative to the predictable delayed contribution, before rollout-step contributions are aggregated.

Internal-DW: Optimal Reliability-Weighted Gradient Transport

Internal-DW changes only backward propagation at residual merges; the forward rollout is unchanged. We first define its residual-route operator and then show that, under an explicit signal–innovation decomposition, its oracle weights are the unique minimum-risk bounded linear weighting of the two backward routes. Figure 2 summarizes the resulting method.

Backward-only routing. For a residual block $y = x + F_{\theta}(x)$, ordinary BPTT maps an incoming gradient $v$ to $v + J_{F}^{\mathsf{T}}v$. Internal-DW instead uses

We assign one pair $w_{k,\ell} = (\alpha_{k,\ell},m_{k,\ell})^{\mathsf{T}}$ to each horizon $k$ and residual layer $\ell$.

Wiener gains. At a fixed $(k,\ell)$, let $v_{k,\ell}$ be the gradient arriving at the output of the residual block. The identity and nonlinear backward routes return

Collect these two route gradients as $d = [d_I\ d_F]$ and decompose them as

where $\mathcal{C}$ denotes the information under which delayed credit is predictable. Thus $s$ is the predictable transported gradient and $n$ its unpredictable innovation. Then $\mathbb{E}[s^{\mathsf{T}}n] = 0$, and

Theorem 1 (Optimal reliability-weighted gradient transport). For $w = (\alpha,m)^{\mathsf{T}}$, let $dw$ be the routed message and $s\mathbf{1} = s_I + s_F$ the fully open predictable target. If $T = P + R \succ0$, then the route risk decomposes as

The corresponding unconstrained and bounded optima are, respectively,

Hence $w^\star$ is the $T$-metric projection of $w_{\mathrm{unc}}$ onto $[0,1]^2$. Since the fully open BPTT operator $w=\mathbf{1}$ is feasible,

with strict inequality whenever $\mathbf{1}$ is not itself optimal.

The first term in (4) measures distortion of predictable credit from closing the routes, while the second measures innovation retained by leaving them open. Thus Internal-DW is the minimum-risk bounded linear trade-off between preserving predictable credit and suppressing transported innovation, rather than a prescribed gradient decay. For a single route, the unconstrained solution reduces to $P/(P+R)$ [55]. The box constraint prevents sign reversal and route-wise amplification. Appendix C.1 gives the proof and Appendix C.4 gives the exact constrained solve.

Plug-in estimation. The oracle moments $P$ and $R$ in Theorem 1 are not separately observed. Fully open total-moment and noise probes produce route messages $d^{\mathrm{tot}}$ and $d^{\mathrm{noise}}$, with the latter obtained by transporting a sample from $\mathcal{Q}_d$. Suppressing $(k,\ell)$, we estimate

The total-moment probe uses the quadratic output loss specified in Appendix C.3, independently of the task loss used for parameter updates. Substituting $(\widehat{P},\widehat{R})$ into (4) and solving over $w \in[0,1]^2$ gives the applied Internal-DW gains. Section 3.3 states the assumptions on $\mathcal{Q}_d$; Appendix Sections C.2–C.5 give the probes, implementation, and constrained solve; and Appendix Sections D.1–D.4 give the sampler constructions.

Noise Identifiability and Sampler Assumptions

Theorem 1 determines the oracle Internal-DW route weights $w^\star$ from the latent decomposition $T=P+R$. However, the open route law does not identify this decomposition by itself. Thus identifying structure is required; Internal-DW supplies it through the noise model $\mathcal{Q}_d$. To see the ambiguity, for any $0 \preceq R \preceq T$, independent variables $s \sim\mathcal{N}(0,T-R)$ and $n \sim\mathcal{N}(0,R)$ produce the same observed distribution $z=s+n \sim\mathcal{N}(0,T)$, while different choices of $R$ can yield different unconstrained optima $w_{\mathrm{unc}}=T(T-R)\mathbf{1}$. Hence even the full distribution of the observed route vector does not, in general, identify its predictable–innovation split, and therefore does not uniquely specify the oracle routing problem without additional assumptions on the dynamics, observations, or covariance structure. Appendix E.1 gives the formal construction. We encode the identifying assumption in $\mathcal{Q}_d$, which generates the output-noise trajectories used to estimate $R$.

DW-Generic. DW-Generic uses centered, lagged training-residual trajectories, with a random trajectory-wide sign, as samples of output noise. Its assumption is correspondingly strong: centered prediction residuals must be a useful proxy for unpredictable output variation. If they contain systematic error that the forecasting model could still learn, DW-Generic incorrectly treats that variation as noise.

DW-Prior. DW-Prior uses either cross-fitted residuals from a domain-supported predictor or samples from a domain-supported covariance model. The predictor-based version treats held-out variation unexplained by the predictor as noise and is appropriate when the predictor captures the conditional mean under the declared conditioning information. The covariance-based version treats zero-mean draws from the fitted covariance model as noise and is appropriate when the selected covariance family captures the relevant dependence. For the spatial samplers used here, this requires approximate translation stationarity and a spectrum that remains stable during calibration.

Exact sampler constructions and implementation details are given in Appendix Sections D.1–D.4. Appendix Table 1 summarizes the dataset-specific assumptions, supporting references, and retained covariance. For each dataset, we select between DW-Generic and DW-Prior using validation loss and then fix that choice for test evaluation.

Data	Assumed structure	Sampler $Q_d$ and retained dependence
DW-Generic	Centered residuals	Lagged residual bootstrap (Equation 17); coordinates and horizons
Mackey–Glass, NARMA, ETTm1, ETTm2	Centered residuals	Lagged residual bootstrap (Equation 17); coordinates and horizons
DW-Prior	Long-range temporal dependence [24]; cross-electrode factor model (ours)	Cross-fitted factor model; electrodes and horizons
iEEG theta	Long-range temporal dependence [24]; cross-electrode factor model (ours)	Cross-fitted factor model; electrodes and horizons
HCP movie fMRI	Shared stimulus response [11], [35]; AR residual model (ours)	Subject cross-fit + AR model; regions and horizons
Shear flow, WeatherBench-2	Approximate translation stationarity, motivated by stationary-flow spectral analysis [2], [51]	Lagged 2-D spectrum; spatial modes, horizons independent

Let $R_0$ denote the oracle noise moment, $R_{\mathcal{Q}}$ the population noise moment targeted by the selected sampler, and $\widehat{R}$ its finite-sample estimate. Their difference decomposes as

Thus data can reduce the second term in (8), but cannot remove the first when the sampler model is misspecified. In particular, DW-Generic can count predictable model error as noise, while DW-Prior can omit predictable structure or covariance outside its chosen representation. The next subsection explains how these two error sources determine the scope of the current estimator.

Scope and Boundaries of Internal-DW

A first scope boundary arises when state history supplies little reliable delayed signal and no other signal source compensates. In this case, the conditional signal moment $P$ is small relative to the oracle noise moment $R_0$, so the Wiener solution assigns little weight to the corresponding routes. Internal-DW then has little reliable delayed credit to preserve.

A boundary arises from sampler-model mismatch. Useful drive-dependent information should contribute to $P$ and favor a more open oracle route. If $Q_d$ does not account for propagated history, observed drive, and their interaction, however, predictable driven variation can instead inflate $R_Q$ relative to $R_0$ along route directions. The applied gain may then close too far and suppress useful gradients. Even under a correctly specified sampler, finite data contribute the estimation error $\widehat{R}-R_Q$. Appendix E.2 gives consistency conditions for the residual-based estimator, while Appendix E.3 shows how route-moment errors perturb the Wiener gain and contribute to excess local route risk. The first case is a low-signal scope boundary; sampler-model mismatch and finite-sample error are boundaries of the current estimator, not of the broader delayed-credit reliability problem.

Experimental Setup

Testbeds. We use one controlled system with known SNR to examine the mechanism and eight application testbeds to evaluate forecasting performance. The controlled system is a stationary, undriven linear-Gaussian AR process. The application testbeds are Mackey–Glass (MG) for autonomous delayed dynamics [30]; NARMA for driven nonlinear control [4]; iEEG and movie fMRI for neural recordings [21]; [53]; ETTm1 and ETTm2 for temporal forecasting [58]; and The Well shear flow and WeatherBench-2 for gridded fluid and weather fields (Ohana et al., 2024; [43]). Appendix F.1 details the controlled system; Appendix Sections G.1–G.2 cover the application datasets and training details, respectively.

Predictive regimes. Figure 3 maps predictive gains from simple held-out probes using observed future drive alone versus drive plus aligned history. Appendix Sections H.1–H.3 detail their construction, interpretation, and validation under controlled drive. The probes identify three predictive regimes and one controlled intervention:

• History-dominated, weak drive: MG, shear flow, ETTm1, and ETTm2.

• Drive-dominated, weak history: NARMA, movie fMRI, and WeatherBench-2.

• Weak drive and weak history: iEEG.

• Joint history-and-drive intervention: driven MG (see Appendix I for the driven-system configuration and intervention results).

These groupings describe predictive structure recovered by the probes, not intrinsic dataset properties, and predict the clearest Internal-DW benefit in the history-dominated, weak-drive regime.

Evaluation and model selection. Estimator variants and checkpoints are selected by validation loss and fixed before test evaluation. For a model trained with horizon $K$, each held-out origin is rolled out through $H_{\mathrm{eval}} = \lceil1.5K \rceil$. For test units $\mathcal{U}$ and origins $\mathcal{O}_u$, the primary metric is

Appendix J gives the origin-selection details.

Baselines and controls. We compare Internal-DW with full BPTT, global gradient clipping (Clip), and forward-Jacobian regularization (JReg). On the four history-dominated, weak-drive testbeds, we additionally compare with TBPTT and a tied static gain $\alpha=m=c$ (Static), analogous to fixed backward damping [19]. The TBPTT segment length, Static coefficient $c$, and choice between DW-Generic and DW-Prior are selected by validation loss and then fixed for test evaluation. Within each dataset, all methods otherwise use the same experimental settings. Appendix Table 1 lists the selected sampler, and Appendix G gives architectures, selection grids, and training details.

Results

The experiments first test whether delayed gradients become less reliable and whether Internal-DW responds appropriately, then determine when this local correction improves forecasting across predictive regimes and training horizons.

Mechanism and Estimator Response

Known-SNR mechanism closure. In the known-SNR system, independently resampling future innovations from a fixed $x_t$ gives the exact conditional-mean gradient and therefore direct measurements of gradient SNR and prefix risk. A decrease in prefix risk means that an added forecast step contributes more conditional-mean signal than innovation; an increase means the reverse. Appendix Sections F.1–F.2 give the construction. Figure 4 follows the resulting chain from gradient amplification to forecasting performance.

Figure 4a–b shows that later gradients grow while their SNR falls. As the included horizon increases, omitted signal decreases but accumulated innovation grows; beyond $K^\star_{\mathrm{grad}}=24$, the marginal innovation cost exceeds the recovered signal and total prefix risk rises. Reliability also fluctuates across intermediate horizons, so neither forecast distance nor a monotone decay schedule identifies which contributions remain useful.

The estimated gains have mean relative error 0.298, whereas fixing all other moments to their true values and estimating only $\hat{R}$ reduces the error to 0.0013. Thus, most gain error comes from finite-sample error in $\hat{T}$, rather than sampler mismatch or finite-sample error in $\hat{R}$. Nevertheless, Internal-DW attains the lowest held-out route risk; shifting the gains across forecast steps increases risk, showing that their placement matters beyond generic attenuation. The transfer-order control also reduces risk, although less than Internal-DW, showing that gains computed from true calibration-route moments transfer imperfectly to test routes. Internal-DW also achieves the lowest forecasting error among full BPTT, clipping, Jacobian regularization, and TBPTT.

Application-data diagnostics. On application data, the true gradient SNR and oracle gains are unavailable, so we use two complementary held-out diagnostics. For each forecast step $k$, $A(k)$ measures the magnitude of that step’s gradient relative to the first forecast step, while $U(k)$ measures whether its direction is useful on a disjoint held-out rollout. Positive $U(k)$ indicates that a small update along the negative step-$k$ training gradient would locally reduce held-out rollout loss, whereas $U(k) < 0$ indicates that it would increase held-out loss. Appendix K gives the exact definitions and evaluation protocol. Figure 5 therefore asks whether gradients that become large at distant forecast steps also remain useful on unseen examples, while panel (b) asks whether the estimated gains close as independent noise increases.

Figure 5 shows that gradient amplification and held-out utility can decouple on application data. Across most testbeds, per-horizon gradients grow even as their held-out utility fluctuates or becomes negative; ETTm1 and WeatherBench-2 provide complementary cases in which growing gradients remain mostly useful. This contrast shows that neither gradient magnitude nor forecast distance alone identifies a transferable delayed contribution. fig:5b supports the corresponding estimator-response claim: both route gains decrease on every testbed as controlled noise increases. The dataset-specific operating points remain informative. WeatherBench-2 stays nearly open, whereas NARMA is nearly closed even under weak perturbation, consistent with its selected sampler treating predictable drive-dependent residual variation as noise. Together, Figure 5 provides application-level support for two linked claims: large delayed gradients need not remain useful on held-out data, and Internal-DW closes its routes in response to increased noise. Because native-data gradient SNR and oracle gains remain unobserved, these diagnostics support the qualitative mechanism and estimator response rather than establishing exact gain optimality.

Forecasting Effects Across Regimes and Horizons

Forecasting across predictive regimes. Forecasting outcomes follow the applicability diagnostic in Figure 3. Figure 6 shows that Internal-DW improves all four history-dominated, weak-drive testbeds over full BPTT, reducing relative $L_2$ by 5.16%–13.76%. It outperforms validation-selected TBPTT on Mackey–Glass, ETTm1, and ETTm2, and the Static control on ETTm1, ETTm2, and shear flow. Thus neither truncation nor a single global damping coefficient consistently reproduces the gains.

At the identified boundaries, Internal-DW provides no consistent benefit. NARMA’s near-closed gains coincide with 7.87% worse forecasting, consistent with the sampler assigning predictable drive-dependent variation to noise. iEEG and WeatherBench-2 remain comparable to full BPTT, while movie fMRI worsens. A controlled Mackey–Glass intervention strengthens this interpretation: adding strong observed drive while retaining substantial history value moves the system across the regime map and reverses Internal-DW from a 5.17% improvement to a 7.71% degradation relative to full BPTT (Appendix I). Detailed runtime measurements are reported in Appendix M.

Fitted optimal training-horizon range. Longer training horizons add both conditional-mean signal and innovation, creating a marginal trade-off that reliability weighting can reshape but not eliminate. We vary the training horizon $K$ while evaluating all checkpoints at a fixed forecast horizon $H$ across MG, ETTm1, ETTm2, and Shear. Figure 7 and the seed-specific curves in Appendix N show that $K^\star_{\mathrm{DW}} \ge K^\star_{\mathrm{Full}}$ for every matched seed and in the seed-averaged fits, with a lower fitted minimum error for Internal-DW on all four datasets. This is consistent with the reliability trade-off in Theorem 1: suppressing transported innovation while preserving predictable delayed credit can increase the marginal value of extending the rollout. Internal-DW therefore extends or preserves the useful training-horizon range rather than improving every horizon pointwise.

Discussion and Conclusion

Internal-DW minimizes a local quadratic route risk rather than the observed rollout loss, so it does not guarantee a better training trajectory. Its gains are only as accurate as the selected noise sampler, and solving each horizon and layer independently ignores interactions after their contributions are combined. These limitations suggest two extensions. Samplers that model innovation conditional on history, observed drive, and their interaction—or architectures that separate these pathways—could broaden the favorable regime. Richer conditional samplers and justified pooling across related channels, horizons, or layers could also reduce model mismatch and finite-sample error.

Long-rollout BPTT can amplify distant gradient contributions without making them reliable. Internal-DW addresses this mismatch by reliability-weighting residual backward routes while preserving the forward rollout and all horizon losses. The noise sampler provides the identifying structure needed to estimate latent route reliability, allowing domain-appropriate samplers to be plugged into a common routing framework. The known-SNR experiment connects gradient amplification, route-risk reduction, and forecasting improvement, while the application results show benefits primarily in history-dominated, weak-drive regimes. Internal-DW can also extend or preserve the useful training-horizon range, although it need not improve every $K$ pointwise and can be neutral or harmful when usable history is weak or the sampler is mismatched to dominant driven information. Thus Internal-DW is not a universal stabilizer; it shows that long-horizon supervision can be retained without trusting every backward contribution equally.

AI USE STATEMENT

We used generative AI tools to help refine hypotheses and diagnostics, assist with figure preparation and improve the organization and wording of the manuscript. The authors checked the resulting derivations, references, figures, and claims against the cited sources and experiments. The authors take responsibility for the final content of the paper and all AI-assisted artifacts.

ETHICS STATEMENT

This work performs secondary analysis of existing human neurophysiological data and includes no new participant recruitment or data collection. The OpenNeuro iEEG release contains publicly shared de-identified recordings for which participants or their legal guardians consented to research use and public sharing, under the approval reported by the originating study ([21]). The HCP movie-fMRI data were collected with informed consent and institutional approval ([53]); we used the open-access release under its data-use terms, did not use restricted identifying variables, and made no attempt to re-identify participants. The experiments study forecasting methodology and are not intended for clinical decision making or subject-level interpretation.

REPRODUCIBILITY STATEMENT

Section 3.2 specifies the backward operator, calibration target, and estimator assumptions. The appendices provide derivations, noise-sampler constructions, implementation details, dataset preprocessing and splits, training hyperparameters, evaluation protocols, and seed-level results. Checkpoints retain the route gains, estimated moments, probe counts, and sampler configuration needed to audit calibration.

Extended Related Work

Sensitivity in autonomous and driven rollouts. Dynamical-systems analyses distinguish asymptotic instability from realized finite-time amplification, which can depend strongly on the starting state, forcing, observation error, and model error [26, 38]. The same Jacobian products govern BPTT sensitivity, but the magnitude of a transported gradient also depends on the residual and its alignment with amplified directions [40, 33]. Related work further shows that forecast errors can concentrate in unstable subspaces and that driven recurrent systems can exhibit generalized synchronization with their inputs [52, 27]. These results characterize the leverage of delayed information, not the reliability of the resulting gradient.

Credit assignment in reinforcement learning and differentiable world models. Long-term credit assignment also arises in reinforcement learning, where much of the literature asks how delayed returns should be attributed to earlier states or actions. PIPPS uses total propagation to combine gradient estimators across pathwise derivative depths, assigning greater weight to lower-variance estimators [39]. RUDDER redistributes delayed rewards through return decomposition, while synthetic-return methods associate past states directly with distant future rewards [3, 42]. Recent LLM-agent training methods similarly decompose long interaction trajectories into training transitions or construct turn-level rewards [28, 50]. Differentiable world models provide a closer computational connection. Stochastic value-gradient methods differentiate policy objectives through learned dynamics, while Dreamer backpropagates value gradients through latent imagined trajectories [18, 17]. Ma et al. further redesign the world-model architecture to provide more direct long-range policy-gradient paths [29]. Collectively, these works alter how delayed returns are assigned or how policy gradients are propagated through model rollouts. Internal-DW instead retains the observed per-horizon targets, forward rollout, and all horizon losses, and estimates the reliability of transported identity and nonlinear route VJPs.

Alternative recurrent-gradient estimators. TBPTT shortens the backward window, whereas ARTBP randomizes the cutoff and reweights surviving paths to obtain an unbiased estimator of full BPTT at the cost of additional variance [48]. RTRL and its scalable variants propagate recurrent gradients online [56, 49, 34]; while activation checkpointing trades recomputation for memory without changing the full BPTT gradient [14, 12, 15]. These methods change how a long-horizon gradient is computed or approximated; they do not estimate the reliability of individual route contributions after transport.

Related methods that weight or modify backward signals. Kalman-inspired neural models learn filtering updates for latent-state estimation: Recurrent Kalman Networks maintain structured latent uncertainty, and KalmanNet learns a Kalman gain under partially known dynamics [5, 44]. Classical adaptive filtering also shows that process and observation covariances cannot generally be recovered from one total covariance without additional dynamical or rank assumptions [32, 6, 37]. This identification principle motivates the explicit sampler contract in Appendix E.1. Synthetic gradients and related local methods instead learn or approximate module-level backward messages [20, 54], while Highway-BP uses residual shortcuts to approximate ordinary backpropagation more efficiently [13]. Internal-DW retains the rollout VJPs and applies covariance-derived reliability weights to routes within each residual merge.

Why Delayed Gradients Accumulate

The sensitivity in (1) obeys the exact forced recurrence

Thus contraction of an old path does not remove the new parameter contribution injected at the next rollout step. To isolate the resulting long-horizon regimes, consider a stationary local linearization and one parameter direction,

Here $\mu$ is the coherent part of the local parameter injection and $\epsilon_k$ is its temporally uncorrelated component. Iteration gives

The second expression is a finite controllability Gramian. If $\rho(J) < 1$, both sums converge and the root-mean-square (RMS) sensitivity approaches a finite plateau. A unit-eigenvalue mode excited by $\mu$ accumulates in the mean, while any semisimple unit-modulus mode excited by $Q$ accumulates in the variance. If $\rho(J) > 1$ and an unstable mode is excited, the corresponding moment grows exponentially. A nontrivial Jordan block on the unit circle can produce still faster polynomial growth. These statements are matrix results; no scalar reduction is required.

The scalar case makes the limits explicit. For a nonnegative gain $a$, coherent injection $\mu$, and innovation variance $\sigma^2$,

with the continuous limits $k$ at $a = 1$. Hence $|a| < 1$ gives saturation, $a = 1$ gives linear accumulation in the coherent component and $\sqrt{k}$ RMS growth in the uncorrelated component, and $|a| > 1$ gives exponential growth. Equal deviations around unity are also asymmetric at long horizons: for $0 < \delta< 1$ and $k > 1$,

so a comparably expansive mode eventually dominates a comparably contractive one. This comparison does not assert that learning makes the two deviations equally likely.

These regimes describe the shared-parameter sensitivity. They carry over to the gradient in (1) when the forecast residual excites the corresponding transported modes. Time-varying Jacobians, correlated injections, and cancellation can change the rate, so these are sufficient regimes rather than a claim that every nonlinear autoregressive model produces large gradients.

Internal-DW: Theory and Implementation

This section first proves the optimality result underlying Internal-DW and then gives the complete training and calibration procedure. Algorithm 1 summarizes how ordinary training steps alternate with fully open calibration probes that update the route moments and gains. The remaining subsections specify the routing placement, probe construction, constrained solve, and calibration cost.

Proof of Theorem 1

Proof. Recall that $d = s + n$, where $s = \mathbb{E}[d \mid C]$ and $n = d - s$. Hence $\mathbb{E}[n \mid C] = 0$, and therefore

For any $w \in\mathbb{R}^2$,

Using this orthogonality,

which gives (4).

Expanding and using $T = P + R$ gives

Since $T \succ0$, $\mathcal{R}$ is strictly convex, and its unique unconstrained minimizer satisfies

Thus

Using $T w_{\mathrm{unc}} = P\mathbf{1}$, completing the square yields

Therefore

so $w^\star$ is the $T$-metric projection of $w_{\mathrm{unc}}$ onto $[0,1]^2$.

Finally, $\mathbf{1} \in[0,1]^2$, so optimality of $w^\star$ gives

From the risk decomposition,

and hence

Because $\mathcal{R}$ is strictly convex, the constrained minimizer is unique; consequently, the inequality is strict whenever $\mathbf{1}$ is not itself optimal.

For a single route, $T = P + R$ is scalar, so the unconstrained solution reduces immediately to $w_{\mathrm{unc}} = P/(P + R)$.

Algorithm 1 Internal-DW training and calibration

Require: Forecaster $f_\theta$, sampler $\mathcal{Q}_d$, and stored $\{w_{k,\ell}, \widehat{T}_{k,\ell}, \widehat{R}_{k,\ell}\}$

1. Use $w_{k,\ell} = (1,1)^{\mathsf{T}}$ during the fully open warm-up

2. for training minibatch $b$ do

3. Roll out $K$ steps and form the task loss

4. if $b$ is a calibration minibatch after warm-up then

5. Form $q_{\mathrm{tot}}$ and $q_{\mathrm{noise}}$ and run two fully open VJP probes

6. Stage updates of $\widehat{T}_b$, $\widehat{R}_b$, and $\widehat{P} \leftarrow\Pi_{\mathrm{PSD}}(\widehat{T}_b - \widehat{R}_b)$

7. Stage new gains by minimizing (4) over $w \in[0,1]^2$

8. end if

9. Backpropagate the task loss with the previously committed gains

10. Update $\theta$ and commit any staged sampler statistics, moments, and gains

11. end for

Where the Route Weights Are Applied

At each forecast step and residual layer, two forward-identity autograd nodes separately weight the identity and nonlinear backward messages, record their unscaled VJPs, and leave the forward value unchanged. The recurrent-state path uses the nonlinear coefficient. Routing is confined to internal residual merges; the outer autoregressive update remains fully open because a single gain there would combine paths represented separately by the two-route estimator.

How the Route Moments Are Calibrated

For calibration sample $i$ and forecast horizon $k$, the total-moment and noise probes are

The quadratic $q_{\mathrm{tot}}$ is the common calibration probe used to estimate route moments and need not equal the task loss used for parameter updates. Both probes are differentiated with all route gains set to one. For the residual block

the recorded two-route message for a scalar probe $q$ is

Accordingly,

These messages update the EMAs in (7). If the block input has dimension $D_{\ell}$, the implementation accumulates

The common factor $1/D_{\ell}$ scales $\widehat{T}$, $\widehat{R}$, and $\widehat{P}$ equally and therefore does not change the Wiener solution.

The probes estimate local route moments under the selected sampler; the resulting Wiener objective is a local MMSE criterion, not a claim of global optimality for the nonlinear training objective.

How the Two Route Weights Are Solved

Each symmetric $2 \times2$ route matrix is stored as three scalars $(II, IJ, JJ)$. Before solving, $\widehat{T}$ and $\widehat{R}$ are symmetrized and projected onto the positive-semidefinite cone, after which

A numerical ridge of $10^{-6}\operatorname{tr}(\widehat{P}+\widehat{R})$ is added only to stabilize the two-dimensional inverse.

The constrained solution is not obtained by clipping the unconstrained one, because correlations between routes can make coordinatewise clipping suboptimal. Instead, the implementation evaluates the feasible unconstrained optimum, the optimum on each edge, and all four corners of $[0,1]^2$, and selects the point with minimum quadratic risk.

When Calibration Runs and What It Costs

Both calibration probes temporarily open every route and use autograd.grad without accumulating parameter gradients. The task backward uses the previously committed gains; newly estimated moments and gains are committed only afterward. This prevents probe noise or an in-place gain change from entering the current task update.

A calibration minibatch adds two VJP traversals. Soft gains do not prune the computational graph, so the current implementation adds periodic computation but does not reduce activation memory. Checkpoints store the gains, moments, probe counts, and sampler configuration.

How the Noise Samplers Are Constructed

Throughout this section, $x_{i,t} \in\mathbb{R}^{D}$ denotes observed sequence $i$ at time $t$, $r$ is a rollout origin, and $k \in\{1,\ldots,K\}$ is a forecast horizon. Sampler $d$ supplies a prediction $\hat{x}^{(d)}_{i,r,k}$ of the common target $x_{i,r+k}$, forms $e^{(d)}_{i,r,k}=x_{i,r+k}-\hat{x}^{(d)}_{i,r,k}$, and returns a zero-mean trajectory $\tilde{e}^{(d)}_{i,r,1:K}\sim Q_d$. Thus $d$ only labels the sampler, while $Q_d$ is the distribution of its final noise trajectories. The same draw enters $q_{\mathrm{noise}}$ in (15); thus every construction below has the common pipeline $x\to\hat{x}^{(d)}\to e^{(d)}\to\tilde{e}^{(d)}$.

Algorithm 2 Constructing a noise-sampler module

Require: Training data, module type, and any predeclared domain structure

1. if DW-Generic then

2. Center lagged $K$-step training residual trajectories

3. Return a whole-trajectory residual bootstrap $Q_d$

4. else if template-based DW-Prior then

5. Split subjects or separated time blocks into folds

6. for fold $f$ do

7. Fit the sampler predictor without fold $f$

8. Store $e^{(d)}_{i,r,k}=x_{i,r+k}-\hat{x}^{(d,-f)}_{i,r,k}$ for held-out origins

9. end for

10. Return a whole-trajectory held-out-residual bootstrap $Q_d$

11. else if spatial DW-Prior then

12. From the previous minibatch, update per-horizon/channel residual means, variances, and 2-D Fourier powers

13. Return the corresponding random-phase Gaussian sampler $Q_d$

14. end if

15. return $Q_d$, which generates $\tilde{e}^{(d)}_{i,r,1:K}$, its fitted state, and the declared assumption

Generic Sampler: Resampling Forecast Residuals

The generic sampler uses the forecaster’s own past errors as examples of plausible output noise. It resamples a whole centered residual trajectory so that dependence across output coordinates and forecast horizons is retained. Here $\hat{x}^{(G)}_{i,r,1:K}$ is the main forecaster’s detached training rollout. Let $\mathcal{B}_G$ be the bank of stored sequence–origin pairs. DW-Generic uses

$(J,S)$ is sampled from $\mathcal{B}_G$, and the independent sign $\xi_{i,r}$ is $+1$ or $-1$ with equal probability. The same sampled pair and sign are shared across horizons, preserving empirical coordinate and cross-horizon covariance. The identifying approximation is the same as in (17): any predictable error remaining in the forecaster residual is counted as noise.

Prior Samplers: Predicting Reproducible Structure First

The two neural samplers first predict structure that should repeat across separated time blocks or subjects, then treat the held-out prediction residual as noise. Cross-fitting ensures that each residual comes from a predictor that was not trained on that held-out block or subject. For a sequence–origin pair $(i,r)$ in fold $f$, sampler $d$ predicts the same targets using a model fit without that fold. Let $\mathcal{B}_d$ be the resulting bank of cross-fitted residual trajectories. Then

$(J,S)$ is sampled from $\mathcal{B}_d$, and the independent sign $\xi_{i,r}$ is $+1$ or $-1$ with equal probability. The same sampled trajectory and sign are used at every horizon. This preserves the empirical coordinate and cross-horizon covariance while making $\widetilde{e}^{(d)}$ zero mean. Each reported bank contains at most 1024 training-only cross-fitted trajectories sampled without replacement.

iEEG: Predicting Shared Temporal Factors

The iEEG sampler assumes that reproducible activity is concentrated in a small set of cross-electrode factors with predictable temporal evolution. It forecasts those factors from past activity and treats the remaining held-out variation as noise. Write each observed theta-envelope sequence as $x_{i,t}\in\mathbb{R}^{D}$. The training sequence–time pairs are divided into contiguous held-out blocks $\mathcal{H}_f$. For fold $f$, the fit set $\mathcal{I}_f$ excludes $\mathcal{H}_f$ and a purge gap on both sides, so no lag used by the predictor crosses the held-out block. From $\mathcal{I}_f$, define

Let $L_{f,\rho}\in\mathbb{R}^{D\times\rho}$ contain the top $\rho$ orthonormal eigenvectors of $C_f$, and let $a_{i,t}^{(f)}=L_{f,\rho}^{\top}(x_{i,t}-\mu_f)$ be the retained factor scores. Each factor is fit independently with a ridge-regularized AR$(p)$ model:

where $\mathcal{I}_f^{(p)}$ contains only targets whose complete lag window is inside one fit segment, $s_{f,j}$ is the mean diagonal entry of that factor’s lag Gram matrix, and $\lambda=10^{-3}$. If the fitted scalar AR has root radius above 0.995, its coefficient vector is multiplied by the largest $\kappa\in[0,1]$ whose companion matrix has spectral radius at most 0.995.

For a held-out rollout origin $(i,r)\in\mathcal{H}_f$, the observed factor history initializes the recursion

with $\hat{a}_{i,r+k,j}^{(f)}=a_{i,r+k,j}^{(f)}$ for $k\leq0$. The component outside the retained factor subspace is not set to zero. It is held at its last observed value

The residual trajectory placed in the common bank is therefore $e_{i,r,k}^{(\mathrm{iEEG})}=x_{i,r+k}-\hat{x}_{i,r,k}^{(\mathrm{iEEG})}$, which (18) converts into $\widetilde{e}_{i,r,1:K}^{(\mathrm{iEEG})}$. The pair $(\rho,p)$ minimizes the held-out MSE averaged equally over horizons and coordinates, with valid origins pooled across the purged folds:

The reported artifact uses five contiguous folds, a 128-step purge gap, $\mathcal{R}=\{8,16,32,64,80\}$, and $\mathcal{P}=\{8,16,32,64\}$; training-only cross-fitting selects $(\rho^\star,p^\star)=(80,64)$. Long-range temporal dependence in neural oscillation envelopes motivates the memory component [24]; the cross-electrode factor model is our identifying assumption.

Movie fMRI: Predicting Shared and Subject-Specific Responses

The movie-fMRI sampler uses different subjects as repeated responses to the same stimulus. It predicts the response shared across subjects together with the temporally predictable part of each subject’s deviation, and treats the remaining held-out variation as noise. Let $x_{i,t} \in\mathbb{R}^{D}$ be the parcellated response of subject $i$ at time $t$, aligned because all subjects view the same movie. We first remove each subject’s static mean and use the centered $x_{i,t}$ throughout the remaining construction:

Subjects, rather than time points, are assigned to cross-fitting folds. For a held-out subject fold $f$, the shared response and subject deviations fitted from the remaining subjects $S_{-f}$ are

For every parcel $c$, a separate AR($p$) model is fit jointly across the fit subjects:

with the same relative ridge $\lambda= 10^{-3}$, Gram-scale $s_{f,c}$, and root-radius cap 0.995 as above. Order $p = 0$ sets the predicted deviation to zero and therefore reduces to the shared response alone. For held-out subject $i \in S_f$ and rollout origin $r$, the observed deviation history initializes the recursive AR prediction $\widehat{\Delta x}_{i,r+k}^{(-f)}$, and

We select

using five subject folds and training data only; the reported artifact selects $p^{\star} = 8$. Equation 18 then returns $\widetilde{\epsilon}_{i,r,1:K}^{(\mathrm{fMRI})}$. Prior work supports the use of cross-subject alignment to capture stimulus-locked responses [11, 35]; the diagonal AR model for predictable subject deviation is our additional identifying assumption.

Spatial Fields: Matching Residual Spatial Structure

The spatial sampler treats a residual map as structured noise rather than as independent pixels. It estimates the residual scale and spatial frequency pattern separately for each horizon and channel, then uses them to shape a white Gaussian field. The Shear and WeatherBench-2 modules construct a different sampler online. Let $B_n$ contain the sequence–origin pairs in calibration minibatch $n$. For $(i,r) \in B_n$, horizon $k$, channel $c$, and grid location $u \in\Omega$, define the detached residual

Its coordinatewise mean and second moment are updated after the minibatch by

with the first observed minibatch used for initialization. Current residuals are standardized using the previous committed moments,

where $\delta_{n-1,k}=10^{-8}\max\{10^{-20},|\Omega|^{-1}C^{-1}\sum_{c,u}M_{n-1,k,c}(u)\}$ is only a numerical floor. Let $\mathcal{F}$ be the orthonormal two-dimensional real FFT and let $w(\omega)=1$ for the DC and, when present, Nyquist columns and $w(\omega)=2$ for the other stored horizontal frequencies. The minibatch power and its unit-spatial-variance normalization are

The stored color spectrum is then

again initialized by the first available normalized spectrum. At the next calibration probe, draw independent white real fields $W_{i,r,k,c}(u)\sim\mathcal{N}(0,1)$ and return

Collecting these fields over $k=1,\ldots,K$ gives the sampler output $\widetilde{\epsilon}^{(\mathrm{sp})}_{i,r,1:K}\sim Q_{\mathrm{sp}}$. The draw is zero mean and preserves coordinatewise residual scale and per-horizon, per-channel spatial power, but not cross-channel or cross-horizon covariance. Sending the same field through both routes retains their induced cross-route covariance. The identifying assumption is approximate translation stationarity on the regular grid ([2], [51]). Using only previously committed moments isolates the probe graph and assumes that parameters vary slowly relative to the EMA window.

How Sampler Outputs Enter Internal-DW

All constructions above provide the same object: a $K$-step output-noise trajectory. Internal-DW passes that trajectory through the same fully open backward graph, so only the assumed noise distribution changes across samplers. Equations 17, 18, and 34 give the corresponding distribution $Q_d$. For every module, a draw $\widetilde{\epsilon}^{(d)}_{i,r,1:K}\sim Q_d$ is inserted into $q_{\mathrm{noise}}$ in (15); the resulting identity/nonlinear VJPs determine $\widehat{R}_{k,\ell}$ through (7). The route probe, PSD projection, and constrained Wiener solve are shared. Blocked out-of-fold prediction can reject the declared sampler assumption when its residual remains predictable, but it cannot prove that the residual is physically irreducible noise.

Why a Noise Sampler Is Needed and How Its Errors Affect the Gains

Proof of Route-Law Non-Identifiability

(4) depends on the decomposition $T=P+R$. Fully open route gradients identify the total moment $T$, but not its signal and noise components separately. Proposition 1 shows that this ambiguity is structural rather than a finite-sample defect of the estimator.

Proposition 1 (The open route law does not identify the split). Let $T\succ0$ be the covariance of an observed two-route vector $z$. For any two distinct matrices $R_1,R_2$ satisfying $0\preceq R_i\preceq T$, define independent Gaussian pairs

Then $z_1$ and $z_2$ have the same complete distribution $\mathcal{N}(0,T)$, while their unconstrained Wiener gains $w_i=T^{-1}(T-R_i)\mathbf{1}$ differ whenever $(R_1-R_2)\mathbf{1}\ne0$. Hence neither the total covariance nor, in general, the complete law of the open route vector identifies its signal and innovation covariances without additional structure.

Proof of Proposition 1. For each $i\in\{1,2\}$, independence gives

Because $s_i$ and $n_i$ are Gaussian, their sum therefore satisfies $z_i \sim\mathcal{N}(0,T)$. Thus $z_1$ and $z_2$ have the same complete distribution despite having different signal–noise decompositions.

For decomposition $i$, the signal moment is $P_i = T - R_i$, so the unconstrained Wiener gain is

Consequently,

which is nonzero whenever $(R_1 - R_2)\mathbf{1} \ne0$. Hence the same open-route law can imply different Wiener gains.

No amount of additional sampling from the open-route distribution resolves this ambiguity, because the indistinguishable decompositions agree at the population-law level. Classical innovation-covariance estimators resolve it only after imposing state-space, observation, repeated-measurement, or covariance-rank assumptions ([32], [6], [37]). Internal-DW makes the analogous domain assumption through $Q_d$ and transports sampled output innovations through the fully open VJPs to estimate the required route-noise covariance.

When a Fixed Sampler Recovers the Intended Noise

For residual-based samplers, we define the population residual targeted by the sampler and state conditions under which its empirical noise moment converges to that target. For domain $d$, let $c_t^{(d)}$ contain the observations available to the noise predictor, and let $\mathcal{G}_d$ be its predeclared predictor class. The best predictor in this class and its residual are

If $\mathcal{G}_d$ contains the conditional mean, this residual is the conditional innovation. Otherwise it also contains predictable structure that the chosen predictor class cannot represent. The guarantee below concerns this explicitly defined residual; it does not claim that model error is physical noise.

Proposition 2 (Consistency under a fixed sampler model). Fix a domain $d$ and basis $B_d$. Suppose that (i) training units are stationary and ergodic, with cross-fitting performed over independent units or asymptotically separated mixing blocks; (ii) the foldwise predictor satisfies $\widehat{q}_d^{(-f)} \to q_d^*$ in $L_2$; (iii) $B_d$ is fixed and bounded; and (iv) the empirical covariance is consistent in the retained finite-dimensional basis. Then the covariance of the out-of-fold residuals converges to

For a fixed bounded linear VJP map applied to the centered retained-basis residuals, the empirical-average route-noise moment $\widehat{R}_{k,\ell}$ also converges to its population counterpart.

Proof. Assumption (ii) makes each cross-fitted residual converge in $L_2$ to $\epsilon^{(d)}$. Boundedness of $B_d$ and the ergodic law of large numbers then give convergence of its first two moments; assumption (iv) gives the covariance limit. For a fixed linear VJP map, each entry of the transported second-moment matrix is a linear function of the residual covariance, so convergence follows directly.

The proposition separates finite-sample estimation error from sampler-model mismatch. More data reduce the former but cannot repair omitted conditioning variables or an inadequate predictor class, because these change the population residual being estimated. For covariance-based samplers, analogous consistency requires consistent estimation of the moments retained by the declared covariance model; misspecification of that model remains a population-level error.

How Sampler Error Changes the Estimated Gains

Why drive-dominated routes can favor the open operator. For one routed direction, write the open message as $d = s + n$, where $s$ is its conditional signal, $n$ is orthogonal noise, $P = \mathbb{E}\lVert s\rVert_2^2$ and $R = \mathbb{E}\lVert n\rVert_2^2$. Applying a scalar gain $c$ gives

Differentiation gives $c^\star=P/(P+R)$. Hence adding useful drive-dependent signal moves the correctly estimated gain toward the open route; degradation requires the applied gain to remain too small.

How moment error perturbs the two-route gain. Let $\widehat{T}=T+\Delta_T$ and $\widehat{R}=R+\Delta_R$, and suppose the unconstrained population solution

lies in the interior of $[0,1]^2$, so that $w^\star=w_{\mathrm{unc}}$. When the perturbations are small enough that the box constraint and PSD projection remain inactive,

This expansion follows by substituting

into the unconstrained plug-in solve and retaining the first-order terms.

Sampler-model mismatch enters through systematic error in $\Delta_R$, while finite task and noise probes contribute estimation error through both $\Delta_T$ and $\Delta_R$. If $b_w=\mathbb{E}[\widehat{w}]-w_{\mathrm{unc}}$ and $V_w=\operatorname{Cov}(\widehat{w})$, then

Thus moment error contributes to excess local route risk through both gain bias and gain variance. This remains a local route-risk statement and does not imply an end-to-end forecasting guarantee.

Controlled Known-SNR Study

This appendix specifies the oracle used in Figure 4. It lets us measure the conditional-mean gradient and innovation noise separately; ordinary data reveal only their sum. The experiment is a frozen-checkpoint probe: it performs no optimizer step and does not modify the model or its forward rollout. Any route gains saved in the checkpoint are overridden by ones, so every $g_k^{(b)}$ below is computed on the fully open BPTT graph. The gradient, gain-identification, and route-risk panels use the same frozen full BPTT checkpoint described below.

Separating Conditional Gradient Signal and Noise

The eight-dimensional process is

with independent Gaussian innovations $\xi_t\sim\mathcal{N}\left(0,\operatorname{diag}(1-a_1^2,\ldots,1-a_8^2)\right)$, where $a_j=(A_{\mathrm{AR}})_{jj}$. The variance choice makes every stationary coordinate have unit marginal variance while retaining distinct positive and sign-alternating memory scales. For an observed state $x_t$, the simulator gives the conditional future exactly:

All future-noise realizations in one comparison begin at this same observed state, so the model prediction and its tangent map are held fixed while only the target innovation changes. Equation 42 is written in raw simulator coordinates; the realized target and its conditional mean are both passed through the checkpoint’s same affine training normalization before the loss is differentiated.

Let $\mathcal{L}_k^{(b)}=K^{-1}\lVert\widehat{x}_{t+k}-x_{t+k}^{(b)}\rVert_2^2$ be the horizon loss for future-noise draw $b$, and let $g_k^{(b)}=\nabla_\theta\mathcal{L}_k^{(b)}$ be its complete full BPTT parameter gradient, including every path from the horizon loss back to the shared parameters. Because the prediction and its derivative are fixed under the conditioning and the loss is quadratic,

Thus the clean full-window diagnostic target is $\mu= \sum_{k=1}^{K} s_k$.

The per-horizon reliability quantities include the SNR in the left panel of Figure 4 and the noise fraction reported in the main text:

The first panel plots total-gradient RMS on the same parameter coordinates at every horizon, normalized by its value at $k = 1$, while (44) reports the signal–innovation split directly.

Risk of Including More Forecast Steps

Define $S_k = \sum_{j=1}^{k} s_j$ and $N_k^{(b)} = \sum_{j=1}^{k} n_j^{(b)}$. The prefix update is $S_k + N_k^{(b)}$, whereas the clean target remains the full-window signal $\mu= S_K$. Consequently,

The cross term vanishes by conditional centering. Full BPTT is the endpoint $k = K$: it has zero missing-signal bias but need not have minimum total risk. This equation is the formal meaning of the middle panel; its finite Monte Carlo closure error is at most 0.082 of the full BPTT risk in the four reported splits.

Probe protocol. The first two panels use $K = 32$, 15 independent test trajectories with one rollout start per trajectory, 64 innovation draws per conditional state, and four independent replicates. All curves within a replicate use the same realized gradients.

Gain Accuracy and Held-Out Route Risk

The third and fourth panels of Figure 4 use a frozen seed-0 full BPTT checkpoint at $K = 32$ and depth four. The diagonal-AR(1) module is fit from 68,541 training transitions without using the true process parameters. Its process-only diagnostic holds the analytic signal moment fixed, whereas the complete pipeline also estimates the finite-sample total route moment. Both are compared with gains obtained from the true process on disjoint calibration routes, producing 128 forecast-step–layer comparisons.

For the held-out risk audit, all gains are fixed before evaluating 15 untouched test trajectories at eight rollout starts with 64 process draws per start. The placement control circularly shifts the estimated gain pairs by half the horizon range within each layer, preserving each layer’s gain multiset. The reference given true $P$, $R$ is calibrated away from the test routes and then transferred to them; it is therefore not a test-conditioned mathematical oracle. Risk is $(w - 1)^{\mathsf{T}}P_{\mathrm{test}}(w - 1) + w^{\mathsf{T}}R_{\mathrm{test}}w$, pooled over all routes and normalized by the fully open operator. Ordinary and route-recording forward paths are bitwise identical for every arm.

Datasets and Training Details

This appendix records the datasets and configurations behind the selected-estimator comparison in Section 5.

Datasets

Table 2 summarizes the dataset sizes after preprocessing. The counted units differ across datasets; their construction and split protocols are described below. iEEG comprises 29 recordings from 16 participants; counts are pooled across participants, with models trained separately for each participant. State dimensions exclude the 512-dimensional CLIP visual inputs. Shear flow uses 32/4/4 source trajectories; clips overlap within each trajectory but never cross trajectories or splits.

Dataset	Counted unit	Train / Val. / Test	State shape	Length
MG	Trajectory	28 / 6 / 6	8	2048
NARMA-5	Trajectory	28 / 6 / 6	8	2048
iEEG	Chunk	464 / 87 / 87	40 or 80	1024
Movie fMRI	Subject	111 / 24 / 24	400	284
ETTm1 / ETTm2	Segment	135 / 45 / 45 each	7	256
Shear flow	Clip	288 / 36 / 36	$4 \times64 \times128$	64
WeatherBench-2	Segment	205 / 22 / 22	$29 \times121 \times240$	64

Mackey–Glass. For each coordinate, Mackey–Glass (Mackey & Glass, 1977) follows

with $\tau= 30$, $\beta= .2$, and $\gamma= .1$. We generate eight independent coordinates with identical dynamics and independent initial histories. The equation is integrated by RK4 (fourth-order Runge–Kutta method) with an internal step of .1, sampled every unit time, and run for a 1000-step burn-in before retaining each trajectory. We generate and cache 40 trajectories of length 2048 using a fixed data-generation seed. For each matched experimental seed, trajectories are randomly divided into 28/6/6 training/validation/test trajectories. Output states are standardized using the mean and standard deviation of the corresponding training trajectories only.

NARMA-5. We use a bounded variant of the NARMA benchmark (Atiya & Parlos, 2000). For each coordinate, bounded NARMA-5 follows

The eight coordinates use independently sampled observed input channels, and the first 200 steps are discarded as burn-in. For both datasets, we generate and cache 40 trajectories of length 2048 using a fixed data-generation seed. For each matched experimental seed, trajectories are randomly divided into 28/6/6 training/validation/test trajectories. Output states are standardized using the mean and standard deviation of the corresponding training trajectories only; the NARMA inputs retain their original scale.

iEEG. The iEEG experiment uses 29 movie-encoding recordings from 16 participants in the dataset of [21]. The neural state consists of macro-contact theta-power signals: 15 participants have 80 channels, and CS62 has 40. Each participant also has aligned 512-dimensional CLIP visual features, supplied as stimulus inputs.

For each participant, we split the movie timeline into 70/15/15 training/validation/test intervals. Repeated viewings use the same movie-time boundaries, preventing the same movie segment from appearing in different splits through different recordings. Repeated recordings are first truncated to their common length. We exclude 256 steps (5.12 s) immediately before each split boundary, not on both sides; no additional 2-s edge trimming is applied. The existing FIF signals preserve recording-wide Morlet z-scoring and noncausal filtering. Zero-phase FIR anti-alias downsampling also precedes the split. Therefore, preprocessing is not entirely split-local or strictly causal. At the final sampling rate of 50 Hz, each model step represents 20 ms. The existing 512-dimensional visual features are sampled at 25 frames/s and aligned by holding the most recent movie frame; audio is not included.

Each processed recording segment is divided into non-overlapping 1024-step chunks (20.48 s each), discarding incomplete trailing chunks. Chunks never cross recording or split boundaries. Multiple recordings from one participant are pooled only within their corresponding splits. Across participants, this yields 464/87/87 training/validation/test chunks, or 475,136/89,088/89,088 retained model steps. These are post-processing counts, not raw sample counts or independent participant counts.

Neural channels and CLIP feature dimensions are standardized separately for each participant using the mean and standard deviation computed from that participant’s training chunks only. The same statistics are applied unchanged to validation and test chunks; no statistics or training data are pooled across participants. Constant dimensions use a unit scale. This training-only claim applies to the second-stage normalization, not to the upstream processing described above. Pooled counts are given in Table 2. Forecasting uses $K = 64$ and test horizons 1–96 (up to 1.92 s); reliability diagnostics use horizons 1–64.

Movie fMRI. Movie fMRI uses the HCP 7T MOVIE1 AP acquisition [53]. We start from the HCP-provided tfMRI-MOVIE1_7T_AP_hp2000_clean.nii.gz volumes in MNINonLinear space, which have 1.6-mm isotropic resolution and a 1-s TR. The 2-mm Schaefer2018 400-parcel, seven-network atlas [46] is resampled to the HCP volume grid using nearest-neighbor interpolation. BOLD values are then averaged within each parcel, and each parcel time series is z-scored over time within the run.

Each subject sequence contains 20 pre-movie TRs, 244 MOVIE1 TRs, and 20 post-movie TRs. The accompanying movie-derived CLIP representation contains 256 1664-dimensional features per TR; these are averaged over the 256 tokens before entering the forecasting model. The packaged data record a six-TR hemodynamic alignment lag, and no additional temporal shift is applied by the dataset loader. Of the 159 subjects with complete inputs, a seed-specific shuffle assigns 111/24/24 subjects to the training/validation/test splits, so no subject appears in more than one split.

ETTm1 and ETTm2. ETTm1 and ETTm2 contain 15-minute measurements from two electricity-transformer records [58]. Each state contains oil temperature and six load measurements, and all seven coordinates are forecast. The observed future drive consists of six deterministic calendar coordinates: sine and cosine encodings of hour of day, day of week, and day of year.

We use the standard chronological partition consisting of the first 12 months for training, the next four months for validation, and the following four months for testing. Each split is divided into non-overlapping 256-step units, producing 135/45/45 training/validation/test units for each dataset. State and calendar coordinates are standardized separately using per-coordinate statistics computed from the training units only.

Shear flow. Shear flow uses the supplied train/validation/test partitions of shear_flow from The Well (Ohana et al., 2024) without reshuffling trajectories. Each state contains the scalar tracer and pressure fields and the two components of the velocity field. Scalar and vector fields are aligned and flattened into four channels before deterministic spatial subsampling by a factor of four, giving $64 \times128$ spatial inputs in our experiments. Within each upstream trajectory, we form 64-frame sequences with stride 16; sequences never cross trajectories or dataset splits. We apply no temporal subsampling, supply no external drive, and apply no additional dataset-level standardization. We use the $Re = 10^{5}$, $Sc = 2$ configuration, containing 32/4/4 training/validation/test trajectories of 200 frames each. With clip length 64 and stride 16, each trajectory yields nine clips, giving 288/36/36 clips in total. Clips overlap within a trajectory, but no trajectory is shared across splits.

WeatherBench-2. WeatherBench-2 [43] is extracted from the official 1.5° ERA5 archive at six-hour resolution on its 121 × 240 grid. The state contains geopotential, temperature, zonal wind, meridional wind, and specific humidity at 1000, 850, 700, 500, and 250 hPa, together with surface pressure, 2-m temperature, and 10-m zonal and meridional winds, for 29 channels in total. The observed future drive consists of sine and cosine encodings of hour of day and day of year.

We use 2007–2015 for training, 2016 for validation, and 2017 for testing, corresponding to 13,148/1,464/1,460 time steps. Per-channel means and standard deviations are computed over all training times and spatial locations and then applied unchanged to every split. No additional spatial subsampling is performed. Each split is divided into non-overlapping 64-frame segments, producing 205/22/22 training/validation/test samples.

Implementation Details for Training

Within each dataset, full BPTT and Internal-DW use the same split, initialization seed, forward model, task loss, optimizer update, validation rule, and evaluation code. Full BPTT keeps every route open; activation checkpointing is used where needed to reproduce the same gradient with lower memory. Internal-DW changes only the backward operator and adds the calibration probes described in Appendix C.

Settings shared across methods. For iEEG, full BPTT and Internal-DW are trained separately for each participant using three matched seeds (0, 1, and 2). Model parameters and training data are not shared across participants; this evaluates within-participant forecasting, not generalization to unseen participants. All matched arms use Adam, train for at most 100 epochs with validation patience 20, and report the minimum-validation-loss checkpoint. The primary arms use global gradient-norm cap 1; stronger clipping is stated with the corresponding control below. Sequence models use four residual Mamba blocks ([16]) with state dimension 16, convolution width 4, expansion factor 2, a linear read-in, and a two-layer read-out, and predict an increment added to the current state. They use a 16-step input window, a 32-step burn-in, and 16 rollout starts. Field models use a depth-4 residual U-Net with channel multiplier 2 and a two-frame input.

Data	Model	$K$/loss	LR/batch	Noise sampler	W/C
MG, NARMA	Mamba-128	32/Rel-$L_2$	$10^{-4}$/32	Generic	8/4
ETTm1/2	Mamba-128	64/Rel-$L_2$	$10^{-4}$/32	Generic	8/4
iEEG	Mamba-256	64/Rel-$L_2$	$10^{-4}$/32	factor Prior	8/4
fMRI	Mamba-4096	64/MSE	$10^{-4}$/32	subject-AR Prior	8/4
Shear	U-Net-32	32/Rel-$L_2$	$3 \times10^{-4}$/4	spectral Prior	8/4
WB2	U-Net-32	48/Rel-$L_2$	$10^{-4}$/8	spectral Prior	2/1

The sequence and WeatherBench-2 runs use weight decay $10^{-4}$ and a step-100 schedule; shear uses zero weight decay and cosine decay to $10^{-6}$. The gains start fully open, with route- and residual-moment EMA coefficients .95 and .99. The Clip control uses cap .1, and JReg uses penalty .1, target norm 1, and finite-difference scale .001.

Selection of TBPTT, Static, and DW variants. TBPTT segment length is selected on seed 0 and then fixed across three seeds: the candidate sets are $\{8,16\}$ for Mackey–Glass and shear and $\{16,32\}$ for ETTm1/2, with selected lengths 8 and 16. The ETT TBPTT runs additionally use cap .1. The Static control ties every routed merge to $\alpha= m = c$, selects $c \in\{.3,.6,.9\}$ on seed-0 validation loss, and then fixes it across seeds. For each dataset, we select one DW-Generic and DW-Prior using validation loss and then fix that choice for test evaluation. Calibration cost and update order are given in Appendix C.5.

Movie-fMRI sampler fitting. The fMRI noise sampler is fitted separately for each seed using only its training subjects. Five-fold subject cross-fitting estimates a time-aligned shared response from the fitting subjects and a coordinatewise ridge-regularized AR model for subject-specific deviations. The AR order is selected from $\{0,1,2,4,8\}$ by cross-fitted multi-horizon MSE and is 8 for all three seeds. Length-64 residual trajectories from the held-out training folds are retained intact as noise templates, preserving cross-parcel and cross-horizon dependence.

How History and Observed Drive Are Measured

Computing the Two Regime Scores

Let $R_{\mathrm{mean}}$ be the held-out risk of a training-mean predictor, $R_D$ that of a predictor using the observed future drive, and $R_{D+H}$ that of the same predictor augmented with history. Figure 3 uses

and reports $\Delta\Psi_H = \Psi_H - \Psi_H^{\mathrm{null}}$, where the null uses circularly shifted history.

Figure 3 first averages the values across linear and nonlinear probes and then clips negative values to zero because the map uses only detected predictive value. Table 4 reports the corresponding unclipped scores.

Data	Linear probe		Nonlinear probe
	$\Psi_D$	$\Delta\Psi_H$	$\Psi_D$	$\Delta\Psi_H$
Mackey–Glass	.000	.442	.000	.370
NARMA-5	.857	−.001	.283	−.001
iEEG	−3.2823	.0770	−.1265	.0032
Shear flow	.000	.859	.000	.544
WeatherBench-2	.219	−.179	.232	.004
ETTm1	−.341	.697	−.355	.110
ETTm2	−.144	.963	−.019	.272
Movie fMRI	.230	.039	.230	.016

Interpreting the Regime Scores

The raw scores measure held-out risk reduction and are not clipped at zero, so they can be negative. $\Psi_D < 0$ means that the drive-based probe performs worse than the training-mean predictor, while $\Delta\Psi_H < 0$ means that aligned history performs worse than the shifted-history null. We interpret either outcome as no detected predictive value under the probe, not as negative information.

For autonomous testbeds, the drive-only predictor is the training mean, so $\Psi_D = 0$ by construction. The drive predictor uses observed future drive for NARMA and WeatherBench-2, calendar encodings for ETTm1 and ETTm2, and a subject-disjoint estimate of the movie-locked shared response for fMRI. A linear ridge probe and a nonlinear local-analog probe select history width on validation data and are evaluated once on held-out horizons $k \ge8$.

Validation of the Predictive-Regime Probes

We perform a controlled Mackey–Glass sweep to verify that the predictive-regime probes respond to increasing observed drive.

For each of the eight independent coordinates, the driven dynamics are

where $\tau= 30$, $\beta= .2$, and $\gamma= .1$. The observed drive is constant within each unit sampling interval, $u(t) = u_j$ for $t \in[j,j+1)$, with

Thus, $\lambda_D$ controls the amplitude of the additive, unit-variance AR(1) drive; $\lambda_D = 0$ recovers the autonomous system. The future drive is supplied to the probes as observed input.

We hold the autonomous dynamics, underlying realizations, and data splits fixed, and vary only the drive scale $\lambda_D$, showing settings from .03 to .34. The corresponding magnitudes of autonomous and driven terms are shown in Table 5. At each setting, we compute the drive-only gain $\Psi_D$ and additional history gain $\Delta\Psi_H$ defined above, using long-horizon averages over $k \ge8$. Table 5 also reports the linear and nonlinear probe scores separately, before averaging or clipping. For the plotted coordinates in Figure 8, the linear and nonlinear probe estimates are first averaged, after which each coordinate is clipped at zero, following the construction used in the main predictive-regime figure.

	Term magnitudes				Linear		Nonlinear
$\lambda_D$	RMS($A$)	RMS($D$)	$\frac{\mathrm{RMS}(D)}{\mathrm{RMS}(A)}$	$\Psi_D$	$\Delta\Psi_H$	$\Psi_D$	$\Delta\Psi_H$
0.03150	0.0420	0.0315	0.75	0.023	0.802	0.040	0.522
0.05925	0.0591	0.0592	1.00	0.144	0.757	0.135	0.533
0.08000	0.0724	0.0799	1.10	0.339	0.680	0.283	0.463
0.13000	0.1039	0.1299	1.25	0.645	0.568	0.504	0.362
0.34000	0.2424	0.3398	1.40	0.930	0.598	0.710	0.086

Figure 8 shows a rightward movement over the displayed range: as $\lambda_D$ increases from .03 to .34, the measured drive-only gain $\Psi_D$ increases from approximately .03 to .82. Meanwhile, $\Delta\Psi_H$ decreases while remaining positive, suggesting that the additional predictive value of history beyond the observed drive weakens. Thus, the probes respond in the intended direction to increasing observed drive while still detecting additional history value.

Controlled Test of Strong Observed Drive

We isolate the effect of observed drive by comparing the autonomous Mackey–Glass testbed with a variant that retains the same delay dynamics and adds a unit-variance AR(1) drive with coefficient .9 and scale .08. For each of the eight independent channels, the observed forcing is an independent stationary AR(1) process with coefficient .9 and unit variance, generated using Gaussian innovations. It enters the Mackey–Glass dynamics additively with scale .08 and is held constant within each unit sampling interval. We retain the autonomous parameters $\tau=30$, $\beta=.2$, $\gamma=.1$, and exponent $10$, and integrate the dynamics using RK4 (fourth-order Runge–Kutta method) with step size .1 after a 1000-step burn-in. We generate 40 trajectories of length 2048 and split them into 28/6/6 training, validation, and test trajectories. The eight future forcing channels are provided to the forecaster as observed future drive. The driven variant continues to have substantial long-horizon history value, but now also has a strong observed drive. Both conditions use $K = 32$, the same matched training protocol, and dense 1:48 relative-$L_2$ evaluation.

condition	Full BPTT	Internal-DW	DW vs. Full
weak drive (autonomous)	$1.05412 \pm.00436$	$.99966 \pm.00958$	5.17% lower
strong observed drive	$.73489 \pm.01064$	$.79157 \pm.00350$	7.71% higher

Internal-DW improves the autonomous condition but is worse on every seed after the observed drive is strengthened. This controlled comparison shows that substantial history alone does not guarantee improvement and provides a concrete stress test of the current estimator’s drive-related boundary.

How Forecasting Error Is Evaluated

Let $\mathcal{U}$ denote the held-out test units and $\mathcal{O}_u$ the forecast origins selected within unit $u$. We roll the model out once from each origin without resets through $H_{\mathrm{eval}} = \lceil1.5K \rceil$ and compute (9). We use at most 64 deterministic evenly spaced origins per test unit, remove duplicate absolute origins, average origins within a unit, and then weight units equally. The sole headline score averages relative $L_2$ over 1:$H_{\mathrm{eval}}$; trained-window, post-window, and terminal summaries are not reported as separate performance criteria.

How Per-Horizon Gradient Utility Is Measured

This appendix gives the protocol behind fig:5a. At a fixed seed-0 full BPTT checkpoint, we pair a calibration example $A$ from the training split with a disjoint example $B$ from the test split, using at most one common-relative-position rollout start from each loader item. For every forecast step $k$, we compute the parameter gradient $g_k^A$ produced by that step’s loss and the full-horizon held-out BPTT gradient $\bar{g}^B$. All horizons and pairs within a dataset use one fixed sample of parameter coordinates with the corresponding full-dimensional scaling. We summarize each horizon by its relative gradient norm and held-out utility:

Here, $A(k)$ measures relative gradient magnitude, while $U(k) = \langle g_k^A/\lVert g_k^A \rVert_2,\bar{g}^B\rangle$ measures alignment with the held-out full-horizon gradient. Under the infinitesimal update $-\epsilon g_k^A/\lVert g_k^A \rVert$, the held-out rollout loss changes by $-\epsilon U(k) + O(\epsilon^2)$; hence $U(k) < 0$ means that the training direction raises held-out loss. The black curve and band in the figure are the median and 10–90% range of $A(k)$, while the purple curve and band are the median and interquartile range of $U(k)$. The displayed percentage is the fraction of forecast steps whose median $U(k)$ is negative. The numbers of disjoint train–test pairs are 6 for Mackey–Glass and NARMA-5, 8 for ETTm1 and ETTm2, 3 per single-recording iEEG participant and 6 per two-recording participant, 4 for movie fMRI and shear flow, and 8 for WeatherBench-2. The iEEG probe uses seed 0 for each of 16 participants. Within-participant medians are aggregated by taking the median across participants at each horizon. Its gray and purple bands are respectively the 10th–90th and 25th–75th percentiles across participant curves. The 45/64 negative-horizon annotation is computed from this displayed median, not by averaging participant-specific negative fractions.

Dense-Horizon Relative-$L_2$ Results

Table 7 reports the dense-horizon test relative $L_2$ values underlying Figure 6. The values are computed using the same evaluation horizons and forecast origins and are reported as mean $\pm$ standard deviation over three matched seeds; lower is better. Because the normalization is dataset-specific, values are intended for within-dataset rather than across-dataset comparison.

Method	MG	ETTm1	ETTm2	Shear flow
Full BPTT	1.05412 ± 0.00436	0.94220 ± 0.02519	0.86842 ± 0.09126	0.16386 ± 0.03538
Clip	1.05923 ± 0.00813	0.93199 ± 0.02141	0.83840 ± 0.02873	0.19576 ± 0.02088
JReg	1.05216 ± 0.00583	0.96169 ± 0.00679	0.86930 ± 0.09105	0.21068 ± 0.01514
TBPTT	1.03440 ± 0.01724	0.96533 ± 0.02596	0.85733 ± 0.06582	0.13222 ± 0.00455
Static	1.00467 ± 0.00727	0.92499 ± 0.01133	0.80427 ± 0.06864	0.23296 ± 0.00567
Internal-DW	0.99966 ± 0.00958	0.87864 ± 0.01758	0.74894 ± 0.04886	0.14523 ± 0.02466
Method	NARMA-5	iEEG	Movie fMRI	WeatherBench-2
Full BPTT	0.87065 ± 0.01501	1.88713 ± 0.03089	0.92206 ± 0.01092	0.61112 ± 0.00052
Clip	0.80678 ± 0.00240	1.91407 ± 0.02655	0.93315 ± 0.00835	0.62510 ± 0.00143
JReg	0.87313 ± 0.01149	1.88745 ± 0.03118	0.96189 ± 0.07436	0.61407 ± 0.00382
Internal-DW	0.93908 ± 0.01382	1.87402 ± 0.03155	0.98418 ± 0.02106	0.61169 ± 0.00294

End-to-end training time

We measure the wall-clock cost of the assigned full BPTT and Internal-DW training configurations on an NVIDIA H200 NVL GPU. Within each matched pair, the two arms use the same dataset, model, training horizon, batch geometry, optimizer, and physical GPU. We discard one warm-up epoch (iEEG uses 40 warm-up training batches) and use the median time of the next three epochs. Table 8 reports the mean and sample standard deviation of these within-run medians over three paired repeats; the execution order is alternated across repeats to reduce order and thermal effects. CUDA is explicitly synchronized at both timing boundaries. The timer covers the training loop only and therefore excludes validation and checkpoint I/O.

Data	Full BPTT (s/epoch)	Internal-DW (s/epoch)	Increase (%)
Mackey–Glass	21.28 ± 0.09	27.43 ± 0.42	28.9 ± 1.7
ETTm1	27.78 ± 0.06	35.32 ± 0.23	27.1 ± 0.9
ETTm2	28.68 ± 0.08	36.31 ± 0.31	26.6 ± 1.0
Shear flow	29.79 ± 0.43	41.55 ± 0.80	39.5 ± 0.7
NARMA-5	22.26 ± 0.13	28.31 ± 0.15	27.2 ± 0.2
iEEG (per subject)	31.64 ± 0.03	44.83 ± 0.39	41.7 ± 1.2
Movie fMRI	498.67 ± 1.98	599.90 ± 3.04	20.3 ± 0.1
WeatherBench-2	32.84 ± 0.39	43.38 ± 0.39	32.1 ± 0.9

The Internal-DW measurements include its online covariance probes, routewise gain solves, and application of the resulting gains. Construction of a train-only domain-prior artifact, when used, is a one-time preprocessing step and is excluded. Thus, this is an end-to-end comparison of the training implementations used in our forecasting experiments rather than an isolated solver microbenchmark.

Across the eight testbeds, the mean paired increase in epoch time ranges from 20.3–41.7%. This additional training cost reflects the online estimation and routewise gradient scaling that distinguishes Internal-DW from full BPTT. The forward model is unchanged, and Internal-DW introduces no inference-time computation.

Per-seed fitted curves

Figure 9 decomposes the averaged training-horizon trends in Figure 7 into individual seeds. For each dataset, seed, and method, we independently fit a least-squares quadratic constrained to have non-negative curvature. A boundary solution may reduce to a linear trend, in which case the fitted optimum lies at an endpoint of the evaluated range. The gray and blue vertical dashed lines indicate the fitted optima for Full BPTT and Internal-DW, respectively.

The Internal-DW optimum lies at or to the right of the Full-BPTT optimum in all 12 dataset–seed comparisons. The shift is strict in 11 comparisons, while ETTm1 seed 0 gives the same boundary optimum for both methods. Thus, the shift toward a longer fitted optimal training horizon is not produced by averaging across seeds.