Tuning the burn-in phase in training recurrent neural networks improves their performance
Authors: Julian D. Schiller
Malte Heinrich
Victor G. Lopez
Matthias A. Müller
Venue: ICLR 2026 Poster (2026)
Keywords: RNN
training
learning theory
optimization
Primary area: learning theory
TL;DR: We establish theoretical bounds on the performance loss when using truncated sequences for RNN training and show how they improve when choosing a suitable burn-in phase.

Abstract: Training recurrent neural networks (RNNs) with standard backpropagation through time (BPTT) can be challenging, especially in the presence of long input sequences. A practical alternative to reduce computational and memory overhead is to perform BPTT repeatedly over shorter segments of the training data set, corresponding to truncated BPTT. In this paper, we examine the training of RNNs when using such a truncated learning approach for time series tasks. Specifically, we establish theoretical bounds on the accuracy and performance loss when optimizing over subsequences instead of the full data sequence. This reveals that the burn-in phase of the RNN is an important tuning knob in its training, with significant impact on the performance guarantees. We validate our theoretical results through experiments on standard benchmarks from the fields of system identification and time series forecasting. In all experiments, we observe a strong influence of the burn-in phase on the training process, and proper tuning can lead to a reduction of the prediction error on the training and test data of more than 60% in some cases.

Published as a conference paper at ICLR 2026
TUNING THE BURN-IN PHASE IN TRAINING RECURRENT NEURAL NETWORKS IMPROVES THEIR PERFORMANCE
Julian D. Schiller
Malte Heinrich
Victor G. Lopez
Matthias A. Müller
Leibniz University Hannover, Institute of Automatic Control, Hannover, Germany
ABSTRACT
Training recurrent neural networks (RNNs) with standard backpropagation through time (BPTT) can be challenging. especially in the presence of long input sequences. A practical alternative reduce computational and memory overhead is perform BPTT repeatedly over shorter segments of the training data set, corresponding truneated BPTT. In this paper, we examine the training of RNNs when using such truncated learning approach for time series tasks. Specifically, we establish theoretical bounds on the accuracy and performance loss when optimizing over subsequences instead of the full data sequence. This reveals that the burn-in phase of the RNN is an important tuning knob in its training. with significant impact on the performance guarantees. We validate our theoretical results through experiments on standard benchmarks from the fields of system identification and time series forecasting. In all experiments, we observe strong influence of the burn-in phase on the training process, and proper tuning can lead to a reduction the prediction error the training and test data of more than 60% in some cases.
INTRODUCTION
Recurrent neural networks (RNNs) are inherently suited for processing sequential data (Hewamalage al., 2021; Lipton al., 2015). For this reason, they have wide range of applications, including natural language modeling (Khan et al., 2023; Mikolov et al., 2011; 2013), text and image generation (Sutskever al., 2011; Gregor et al.. 2015), and computer vision (Buhrmester et al., 2021), as well as classical time series tasks such as modeling dynamical systems (Pillonetto et al., 2025; Beintema et al., 2023; Forgione& & Piga, 202 abb, Schussler et al., 2019: Wang & Chen, 2006) and time series forecasting (Hewamalage et al., 2021: Elsworth & Güttel, 2020; Bianchi al., 2017). While the Elman network (Elman, 1990) represents one of the most classic RNN architectures, prevailing approaches in practice combine the vanilla RNN with additional gating mechanisms, such Long Short-Term Memory (LSTM) (Hochreiter & Schmidhuber, 1997) and Gated Recurrent Units (GRU) (Cho et al., 2014). RNNs have recently become very active field of research again, driven also by the impressive performance that transformers have recently achieved for long-range reasoning (Kim al., 2025). Modern architectures include deep state-space models (SSMs) such as S4 (Gu et al., 2022) and Mamba (Gu & Dao. 2024; Wang et al.. 2025), as well as linear recurrent units (LRUs) (Orvicto al., 2023; Zucchet et al., 2023). Other recent works involve Koopman theory (Li et al., 2025), sparse deep learning (Zhang et al., 2023), and probabilistic methods (Coscia et al., 2025).
RNNs are traditionally trained using backpropagation through time (BPTT) (Werbos, 1990); namely, in an iterative fashion using steepest descent. This requires to perform, at each iteration, forward pass to compute the loss and backward pass to compute the gradients, each running over the whole training sequence. However, in case the training sequence is large, many numerical problems can arise: the loss function landscape becomes increasingly complex (Allen-Zhu et al., 2019; Ribeiro et al., 2020), the computations of the forward and backward pass become computationally intense, and they may even blow up due to vanishing or exploding gradients (Bengio et al., 1994; Pascanu al., 2013; Zucchet & Orvieto, 2024). A standard approach ensure numerical stability, efficiency, and faster training due to less computational and memory overhead is to simply truncate the forward and backward passes, realizing truncated backpropagation through time (TBPTT) (Williams & Peng.
*Corresponding author (c-mail: schiller iriiller@iit..uiii (ia.erriiiiiiimmmeeeeeeeell code is publicly available in the GitHub repository at https: / / coub.comm jdachi 11 11e/ /1e/rnn-eeeeeeeee wwe/rnnntraiinin rnn-training

---

Published as a conference paper at ICLR 2026
1990). This can be leveraged in modified, stochastic version of the parameter update law (i.e., using stochastic gradient descent, SGD), where stochasticity comes from the fact that the loss gradients computed over (a batch of) random segments of the training data set. Overall, TBPTT can be regarded as a practical compromise between learning effectiveness and computational effort. is a standard method for RNN training, particularly in the presence of long training sequences (Mikolov et al., 2011; Tang & Glass, 2018; Aicher et al., 2020; Xu et al., 2023).
Computing the forward pass in TBPTT requires initializing the hidden state, which is usually set to zero or chosen at random. This may hinder the learning of long-term dependencies. Alternatives include propagating the hidden state across iterations (Xu et al., 2023; Elsworth & Güttel, 2020; Katrompas & Metsis, 2021) or estimating it using additional trainable autoencoders (Mohajerin & Waslander, 2017; Masti & Bemporad, 2021; Forgione et al., 2022; Beintema et al., 2021; 2023). While beneficial when long-term dependencies are crucial, zero initialization prevails in applications because the implementation is much simpler and renders the overall training process more efficient, numerically stable, and robust with respect to local minima (Xu et al., 2023; Allen-Zhu et al., 2019).
In this work, we consider the training of RNNs for time series tasks, relying on TBPTT with zero initialization to ensure a simple implementation and efficient training. To cope with internal transients resulting from initialization, we introduce the burn-in phase of the network as an additional parameter of the training process, effectively excluding them from the loss function. This was also previously suggested by Jaeger (2002); Bonassi et al. (2022); Beintema et al. (2021), referring to it (inconsistently) as "wash out" or "burn time"; however, without analyzing its influence the resulting performance or providing systematic tuning guidelines, leaving it a rather unnoticed heuristic.
In this paper, we provide answers the following questions.
1. How much do we lose in the training by employing TBPTT with zero initialization?
2. How can we suitably modify the training process to reduce this loss?
We establish novel theoretical results on the performance of RNN training when optimizing over subsequences of the training data using truncated learning approach. More specifically, we estimate the performance gap (regret) with respect a challenging benchmark RNN, trained by optimizing its hidden states on all of the training data. Regret is a standard measure for assessing algorithms in reinforcement learning (Jaksch al., 2010; Agrawal al.. 2021), online optimal control (Agarwal et al., 2019; Li et al., 2019; Nonhoff & Müller, 2022; Martin et al., 2025), and state estimation (Brouillon et al., 2023; Schiller al., 2025). Our technical analysis based on viewing the training problem through the lens of optimal control, focusing in particular on the turnpike phenomenon. This property generally implics that optimal trajectories most the time stay close to an optimal time- varying reference, which is regarded as the turnpike, see, e.g., Trélat & Zuazua (2025); Faulwasser & Grüne (2022); Zaslavski (2006); Carlson et al. (1991) for further details on this topic. Turnpike- related control
arguments have generally proven useful for analyzing the performance of optimization- and estimation techniques (Faulwasser et al., 2018; Grüne & Pirkelmann, 2019; Schiller et al., 2025).
The derived regret estimates depend the hyperparameters of TBPTT, particularly the interplay between the internal forgetting rate of the network and the burn-in phase. These novel insights can be leveraged to select suitable burn-in phase during training, ensuring that the resulting performance is close to that of the benchmark RNN. We validate our theoretical results in experiments using real world datasets from the fields of system identification and time series forecasting. In all experiments, we observe a strong influence of the burn-in phase on the training process, and proper tuning can lead to reduction of the prediction error on the training and test data of more than 60% in some cases. The main findings of this paper for RNN training using TBPTT can be summarized as follows:
Theory: We develop the concept of the burn-in phase from a heuristic to a theoretically grounded method.
Practice: Tuning the burn-in phase can yield significant performance increase.
2 RECURRENT NEURAL NETWORKS
The core of any recurrent network is the internal dynamics of its hidden states. Given an input sequence { of certain length T' > 0, the evolution of the hidden state ht for some pre-
2

---

Published as a conference paper at ICLR 2026
1/1
3/3
202
274
RGG
9())
layers
3/1
3/2
3/x
23
2T
20
222
<33
2
Figure I: General structure of an RNN (left) with its internal architecture unfolded in space and time (middle). The network consists of L recurrent layers and processes input sequence of length T. Black arrows denote instantaneous spatial processing, red arrows denote temporal processing. Each arrow represents learnable linear or nonlinear projections and/or layers. potentially with additional skip connections, which can generally be captured by the functions I' and g. The right figure shows an exemplary RNN for time series forecasting applications, producing a forecast of F =3 future values of sequence (rx) ) using a lookback window of N = 4.
specified initial state ho can generally be described by the equation
h (he-11 0n), == -....,,
(I)
where f is function parameterized by the network parameters On E RR"., The output y/ of the network at time t is produced through an output layer, which can be characterized in terms of a function g such that
y/ = g(h,, me; 0y); 1= 1=1,.,,,
(2)
where By E By corresponds to the parameters of the output layer. For convenience, throughout
the following we consider the combined parameter vector which consists of all the learnable parameters of the model in (1) and (2). Here, want to emphasize that we consider a general multivariate setup, where x; E RR hr E Rdd,, and yt E R.d.
[0] 0_1 E IR" with n = 7, + Ny
Overall, the output y/ at time 1,..., T depends the initial hidden state hy, the input sequence {15})==1 and the model equations (1)-() for particular network parameter 0. To make this explicit, we sometimes express the output y/s in terms of the mapping ya (ho, 0, {-,}Ş=1> or, for convenience, using the full input sequence y4 (ha. 8, {2, } (which is equivalent, as (1) ensures causality).
Note that the model in (1)-2) constitutes a very general RNN description, which includes all com- monly used network architectures, such classical Elman RNNs (Elman, 1990), LSTMs (Hochreiter & Schmidhuber, 1997), GRUs (Cho et al., 2014), and echo state networks (Jaeger, 2002), but also recently proposed architectures such LRU networks (Orvieto al., 2023). In deep RNN structures, the hidden state ht corresponds to the stacked hidden states of all individual layers, and the function f embodies the dynamics of the individual layers together with possibly additional nonlinear projec- tions, gated linear units, and/or normalization operators between the layers, see also Figure 1. For the most simple case of a one-layer Elman RNN, the functions f and g and the parameter vector specialize to
f(h, x; 0) = o( W + WA. + bu),
xhh, = Why + by. (ven veo Un vec(Why.) _.
Here, the network dynamics (1) involve the weight matrix Wit for the recurrent connection, the weight matrix Wxh between the input and the hidden layer, the bias vector bn, and the activation function & (e.g., tanh or ReLU);; the output (2) is formed by a linear layer with weight matrix Why and bias vector by. All the weight matrices and bias vectors can be collected in the parameter vector O.
Note that we define the network output in (2) for all f = 1, ..., T. While this is natural in the context of dynamical systems, in many standard RNN applications such as classification tasks and time series forecasting, the output is usually only formed after the entire input sequence has been processed, that is, at time T. However, g still be evaluated earlier times 1,..., T - 1, see Figure 1. Finally, note that the RNN (1)-() covers any parameterized nonlinear state-space model; the subsequent analysis therefore directly also covers elassical system identification problems (L.jung, 1999).
3

---

Published as a conference paper at ICLR 2026
D1
m
171
D
M
D
Figure 2: Data segmentation in TPBTT learning. Depicted are the full data sequence D of length T (green) and four (out of S) subsequences of length (red). The subsequences D; and D;-1 start at samples and 8;-1, respectively. and result in the overlap o; compare (5). Blue highlighted segments appear in the loss function (0; D;) in (7), which is determined by the burn-in phase length m.
TRAINING RNNS USING TBPTT AND MINI-BATCH SGD
In this section, revisit the overall TBPTT learning algorithm to properly introduce our setup and notation. In particular, we illustrate the data segmentation process inherent to TBPTT (Section 3.1) I), define the segment-wise loss function (Section 3.2), show basic mini-batch training scheme (Section 3.3), and present a suitable measure of training performance (Section 3.4). In this work, we consider the standard mean squared error (MSE) loss, where we treat the burn-in phase of the RNN as an additional hyperparameter of the training process (see Section 3.2). In Section 4, we establish novel regret guarantees for truncated learning, which exhibit strong dependence on the burn-in phase.
3.1 DATA SEGMENTATION IN TBPTT
Let an input-output time series pair D = (x 2/ ) of length be given and denote the input data sequence as X := {= [==1" TBPTT essentially implements the BPTT algorithm on short segments of the training data set. To this end, we divide the training sequence D into S subsequences, each of which has length N < T. We define the index sets for the subsequences
S == {1,..., S} and for the elements within a subsequence as N {1,... N}. Moreover, we denote the input-output data points of subsequence E S that starts data sample 8; [1, N by 29 2 -
and 0) == Va for all j E N. For convenience, these are grouped as follows:
D, = {(=))) 1) XT (+ (1))) i E
In order make of all available data during training, in the following restrict ourselves the case where the subsequences D; cover the complete data set D, and furthermore, assume that the start samples 8;, i arranged in ascending order. Hence, subsequence D; may share some data points with the preceding subsequence D;-1, which can be characterized by the overlap o, as
0, N (s: 8;-1), 1), 2,... S.
(5)
Here, we note that o takes values in the interval [0, 1]. Considering all S subsequences, we define the minimum overlap according to
Omin min N 25:55
(s: 8;-1)-
25:55
The overall data segmentation process is also illustrated Figure The variables N, and s; i E S, represent important hyperparameters of the TBPTT algorithm that have a strong influence on the training process and the resulting performance the trained RNN. In general, suitable values may be chosen by leveraging additional information about the given data set and the underlying application (such as specific patterns or seasonal components of different length). If no such information is available, for fixed value of N, a practical choice is to select S = T N + 1 (i.e., use as many subsequences of length N as possible), which implies that 8; i for all E and yields the maximum overlap o N - Omin uniformly for all 2 < i <i<a S.
3.2 SEGMENT- WISE LOSS FUNCTION WITH BURN-IN PHASE
In this paper, we consider the MSE loss function a standard choice for time series tasks. Extensions other loss functions suitable for, e.g., classification problems interesting topic of future work.
4

---

Published as a conference paper at ICLR 2026
For partieular parameter 8, we define the loss function involving the data subsequence D, as
L(O: Dj) N
(0,, - u/ E
j=
The variable m appearing in the lower bound of summation in (7) refers to the burn-in phase of the network, see Figure 2. an additional design variable to compensate for network outputs that might be negatively affected by internal transients resulting from the zero initialization. It can be chosen from the interval [o, N V- 1], where m recovers the standard MSE loss (without burn-in phase). Similar approaches have been used previously (Jaeger, 2002; Bonassi et al., 2022; Beintema et al., 2021), albeit without performing any detailed analysis or providing tuning guidelines. Note that N 1 captures scenarios in which the output is generated by default only after processing the entire imput subsequence X; of length N, as is the case with time series forecasting, compare Figure (right). Nevertheless, in the following sections show both theoretically and experimentally that even here, a choice of 0 < m < should generally be considered a hyperparameter of RNN training.
can be beneficial for overall performance, suggesting that m
3.3 BASIC TRAINING ALGORITHM
The network parameter 0 are usually trained iteratively with SGD using mini-batches (Lipton et al., 2015; Hewamalage et al., 2021; Tian et al., 2023). More specifically, each iteration k, we choose random subset B3 C of size [Bx] =: b S to update 0 according
0k+1 02 - ~0k
AOA = b
E VoL(On; D.). IEB
Here, A0p represents the batch-averaged gradient. The learning 7is is usually chosen to be a small constant or adapted with time (Tian et al., 2023). The learning procedure is summarized in Algorithm 1. Here, Step 4 actually corresponds subsequence i contained in the current mini- batch B₁. forward pass to compute the segment- wise loss L(#μ; Dz), and a backward pass to compute the corresponding gradient. However, these computations are completely independent for each i 83 and can therefore be performed in parallel (possibly on a GPU), which significantly accelerates the training. The basic algorithm can be adapted, c.g., by employing SGD variants such as AdaGrad, RMSProp. Adam, etc. (Tian et al., 2023; Bottou & Bousquet, 2007; Kingma & Ba, 2015), all of which are part of standard deep learning libraries such as PyTorch. Recent results on the global convergence and stability of SGD algorithms in the context of machine learning can be found in the work Patel et al. (2022); Allen-Zhu et al. (2019).
an additional loop that performs, for each
Algorithm 1 RNN training using TBPTT with mini-batch SGD
Require: data set D. learning rate 176
1: Initialize 0g E IR", k = O.
2: while algorithm not converged do
3: Sample B3 C S
>> Select current mini-bater b Compute batch-averaged gradient Update parameter
4: A0, 5 Lies, (oI D,)
5: 0k+1 0/ n&AAR
6:
k=
7:
end while
Remark 1. Algorithm 1 employs zero initialization of each segment-wise forward pass, which is the most classical form of TBPTT. slight modification is to preserve the chronological order of the subsequences when selecting the current batch BR in Step 3 and pass the hidden states from iteration to k + 1, which is sometimes referred to as "stateful" training (Katrompas Metsis, 2021; Elsworth & Güttel, 2020; Xu et al., 2023), Although this can be advantageous when long-term dependencies are crucial, the former approach prevails in deep learning applications because the randomness of the training algorithm can help escape from local minima and generally renders the training process more numerically stable, robust, and efficient (Xu et al., 2023; Allen Zhu et al., 2019). Still, extending our theory stateful training could be an interesting topic for future work.
5

---

Published as a conference paper at ICLR 2026
3.4 TRAINING PERFORMANCE
Recall that TBPTT a practical approximation of BPTT using the full training data set D. Hence, seems meaningful measure the performance of the training algorithm by re-evaluating the network after successful training, processing the full training input sequence X and computing the MSE between the generated outputs and the ground truth data y Specifically, for given and initial hidden state ho. we define
Σ lvo 8, xxx 1/11?. P(ho, 8; B; T m t=++1
(9)
To ensure a fair comparison between parameter sets when using different network initializations, we in fact consider the truncated MSE in (9), where, similar to the segment- -iise loss in (7), we exclude internal network transients by incorporating the burn-in phase m.
REGRET GUARANTEES FOR TRUNCATED LEARNING
This section contains theoretical results, providing answers the research questions stated in Section 1. Our results indicate how the design parameters S (number of the subsequences). N (length of the subsequences), and burn-in phase length of the network influence the training performance of RNNs when pursuing truncated learning approach. Here, we assess the training in terms of the achieved performance gap (i.e., the regret) with respect challenging benchmark.
4.1 AN OPTIMIZATION PERSPECTIVE
The averaged gradient A02 in (8) can be interpreted as an unbiased estimate of the gradient of the full-batch loss function, that is, the right-hand side of (8) with B and b replaced by S and S, respectively, compare Tian et al. (2023). Hence, TBPTT as implemented by Algorithm actually aims solving the following optimization problem:
V min
1 E L(0, Di) min S(N m) 111=+1+1 S
1
Σ
Σ
Σ
10, - 1/ (10)
The optimization objective in (10) involves hidden state sequences of length N ++1, each initialized with zero, coupled by the same parameterization 0. This is conceptually different from the definition of the performance metric in (9), which assesses the training performance based on single hidden state sequence of length T + 1. This naturally raises the question: What if we leverage this additional knowledge during optimization? This motivates considering the following benchmark problem:
bb
min (be) S(N - m)
I1y (he 0, x 11
(11
(=1 j=m+1
subject to hr 1(11-1,21 0),
T.
(11b)
Note that the benchmark problem in (11b) uses the original TBPTT optimization objective from (10), i.c. aims minimizing the averaged MSE of S output sequences length N. However, the additional constraint in (11b) ensures that the underlying hidden state sequences actually correspond to different subsequences of single long hidden state sequence (h:) Comparing the solutions (10) and (11) therefore allows to investigate and quantify the effects of using zero initialization in TBPTT instead of the optimal initializations. Since both problems share the same optimization objective, we can perform this comparison in terms of: 1) the training regret, quantifying the optimality gap V (see Theorem 1); and 2) the performance regret, where we apply the parameterized networks to the full training imput sequence and evaluate the performance metric (9) (see Theorem 2).
In the following, we refer to the solution of the TBPTT problem (10) with the parameter 0*, and to the benchmark problem in (11) with the hidden state sequence ( T and network parameter @b.
4.2 MAIN RESULTS
In this section, we present novel accuracy and regret guarantees for truncated learning. To provide insights into the training process and its tuning knobs that are also valuable practitioners, we focus

---

Published as a conference paper at ICLR 2026
on the key messages and defer most technical details to the appendix; moreover, we formulate theorems rather informally in the main body of the paper, while they are made precise in Appendix A.
In line with standard learning techniques, that (a7 , 3/) ) E x²d x yd for all f = 1,., T, where the sets X' C. Rd. and yrd C. R are compact. Usual design choices include the intervals [0. 1] or [-1, 1] for the elements of each data point. We focus on parameterizations that result in the RNN having the following output stability property.
suppose that the training data is normalized in the sense
Assumption 1. There exists a (non-emptyy set O C TR such that, for all 0 E Q the network (1)-(2) is exponentially incremenallyy output stable along input sequences of the training data, that is, there exist constants C > and X E (0, 1) such that
l1 llychh (1) X) - yc( (2) (0,X)II S CX
hn 11 11. t- 1.2...... T"
(12)
for any two initial hidden states h10) (/2)) E " and any sequence X of length T" E [0,T] extracted
from the training input sequence X
Assumption essentially requires that the RNN "forgets" its initialization over time, so that for the same input sequence X, it always produces output sequences that converge each other, despite different initial hidden states. This is a standard condition (I.jung. 1978), and moreover, generally represents a desirable feature in the context of RNNs. Indeed, when applying an RNN in practice, should always generate meaningful outputs that are not overly sensitive to its initialization. We want to emphasize that Assumption 1 does not refer to a property of some underlying (unknown) data- generating system-ii it is rather property of the trained RNN model. Consequently, one could directly enforce it during training by normalizing and clipping certain RNN weights when performing the gradient updates such that the RNN forward dynamics essentially become stable (i.e., contractive), see Miller & Hardt (2019) for more details and corresponding methods (applicable to, e.g., LSTMs) and compare also Bonassi et al. (2021); Kolter & Manek (2019); Kozachkov et al. (2022). Nevertheless, invoking these sufficient conditions is generally conservative, potentially limiting expressiveness of the RNN. However, note also that contractive dynamics not necessary for Assumption 1 1-rather, it suffices that the trained model satisfies the property in (12) only along input sequences of the training data Xd, which can be easily tested numerically after training, compare Section 5.
Our first result bounds the gap between solutions to the TBPTT problem (10) and the benchmark (11). Theorem 1 (informal). Let the TBPTT solution 0* of the problem in (10) and the benchmark solution 9 of the problem in (11) be given and assume that both parameters render the respective RNN incrementally output stable in the sense of Assumption
Under an additional technical condition, there exist constants C2, C2 >> 0 such that the following properties hold:
1. Outpul accuracy:
Σ
11 (0, 0", x us (h,, 06, x,DiI S CiX"
j=m+1
xm
2. Training regret:
V
Vb S
C2 N m
Here, X E (0, 1) is from Assumption 1.
We provide the proof of Theorem in Appendix A.3 (Corollary 1). The first property in Theorem 1 states that the sum of the averaged squared errors between the predictions generated by the TBPTT solution 0 and the benchmark solution 0 along all subsequences S is bounded by a constant. As this constant does not depend on N, this sum must become smaller if N is increased, which necessarily implies that the TBPIT outputs get closer to the benchmark outputs. Moreover, the bound gets smaller for larger values of m (as the outputs are less influenced by transients) and for smaller values of X (as the outputs converge faster). The second property given in Theorem reveals that the optimal value for m when aiming for a small regret bound is in fact highly dependent on the convergence rate X; for slow convergence (larger values of A), smaller values of m are preferable, and vice versu.
Remark 2. The statement of Theorem I (and Theorem 2 below) involves "an additional technical condition", " As we show in detail in Appendix A, we need to exclude cases where the optimal output sequences fulfill a certain geometric relation, leading to the same value functions despite being different. However, this is rather of formal interest. since 1) it can be avoided in advance by choosing a suitable RNN architecture; 2) we conjecture that our theory can be generalized to these cases as well while preserving the key implications. For further details, see Remark 4 in Appendix A.2.
7

---

Published as a conference paper at ICLR 2026
Our second result assesses the trained RNN when processing the full training input sequence X Theorem 2 (informal), Let the TBPTT solution 0 of the problem in (10) and the benchmark solution 0 of the problem in (11) be given and assume that both parameters render the respective RNN incrementally output stable in the sense of Assumption 1. Let the burn-in phase parameter m satisfy m Omiu, where owin is the minimum overlap between subsequences as defined in (6). Under an additional technical condition, there exist some E1, Eg > such that the following properties hold:
T
(S - 1)A2 + S
Σ 10,, ,x) - (hr( 0b, )22 S E
1. Outpul accuracy:
T- 7TB 1=m+1
T m
(S - 1)) + SAm
2. Performance regret:
P(O, D) - P(hj. 0h; D) - Fg
--m
Here, X E (0, is from Assumption I and P corresponds to the training performance defined in (9).
Theorem 2 is proven in Appendix A.4 (Theorem 3). It reveals similar dependencies of the tuning parameters on the training performance as Theorem 1, ensuring that the output sequence generated by the TBPTT solution is close to that of the benchmark, also when processing the full input sequence. However, there are also subtle differences: First, there is an additional condition requiring m S Ominn- This ensures that the full-batch loss functions used in (10) and (11) involve all data samples [y Second, there an additional influence of the number of training subsequences S. Still, the regret bounds are mainly influenced by the parameters m and X, whereby smaller values of m are preferable for > and vice versa. If an upper bound for À is enforced during training, the optimal value of m minimizing the regret estimates from Theorems and 2 could, in principle, be computed explicitly. However, the estimates are generally conservative, as Assumption (on which they are based) also covers the slowest possible decay and is therefore inherently conservative in quantitative terms. Our theoretical results should therefore be understood as qualitative tuning guidelines.
EXPERIMENTS
We validate our theory on various time series tasks (the code is available online at https: //github. com/ httph/dschhheeeeeeeeeeaeeaeeeeee caale/rnntttte training). We first consider a small-scale problem of academic nature in Section 5.1, where we rely on synthetic data. This allows us to solve the bench- mark problem (11) using direct nonlinear optimization methods and hence to illustrate Theorems 1 and 2. In Section 5.2, we show that the obtained insights carry over to RNN training applied public data sets from the fields of system identification and time series forecasting, using classical TBPTT and mini-batch SGD methods. Here, we numerically tested the TBPTT-trained models for compliance with Assumption 1, which was always fulfilled with À = 0.99, ..... 0.9999 < 1. However, as the benchmark problem (1I) can no longer be directly solved due the size of the problem, we assess RNN performance by considering the MSE of the generated outputs with respect to training and test data and additionally compare it with BPTT and stateful (see Remark 1) TBPTT training results.
5.1 SYNTHETIC DATA
We aim to model an unknown linear single-input single-output system, for which noisy training data set of length T = 100 is available (see Figure 5 in Appendix B.1). To this end, we design a lincar RNN (Proca et al., 2025) with one-dimensional hidden state, corresponding to the model from (3)
with <(x) = x and by = by- = O, where we additionally enforce Whall ll < 0.999 during optimization to ensure Assumption 1. For training,
select N 21, set S T N + 1 80 (compare Section 3.1), and compute the solutions to the TBPTT problem (10) and benchmark problem (11) for different values of m using the nonlinear programming solver IPOPT (Wächter & Biegler, 2005) in MATLAB. The key training results depicted in Figure 3 are in line with our theoretical results provided in Section 4.2; specifically, the left plot shows that the j-step batch-averaged output error ej : > Ilv, (0, , x - y (hb; 1bb,,0 "b, ,X') converges to neighborhood of zero, the size of which decreases with increasing m. Hence, larger values of allow for longer transient phases, and the RNN effectively exploits this during training to get closer to the benchmark outputs. As indicated by the right plot of Figure 3, this behavior is also reflected in the resulting performance achieved over the training and an additional test data set (compare Appendix I. Here, the regret exponentially converges to zero with increasing m. and almost zero regret can be achieved when using m ' 12.
8

---

Published as a conference paper at ICLR 2026
10
8V8. output eITor
10
10
performance
10
test data
12
10-3
training data
10
m
16
10-9
m
m
20
10
X
6
11
16
21
8
16
time j
burn-in phase length
Figure 3: Training results for the synthetic data experiment for different values of the burn-in phase m. Left: batch-averaged output error ej. Right: performance P (the truncated MSE from (9)) achieved by solving the TBPTT problem (10) (*o' markers) and the benchmark problem (11) ("x' markers).
Silver-Box
Electricity
1
10
8823
8808
23838
21255
34444
- N 48
0.2
N -96
N 1192
10
10
OC
0.15
8 60
Xx
10-4
==100
- =5000
0.1
10a2
N 200
- -1000
Silver-Bi Box
W-H
RLC
Electricity
Traffic
Solar
10
102
103
50
100
150
200
burn-in phase
burn-in phase
Figure 4: Training results. Left and middle: MSE of RNN predictions for training ('o") and test data (*x") after (re- )training; gray lines correspond to the BPTT-trained RNN evaluated on the training (solid) and test data (dashed). Right: Averaged TBPTT training times relative to BPTT.
5.2 DATA SETS FROM SYSTEM IDENTIFICATION AND TIME SERIES FORECASTING
In this section, we consider the training of RNNs on public data sets from the fields of I) system identification: Silver-Box (Wigren & Schoukens, 2013), Wiener- Hammersennni (W-H) (Schoukens et al., 2009), RLC circuit (RLC) (Forgione & Piga, 2021a); and 2) time series forecasting: Electricity, Traffic, Solar EoleryEnergy (Lai et al., 2018; Wang et al., 2025).
For system identification,, we aim to design an RNN as a black-box model of the state-space represen- tation of the respective underlying dynamical system. Here, choose LSTM cell with hidden state dimension d₁ = 8 and a linear output layer, see Appendix B for further details. For each data set, we train multiple RNN instances using TBPTT and various combinations of window length N and burn-in phase m. In Figure 4 (left) and Figure 6 in Appendix B, we observe a strong influence of on the MSE achieved the training and test data. Interestingly, the results show very similar qualitative behavior when varying m, regardless of N. The training results are summarized in Table I (upper half), where we chose m² based on an additional validation data set and refer m baseline choice. Consequently, proper tuning of m can yield significant reduction of the MSE on both training and test data compared to baseline TBPTT, narrowing the gap to BPTT performance on the training data while avoiding the risk of overfitting and reducing the training time (see Figure 4).
the
For time series forecasting, we design an RNN to predict future values of the time series contained in the lookback window of length N. compare Figure I (right).. Here, choose F = 96 and consider lookback windows of different lengths N 48, 96, 192, compare Wang et al. (2025). To simplify the training, restrict ourselves to univariate forecasting and consider one LSTM cell with hidden state dimension d,, = 96 and a linear output layer. For each data set, we train multiple RNN instances for different values of N and m. To this end, split the available time series into training, validation,

---

Published as a conference paper at ICLR 2026
Table I:
Training on public system identification (sysid) and time series forecasting (TSF) data sets.
MSE training data
MSE test data
N
m
TBP m
TBP m
TBP sf
BP
TBP 17
TBP m
TBP sf
BP
W-H Silver-Box
100
94
0.242
0.042 (-83%)
0.059
0,003
0.291
0.142 (-51%)
0.222
0.251
200
94
0.050
0.030 (-41%) ))
0.048
0.003
0.140
0.135 (-4%)
0.221
0.251
500
94
0.030
0.019 (-37%)
0.051
0.003
0.121
0.122 (+1%)
0.185
0.251
1000
94
0.026
0.012 (-54%)
0.032
0.003
0.117
0.112 (-4%)
0.146
0.251
100
21
0.220
0.153 (-31%)
0.541
0.063
0.511
0.488 (-4%)
1.014
0.445
200
36
0.168
0.138 (-18%)
0.514
0.063
0.454
0.476 (45%)
0.961
0.445
500
0
0.135
0.135 (#0%)
0.465
0.063
0.397
0.397 (+0%)
0.875
0,445
1000
0
0.116
0.116 (#0%)
0.391
0.063
0.372
0.372 (+0%)
0.737
0.445
100
83
2,290
1.122 (-51%) %)
1.101
0.107
3.527
1.093 (-69%)
1.183
0.176
RLC
200
110
0.971
0.309 (-68%)
1.057
0.107
1.462
0.307 (-79%)
1.186
0.176
500
130
0.442
0.186 (-58%)
0.716
0.107
0.538
0.228 (-58%)
0.932
0.176
1000
110
0.306
0.176 (-42%)
0.674
0.107
0.402
0.362 (-10%)
0.871
0.176
Electr.
48
10
0.107
0.101 (-6%)
0.165
0.103
0.150
0.146 (-3%)
0.418
0.227
96
14
0.098
0.089 (-10%)
0.160
0.103
0.159
0.149 (-6%)
0.279
0.227
192
20
0.099
0.083 (-16%)
0.152
0.103
0.146
0.144 (-1%)
0.219
0.227
Traffic
48
0
0.297
0.287 (-3%)
0.321
0,240
0.343
0,316 (-8%)
0.728
0,498
96
75
0.228
0.234 (+3%)
0.331
0.240
0.305
0.280 (-8%)
0.813
0.498
192
20
0.251
0.205 (-18%)
0.370
0.240
0.284
0.276 (-3%)
0.359
0.498
Solar
48
31
0.300
0.312 (+4%)
0.312
0.152
0.250
0.257
(+3%)
0.772
0.110
96
75
0.068
0.069 (+2%)
0.203
0.152
0.227
0.184
(-19%)
0.349
0.110
192
141
0.020
0.020 (-1%)
0.185
0.152
0.211
0.112 (-47%)
0.092
0.110
RNNs are evaluated over the interval [1, 7) (sysid) lookback windows of length N (TSF); "TBP m":
baseline TBPTT training using m-m=0 (sysid) and m-m-N -- (TSF); "TBP m TBPTT training with m=m", where values in parentheses denote relative change with respect the respective bascline; "TBP sf stateful TBPTT training using m=n; "BP": full BPTT training using m=00. The MSE values corresponding system identification datasets (Silver-Box, W-H, RLC) scaled by a factor of 10².
and test segments, where we first (pre-)train the RNN on the training segment, tune m² by additionally considering the validation segment, re-train the RNN using the combined training and validation segments, and ultimately evaluate the forecast performance on the test segment, Appendix B.2 for further details. Figure 4 (middle) and Figure 6 in Appendix B.2 show the forecast MSE after re- ee-training the RNN using the considered combinations of N and m. Here, we recall that N corresponds to the standard configuration in time series forecasting (compare the discussion at the end of Section 3.2), which we consider the the baseline TBPTT with corresponding value m. In the bottom half of Table 1, we compare the MSE achieved after training using m = m with the baseline TBPTT, stateful TBPTT (i.c.. passing hidden states over lookback windows during training), and full BPTT using m = 0 (which corresponds to inferring the forecast F using all available past data). Again, we find that tuning m² is beneficial overall, we can achieve a significant reduction in the training and/or test MSE in the majority of cases. Choosing m < m. can be interpreted as additional regularization, which turns out to be numerically advantageous when performing the gradient updates during training. The fact that TBPTT with zero initialization outperforms stateful TBPTT and BPTT training in some cases in Table is due to the positive effect of randomness in the training process (the data batches are shuffled, see Remark 1). which renders the training more numerically stable and yields faster convergence (note that we use the same number of iterations for all training methods).
CONCLUSIONS
In this paper, we developed the concept of burn-in phase from a heuristic to a theoretically grounded method for RNN training. Our results indicate how the training hyperparameters in particular, the length of the burn-in phase- phase-influence influence the training performance. In experiments on public data sets, showed that tuning the burn-in phase improves RNN performance compared to standard choices.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work was supported by the Deutsche Forschungsgemeinschaatt (DFG, German Research Founda- tion) as part of the Research Unit "Active Learning for Systems and Control (ALeSCo)" project number 535860958.
REFERENCES
Naman Agarwal, Brian Bullins, Elad Hazan, Sham Kakade, and Karan Singh. Online control with adversarial disturbances. In Proc. Int. Conf. Mach. Learn. (ICML), volume 97 of PMLR, pp- 111-119, 2019.
Priyank Agrawal, Jinglin Chen, and Nan Jiang. Improved worst-case regret bounds randomized least-squares value iteration. Proc. AAAI Conf. Artif. Intell., 35(8):6566-6573 2021.
Christopher Aicher, Nicholas J. Foti, and Emily B. Fox. Adaptively truncating backpropagation through time to control gradient bias. In Proc. 35th Conf. Uncertain. Artif. Intell, volume 115 of PMLR. pp. 799-808, 2020.
Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. On the convergence rate of training recurrent neural networks. In Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 32, 2019.
Gerben I. Beintema, Roland Tóth, and Maarten Schoukens. Nonlinear state-space identification using deep encoder networks. In Proc. 3rd Conf. Learn. Dynamics Control (L4DC), volume 144 of PMLR, pp. 241-250, 2021.
Gerben I. Beintema, Maarten Schoukens, and Roland Tóth. Deep subspace encoders for nonlinear system identification. Automatica, 156:111210 2023.
Yoshua Bengio, Patrice Simard, and Paolo Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE Trans. Neural Netw., 5(2):157-166, 1994.
Filippo Maria Bianchi, Enrico Maiorino, Michael C. Kampffmeyer, Antonello Rizzi, and Robert Jenssen. Recurrent Neural Networks for Short-Term Load Forecasting: An Overview and Compar- ative Analysis. Springer, Cham, Switzerland, 2017.
Fabio Bonassi, Marcello Farina, and Riccardo Scattolini. On the stability properties of gated recurrent units neural networks. Syst. Control Lett., 157:105049, November 2021.
Fabio Bonassi, Marcello Farina, Jing Xie, and Riccardo Scattolini. On recurrent neural networks for learning-bast based control: Recent results and ideas for future developments. Process Contr, 1 114: 92-104, 2022.
Léon Bottou and Olivier Bousquet. The tradeoff's of large scale learning. In Adv. Neur Inf. Process. Syst. (NeurIPS), volume 20, 2007.
Jean-Séhastien Brouillon, Florian Dörfler, and Giancarlo Fernari-Trecate. Minimal regret state estimation of time-varying systems. IAAC-PapersOnLLiee, 56(2) 2023.
Vanessa Buhrmester, David Münch, and Michael Arens. Analysis of explainers of black box deep
neural networks for computer vision: a survey. Mach. Learn. Knowl. Extr., 3(4):3664))))))) 989, 2021.
Dean A. Carlson, Alain B. Haurie, and Arie Leizarowitz. Infinite Horizon Optimal Control. Springer, Berlin, Heidelherg, Germany, 1991.
Kyunghyun Cho, Bart van Merrienboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using RNN encoder-decoder for statistical machine translation. In Proc. Conf. Empir. Methods Nat. Lang. Process. (EMNLP), pp. 1724-1734, 2014.
Dario Coscia. Max Welling, Nicola Demo, and Gianluigi Rozza. BARNN: A bayesian autoregressive and recurrent neural network. In Proc. Int. Conf. Mach. Learn. (ICMI.), volume 42, 2025.
11

---

Published as a conference paper at ICLR 2026
Jeffrey L. Elman. Finding structure in time. Cognitive Sci., 14(2):179-21 1990.
Steven Elsworth and Stefan Güttel. Time series forecasting using LSTM networks: A symbolic approach. arXiv:2003.05672 05672,, 2020.
Timm Faulwasser and Lars Grüne. Tumpikk properties in optimal control. In Emmanuel Trélat and Enrique Zuazua (eds.), Handbook of Numerical Analysis, volume 23, pp. 367-400. Elsevier, Amsterdam, The Netherlands, 2022.
Timm Faulwasser, Lars Grüne, and Matthias A. Müller. Economic nonlinear model predictive control. Found. Trends Syst. Control, 5(1):1 2018
Timm Faulwasser, Lars Grüne, Jukka-Pekka Humaloja, and Manuel Schaller. The interval tumpike property for adjoints. Pure Appl. Funct. Anal.,
2022.
Marco Forgione and Dario Piga. Continuous-time system identification with neural networks: Model structures and fitting criteria. Eur. J. Control, 59:69-88 81, 2021a.
Marco Forgione and Dario Piga. dynoNet: A neural network architecture for learning dynamical systems. Int. J. Adapt. Control Signal Process, 35(4))612-62) 2021b.
Marco Forgione, Manas Mejari, and Dario Piga. Learning neural state-space models: Do we need
state estimator? arXiv:2206. 12928, 2022.
Karol Gregor, Ivo Danihelka, Alex Graves, Danilo Rezende, and Daan Wierstra. DRAW: A recurrent neural network for image generation. In Proc. Int. Conf. Mach. Learn. (ICML), volume 37 of PMLR, pp. 1462-1471, 2015.
Lars Grüne and Matthias A. Müller. On the relation between strict dissipativity and turpike properties. Syst. Control Lett., 90:45-55 53, 2016.
Lars Grüne and Simon Pirkelmann. Economic model predictive control for time-varying system: Performance and stability results. Optim. Control Appl. Methods, 41(1):42-64 2019.
Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces. 12200:23/2.00752. 2024.
Albert Gu, Karan Goel, and Christopher Re. Efficiently modeling long sequences with structured state spaces. In Proc. Int. Conf. Learn. Represent. (ICIR), 2022.
Hansika Hewamalage, Christoph Bergmeir, and Kasun Bandara. Recurrent neural networks for time series forecasting: Current status and future directions. Int. J. Forecasting, 37(1):388 -2443888877,, 2021.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Comput., 173: 1780, 1997.
Herbert Jaeger. Tutorial on training recurrent neural networks, covering BPTT, RTRL, EKF and the echo state network approach. German National Research Center for Information Technology, GMD Report 159, 2002.
Thomas Jaksch, Ronald Ortner. and Peter Auer. Near-optimal regret bounds for reinforcement learning. J. Mach. Learn. Res., (1(51):1563-166 2010.
Alexander Katrompas and Vangelis Metsis. Enhancing LSTM models with self-attention and stateful training. In Proc. SAI Intell. Syst. Conf. pp. 217-235, 2021.
Wahab Khan, Ali Daud, Khairullah Khan, Shakoor Muhammad, and Rafiul Haq. Exploring the frontiers of deep learning and natural language processing: A comprehensive overview of key challenges and emerging trends. J. Nat. Lang. Process., 4:100026, 2023.
Jongseon Kim. Hyungjoon Kim, HyunGi Kim, Dongjun Lee, and Sungroh Yoon. A comprehensive survey deep learning for time series forecasting: architectural diversity and open challenges. Artif. Intell. Rev., 58:216, 2025.
12

---

Published as a conference paper at ICLR 2026
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Proc. Int. Conf. Learn. Represent. (ICLR), volume 3, 2015.
J. Zico Kolter and Gaurav Manek. Learning stable deep dynamics models. In Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 32, 2019.
Leo Kozachkov, Michaela Ennis, and Jean-Jacques Slotine. RNNs of RNNs: Recursive construction of stable assemblies of recurrent neural networks. In Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 35, 2022.
Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu. Modeling long- and short-term temporal patterns with deep neural networks. In Proc. 41st Int. ACM SIGIR Conf. Res. Dev. Inf. Retr, pp. 95-104. ACM, 2018.
Yingying Li, Xin Chen, and Na Li. Online optimal control with linear dynamics and predictions: Algorithms and regret analysis. In Adv. Neural Inf. Process. Syst. (NeurIPS), volume 32, 2019.
Zhaoyang Li, Minghao Han, and Xunyuan Yin. MamKO: Mamba-based Koopman operator for modeling and predictive control. In Proc. Int. Conf. Learn. Represent. (ICLR), volume 13, 2025.
Zachary C. Lipton, John Berkowitz, and Charles Elkan. A critical review of recurrent neural networks for sequence learning. arXiv:/506.000/// 2015.
I.. Ljuug. Convergence analysis of parametric identification methods. IEEE Trans. Autom. Control, 23(5):770-783, 1978.
Lennart Ljuug. System Identification: Theory for the User. Prentice Hall, Upper Saddle River, NJ, USA, edition, 1999. 14th printing.
Andrea Martin, Luca Furieri, Florian Dörfler, John Lygeros, and Giancarlo Ferrari- Treeatt. On the guarantees of minimizing regret receding horizon. IEEE Trans. Autom. Control, 70(3): 1547-1562, 2025.
Daniele Masti and Alberto Bemporad. Learning nonlinear state-space models using autoencoders. Automatica, 129:109666, 2021.
Tomas Mikolov, Stefan Kombrink, Lukas Burget, Jan Cernocky, and Sanjeev Khudanpur. Extensions of recurrent neural network language model. In Proc. IEEE Int. Conf. Acoust. Speech Signal Process., pp. 5528-5531, 2011.
Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. Distributed representations of words and phrases and their compositionality. In Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 26, 2013.
John Miller and Moritz Hardt. Stable recurrent models. In Proc. Int. Conf. Learn. Represent. (ICLR), 2019.
Nima Mohajerin and Steven L. Waslander, State initialization for recurrent neural network modeling of time-series data. In Proc. Int. Jt. Conf. Neural Netw., pp. 3330-2337, 2017.
Marko Nonhoff and Matthias A. Müller. Online convex optimization for data-driven control of dynamical systems. IEEE Open Control Syst., 1:180-193, 2022.
Antonio Orvieto, Samuel L. Smith, Albert Gu, Anushan Fernando, Caglar Gulcehre, Razvan Pascanu, and Soham De. Resurrecting recurrent neural networks for long sequences. In Proc. Int. Conf. Mach. Learn. (ICML), volume 202 of PMLR, pp. 186670226695 2023.
Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In Proc. Int. Conf. Mach. Learn. (ICML), volume 28 of PMLR, pp. 1310-1318, 2013.
Vivak Patel, Shushu Zhang, and Bowen Tian. Global convergence and stability of stochastic gradient descent. In Adv. Neur Inf. Process. Syst. (NeurIPS), volume 36, pp. 36001-36622, 2022.
13

---

Published as a conference paper at ICLR 2026
Gianluigi Pillonetto, Aleksandr Aravkin, Daniel Gedon, Lennart Ljung, Antônio H. Ribeiro, and Thomas B. Schön. Deep networks for system identification: A survey. Automatica, 171:111907,
2025.
Alexandra Maria Proca, Clémentine Carla Juliette Dominé, Murray Shanahan, and Pedro A. M. Mediano. Learning dynamics in linear recurrent neural networks. In Proc. Int. Conf. Mach. Learn. (ICML), volume 42, 2025.
James B. Rawlings, David Q. Mayne, and Moritz Dichl. Model Predictive Control: Theory, Com- putation, and Design. Nob Hill Publish., LLC, Santa Barbara, CA, USA, edition, 2020. 3rd printing.
António H. Ribeiro, Koen Tiels, Jack Umenberger, Thomas B. Schön, and Luis A. Aguirre. On the smoothness of nonlinear system identification. Automatica, 121:109158, 2020.
Julian D. Schiller, Lars Grüne, and Matthias A. Müller. Performance guarantees for optimization- based state estimation using turmpike properties. arXiv:2501. /8385, 2025. Provisionally accepted for publication as full paper in IEEE Trans. Autom. Control.
Joannes Schoukens, Johan Suykens, and Lennart Ljung. Wiener-Hammersienr benchmark. 15th IFAC Symp. Syst. Ident. (SYSID), 2009.
Max Schussler, Tobias Munker, and Oliver Nelles, Deep recurrent neural networks for nonlinear system identification. In Proc. IEEE Symp. Ser. Comput. Intell., pp. 448-454, 2019.
Ilya Sutskever, James Martens, and Geoffrey Hinton. Generating text with recurrent neural networks. In Proc. Int. Conf. Mach. Learn. (ICML), pp. 1017-1024, 2011.
Hao Tang and James Glass. On training recurrent networks with truncated backpropagation through time in speech recognition. In Proc. IEEE Spok Lang. Technol. Workshop, pp. 48 48-55, 2018.
Yingjie Tian, Yugi Zhang, and Haibin Zhang. Recent advances in stochastic gradient descent in deep learning. Mathematics, 11(3):682, 2023.
Emmanuel Trélat and Can Zhang. Integral and measure- turmpike properties for infinite-dimensiona optimal control systems. Math. Control Signals Syst., 30(1), 2018
Emmanuel Trélat and Enrique Zuazua. Turniike in optimal control and beyond: survey. arXiv:2503.20342, 2025.
Jeen-Shing Wang and Yen-Ping Chen. A fully automated recurrent neural network for unknown dynamic system identification and control. IEEE Trans. Circuits Syst. 1, 1366)))33331137,, 2006.
Zihan Wang, Fanheng Kong, Shi Feng, Ming Wang, Xiaocui Yang, Han Zhao, Daling Wang, and Yifei Zhang. Is Mamba effective for time series forecasting? Neurocomput, 619:129178, 2025.
Paul J. Werbos. Backpropagation through time: what it does and how to do it. Proc. IEEE, 78(10): 1550-1560, 1990.
Torbjorn Wigren and Johan Schoukens. Three free data sets for development and benchmarking in nonlinear system identification. In Proc. Eur. Control Conf. (ECC), pp. 2933-2938, 2013.
Ronald J. Williams and Jing Peng. An efficient gradient-based algorithm for on-line training of recurrent network trajectories. Neural Comput., 2(4):490-501 1990.
Andreas Wächter and Lorenz T. Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program., 106(1):25-57) 2005.
Shaoming Xu, Ankush Khandelwal, Xiang Li, Xiaowei Jia, Licheng Liu, Jared Willard, Rahul Ghosh, Kelly Cutler, Michael Steinbach, Christopher Duffy, John Nieber, and Vipin Kurnr.. Mini-batch learning strategies for modeling long term temporal dependencies: A study in environmental applications. In Proc. SIAM Int. Conf. Data Mining, pp. 649-657, 2023.
Alexander J. Zaslavski. Turnpike Properties in the Calculus of Variations and Optimal Control. Springer, New York, NY, USA, 2006.
14

---

Published as a conference paper at ICLR 2026
Mingxuan Zhang, Yan Sun, and Faming Liang. Sparse deep learning for time series data: Theory and applications. In Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 37, 2023.
Nicolas Zucchet and Antonio Orvieto. Recurrent neural networks: vanishing and exploding gradients are not the end of the story. In Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 38, 2024.
Nicolas Zucchet, Robert Meier, Simon Schug, Asier Mujika, and Joao Sacramento. Online learning of long-range dependencies. Adv. Neur. Inf. Process. Syst. (NeurIPS), volume 37, 2023.
A THEORETICAL ANALYSIS
This section contains proofs and further technical details of Theorems and 2 stated in Section 4.2. Our analysis builds on a generalized optimization problem (Section A.1), which captures the TBPTT problem from (10) and the benchmark problem from (11). In Section A.2, show that this generalized problem exhibits an output turppike property, which is used in Sections A.3 and A.4 to prove Theorems I and 2, respectively.
A.1 THE LIFTED OPTIMIZATION PROBLEM
To allow for a precise technical analysis, re-formulate the problems (10) and (11) in multiple shooting fashion (see, e.g., Section 8.5 of the book by Rawlings et al. (2020) for further details on this topic). More specifically, implement the output mapping used in the respective loss functions as part of the problem constraints and introduce additional decision variables.
To this end, we define the containers of the predicted hidden state and output sequences as
H (hyje hijeNo.iies. Y ==
(13)
where No == U (O). Furthermore, we define the stage cost Z/p1 IRd - Rgo associated with stage N and subsequence S according
()((()) = lly "
(14)
Now, we can express the loss functions used in (10) and (11) as
N
Us S(N m)
(n(()))))
(15)
==1 j=m+1
and formulate the following optimization problem:
min Js(Y: D) HYY.
(16a)
subject to
M1: (16:0) 0), Vje E, Vies,
(16b)
8z1z g(hb)s 2/11, Vj N, Vies,
(16c)
(hac) E N+!, fune VV,, ViES
(16d)
E θ R(H) 0.
(16e)
(16f)
The constraints in (16b) and (16c) imply that over each subsequence E S, the corresponding hidden state and output sequences satisfy the network equations from (1) and (2) for all j E N. Condition (16d) ensures that these sequences have elements from the sets H c Idd and ) 18 respectively, where we assume that H. and y are compact. Note that this is required for technical reasons in the proof of our results below and in particular does not need to be incorporated in, e.g., Algorithm 1 (as suitable sets H. and y can be found in any case after successful training, compare also Remark 3 below). The constraint in (16e) ensures that the resulting parameter 0 leads network model that is incrementally output stable in the sense of Assumption 1.
Depending on the last constraint in (16f), consider three different settings:
15

---

Published as a conference paper at ICLR 2026
1. The TBPTT case:
R(H) [Ao]] hojss
(17)
This results in each initial hidden state hois E S, being initialized to zero. From theoretical perspective. the problem (16) under (17) fully equivalent the original TBPTT problem in (10) (provided that the solution to the latter satisfies Assumption I, i... H" E (). We denote the solution to (16) under (17) as (H", Y Y*,
2. The coupled case:
R(H)= [hoj2
hN- 1 ..., hays hy -oo-osss3 [S-1]
(18)
This condition implies that for each i E S with i 2 2, the corresponding initial hidden state hojs is equal to the hidden state -o where o; corresponds to the overlap between the subsequences and i as defined in (5). As consequence, all hidden state sequences contained in H are subsequence of the same global hidden state sequence of length T + 1, where the initial hidden state state hoj1 is left as a free optimization variable. We denote the solution to (16) under (18) as (Ho, Y 0%). Again, from theoretical point of view, the problem formulation in (16) under (18) is fully equivalent to the benchmark problem in (11) (provided that 0 E +). In particular, their solutions sets coincide in the sense that 0b 0e and
hb - E Vi
(19)
3. The unconstrained case:
R(H) =0 0.
(20)
Here, we do not enforce any additional constraints on the hidden state sequences contained in H (except compliance with the model equations (16b) and (16c)). We denote the solution to (16) under (20) as (Hxyyx, )²))
Remark 3. Without loss of generality, we assume that the compact set y used in (16d) is large enough such that yj 00 .x() y and y, 0s, X,i) E y for all j N and all i E S, and furthermore, ye(0, 0r, Xd)= E Y and ve(hi uo(ho 09 xa) E for all times E [1,T]
A.2 THE OUTPUT TURNPIKE PROPERTY
In this section, establish relation between (H*,"* Y O*) and the unconstrained solution (H y ,=,,0)))) 0²)). Note that this is an intermediate step which enables us to derive the desired relation between (H** Y*, 0*) and the coupled solution (HS, Ye, 0²) in Sections A.3 and A.4. Our analysis is inspired by classical turnpike theory, which is an important concept in the context of optimal control (Trélat & Zuazua, 2025; Faulwasser & Grüne, 2022; Zaslavski, 2006; Carlson et al., 1991). In particular, we prove a turnpike property of the optimization problem (16), characterized in the output space. For the sake of clarity, we provide the following roadmap:
1. We first establish an upper bound on the difference of the value functions resulting from solving (16) under (17) and (20) (see Lemma 1). When interpreting the parameter 0 as the (constant) control input to the dynamical system (I), this result can be interpreted as a stabilizabilltty property of the latter, characterized in terms of the cost function (15): When increasing the prediction horizon N, the bound on the cost difference stays the same, which implies that the optimal input 8* effectively stabilizes the "output reference" y
2. In Lemma we establish structural relation between the outputs Y* and y using the fact that both correspond to optimal solutions of the same optimization problem (16).
3. This result is used in Lemma 3 to prove coercivity property of the difference of the value functions resulting from solving ( (16) under (17) and (20): The difference is positive definite in the sum of the differences of the corresponding output predictions.
4. Finally, by employing Lemmas and 3, we can characterize the desired turppike property of the problem (16) (see Proposition 1), which states that Y must be close to y (i.c., the output tumpike) for most of the time. We provide a detailed interpretation of this result below Proposition 1.
16

---

Published as a conference paper at ICLR 2026
Lemma (Stabilizability). Let Assumption hold. Then, there exists constant C such that
ΣΣ j=1j=+1
Σ
- -1y S ČSA";
(21)
where À E (0, 1) is from Assumption 1.
Proof. The claim follows from the fact that 8 θ by (16e) and Assumption In particular, let us consider the candidate solution (H, Y, ) such that Øzソ: = y, (0, 000 xx.) for all j W and E S. Here, y contains S output sequences of length N, each of which is generated by the RNN with 6x initialized at zero, and being driven by the respective subsequences of the training input data x; assume that the compact set y is such that
i E S. Note also that, without loss of generality, Gile E for all j E N and all E S, compare Remark 3. Consequently, (H, Y, 8%)) satisfies the constraints (16b)-(16f) and (17) and thus constitutes a feasible candidate solution to the problem (16) under (17). By optimality of y we have that (Y(""; D) S JY(Y; D), which yields
Σ
(n(0)p) ()(())))
(22)
Because both 0zjc and Vije are elements of the compact set y for each j E N and i E S and the stage cost (g)(()) is quadratic in y, there exists Lipschitz constant Li > 0 such that
61000, .) (1(()))) & (0))))))) - (1))0) S L 0g 9/
(23)
uniformly for all Úyte Dif E y and all j N and i E S. (In fact, I, = 2(max [lvll max [IBll) is a valid choice, recalling that y and yd are compact). Moreover, due to the fact that the output predictions Øziz and Vif are generated using the RNN with the same parameter g E θ but different initial hidden states, can apply the property in (12) from Assumption for each d using (1) = and h (2) = Gfje which lets us infer that
1 o ll s Ca S ChN j JEN, IES,
(24)
where ì == ma |h] using (1) and the fact that H is compact. The combination of (22)-(24) together with the application of the geometric series then leads to
N
Σ
Σ
-=
Σ x= L,C Σ nm d N-++1
A
< C,
Σ
1-X T,Ch Σ
Am
1-X
-1j=-11
By defining := TCchA((1 - A) we obtain (21), which finishes this proof.
Lemma 2. Consider the solutions (H* Y 0²) and (Ho0 yo y00,000) There exists a constant e 2 such that
N
5
N
Σ
Σ +=
2 - 11 (U5 - 2
llu - 17.
(25)
1=1 j=m+1
Proof. The property in (25) follows from optimality of (H1*, Y*, 0*) and (/=,,>,00 In partic- ular, since (H* Y*, 0*) represents a feasible candidate solution the unconstrained problem (16) under (20), we know that Jy (0) < JY(Y"), which leads to
Σ Σ
(1((0gH)) - 4, x 2 0.
(26)
17

---

Published as a conference paper at ICLR 2026
Due to the fact that (yi((y) is quadratic in y E V, it is equivalent to its second-order Taylor polynomial. Consequently, we can express (11)(15))) in terms of [o)((0 as
(H((5))) = ((((()())) + Alll
dy (ujfe)
(ujj= wうこ + (0] vうリe) T -
= = + 2(0 v)(e) "Wfie - )タ + l1bji 0
where is the Hessian of lyt with respect to y E y and simply reduces 5 21p- Thus, the left-hand side of (26) can be re-written as
8
N
N
Σ
Σ
(115) u + 2(15ii =).a (v]i - vz())
1=1 j=m++1
1=1
j=m+1
(27)
The inequality in (26) then implies that
s
N
N
Σ
E
2(0
Σ
Σ
sfill,
1=
(=1 1=1++
which establishes the desired property in (25) for e and hence concludes this proof.
Lemma 3 (Coercivity). Let the property in (25) hold for some e >
Then, the solutions (H*, Y 0*)
and (H y 800) satisfy
N
N 9/
Σ
Σ
b((05p) - ()((())))))
Σ
(28)
--1=++1
=1j=++1
Proof. The claim follows by applying the property (25) to (27) in the proof of Lemma 2.
Remark 4. Requiring <> 1 in Lemma 3 is rather mild technical condition on Y with respect to Yo. In particular, using K 1/c E (0,1), the result fails only if. for Ei= Σn Ej=m+1 llu//
Siall 122 + 0, the following holds:
N Ej=m+1 2(U 1/ (V/ v9t) < n VK (0,1),
-15
ES E j=m+1 1155 - 0fjil?
which is the case if and only if
9
N
5
N
ΣΣ
VI 1z()) =
Σ j-1++
Σ
I|55 - 35 + O.
(29)
jij=m+1
i=1
Solutions to this equation correspond to the case where Js(Y";; D) Js= D) but - + 0 for some j E N with j 2 m + and i E S. This may occur y (specifically, its restriction to the interval [m + 1, NN) is non-unique; in particular, the above implies that Y then also corre- sponds to an optimal solution to the unconstrained problem. But this means that we can simply set (H , y²,,00===----- (H*, Y*, O'), rendering condition (25) again valid for any e >> 1. Alternatively, (29) could be avoided in ad suvvane by proper design of the RNN, ensuring that the outpul turnpike Y unique (which does not require uniqueness of the solution with respect the hidden state trajectory H or 800 compare also the discussion below Proposition 1). Moreover, we also conjecture that our analysis can be extended to explicitly cover the case where Y is non-unique using arguments inspired by Trélat & Zhang (2018); namely, by characterizing the turnpike property using turnpike set T containing all possible output sequences achieving the minimum cost (instead of one particular output sequence) and by deriving the results in terms of the point-to-set distance of the solution Y* to the turnpike set T. Finally, we would like to point out that in our synthetic experiment conducted in Section 5.1, where were able explicitly check this condition, it was always possible choose e >> 1 such that property (25) is satisfied (even when changing the network architecture and using different nonlinear activation functions).
18

---

Published as a conference paper at ICLR 2026
Proposition (Output turpiike property). Let Assumption 1 be satisfied. Consider the solutions (H", YY, , and (H Y _0)) and suppose that the property in (25) holds for some > Then,
there exists constant K >> 0 such that
N
8
Σ
Votall KAT,
(30)
j=++1
where X E (0, 1) is from Assumption
Proof. The statement follows by considering the difference of the loss functions
Σ
Σ
G1((5(1)) - - (1)()
1=1 f=m=1
and applying the lower bound from Lemma and the upper bound from Lemma 1. The definition of K := Ce/(e 1) with C from Lemma then leads to (30), which finishes this proof.
Some remarks are in order.
Proposition provides characterization of the j-step prediction error, averaged over all subsequences:
10 DR JEN.
Specifically, property (30) implies that the of efe over [m + 1, N] is uniformly bounded independent of the length of the subsequences N. For a fixed value of m. this must imply that ej must become smaller for larger values of N for most of the times j [m + NJ. By non-negativity of the norm, we can conclude similar statement for the output predictions corresponding to each subsequence E S: Vj i gets closer to Vji; for most of the times j [m + 1, .] when increasing the value of N. Hence, the latter can be regarded as the turnpike for the former, and the larger we choose the length N, the more elements of YyІ= are close V/fK'
The fact that the right-hand side of (30) is independent of N is in line with intuition: Differ- ences in the output predictions (and thus in the loss functions of the optimization problems under consideration) are caused solely by the burn-in phase of the RNN resulting from zero initialization, whose influence is limited and, for fixed network parameterization 0, does not depend on N.
The main difference to the tumpike-rrlated literature is that we formulate the turnpike property with respect to the outputs-without as usual, characterizing the behavior of optimal states, controls, and possibly adjoints, compare Trélat & Zuazua (2025); Faulwasser & Grüne (2022); Zaslavski (2006); Carlson et al. (1991). This is natural in our setting, as the set of hidden states and network parameters (i.e., the "controls"") that generates specific output sequence typically not a singleton (the same output sequence can usually be generated by variety of hidden states and parameters without significantly restricting the network architecture).
The turnpike characterization in (30) is conceptually related to the interval turnpike property from Faulwasser et al. (2022). Note that (30) can be easily reformulated in terms of measure/cardinality turnpike property by estimating the number of elements of Y that lie outside a certain e-neighborhood of the tumpike y Similar characterizations are promi- nent in the most recent mempike-rrllttd literature (Trélat & Zuazua, 2025; Faulwasser & Grüne, 2022; Grüne & Pirkelmann, 2019; Faulwasser al., 2018); however, the quantitative nature of the bound in (30) turns out to be more useful when establishing the results in Sections A.3 and A.4 below.
The fact that (Lemma 1) and coercivity (Lemma 3) is in line with standard turnpike results (compare, e.g., Trélat & Zuazua (2025); Trélat & Zhang (2018); Grüne & Müller (2016)). There,
establish the turnpike property (Proposition 1) using stabilizability
19

---

Published as a conference paper at ICLR 2026
however, coercivity is usually first established as consequence of additional dissipativity assumption on the problem, mainly due to the fact that the results are tailored to general cost functions. In contrast, here can leverage the quadratic nature of the MSE loss and thus avoid an additional dissipativity condition. Nevertheless, we would like to point out that, one might expect, Lemma 3 directly implies dissipativity property, where the storage function can be chosen as a constant and the supply rate involves a strictness term capturing the output differences.
A.3 ACCURACY AND REGRET OVER SUBSEQUENCES OF THE TRAINING DATA
In this section, we leverage the turnpike property provided by Proposition to establish accuracy and regret guarantees for the solution (H*, Y*,"*) with respect to the benchmark (H" Y 8²) The following should be noted in this regard.
Remark 5. The previous results (in particular, Lemmas 1-3 and Proposition 1) also hold true if (H*, Y ) in the statement of the results is replaced by the coupled solution (Ho, Y 7.) The respective proofs can be straighfforwardly adapted using the fact that (Ho y2,0%) is a solution to the optimization problem (16).
Hence, by combining Proposition and Remark 5, we can characterize the relation between (H*, Y+, (*,0)) 8*) and (H",Y "Y,,8" as shown in the following corollary.
Corollary 1. Let Assumption I be satisfied. Suppose that there exists some > 1 small enough such that the property in (25) holds true for both (H", YO,, following statements hold:
"x,00) ,0)) and (H y,,0),, see Remark 5. Then, the
1. Output accuracy: The learned output sequences contained in Y and Y satisfy
ΣS
Σ S
Σ 1151 - 4KX",
(31)
j=m+1
where K 0 is from Proposition and X (0, 1) is from Assumption 1.
2. Training regret: Y* and Y satisfy
A
Js(Y" - Js("; SC N
(32)
where C 0 is from Lemma and X e 1) is from Assumption 1.
Proof. The first statement is consequence of Proposition Remark and application of the Cauchy-Schwartz and Young's inequality. In particular, (31) follows by noting that
Σ
1b/|=
j=m+1
S
2
|1 v +2
Σ j=m-1
Σ
1b
5
$
j=++1
1=1
< 4m.
The second statement can be established by applying similar steps as in the proof of Lemma 1. In particular, using the candidate solution Kزرs M (0, 0= X(²) E V. E N, i E S, have that JY D) - (Y'; D) < JÝ: D) - JYs; ; D) by optimality of Y*. Now, we can apply the same arguments that followed (22), which leads to (32) and thus concludes this proof.
O
Remark 6. The estimate the training regret in (32) in fact also holds true for e 2 1, as it does not require the coercivity property established in Lemma 3.
Recalling the equivalence between the coupled solution (H Y 88) and the actual benchmark consisting of the hidden state sequence {hh}T and the parameter 6² as specified in (19), one can directly see that Corollary formalizes (and thus proves) Theorem 1.
20

---

Published as a conference paper at ICLR 2026
A.4 ACCURACY AND REGRET OVER THE FULL TRAINING DATA
In this section, we establish accuracy and regret estimates for the RNN parameterized by 0* when processing the full training input sequence X The estimates are given with respect the accu- racy achieved by the benchmark (hg. 8%) (that is, the the solution to (11)), where we measure the performance using the truncated MSE loss as defined in (9).
The following result constitutes the formal version of Theorem 2.
Theorem 3. Let Assumption 1 be satisfied. Suppose that there exists some > 1 small enough such that the property in (25) holds true for both (H* Y",0" and (Ho Y )) ) cf. Remark 5. Then, the following statements hold:
1. Averaged oulput accuracy: There exists constant Ej > such that
11b(0, 0", x²) ve 0b, x2) S Fg (S - 1)A20mm + SA
(33)
T t=++1
m
2. Performance regret: There exists a constant Fgg > O such that
P(O,0;;)) - P(hi, 0"; D) E2 Es v 1)A20min + SAT
(34)
m
Proof. We start with some useful definitions and observations. Let {hī ) be the hidden state sequence that generated by the RNN from (1)-2) parameterized by () when being evaluated over the full training input sequence X and initialized with bツ = 0. Define the corresponding output sequence as yr == (0,00* xx), t [1,7]. Similarly, define the benchmark output sequence g6 = (h 0h, ).), J. [1, Recall that (h) 0h) = (h 0*) by (19), which yields
Vー,+j-1 Vff EJEN, V.ES.
(35)
Without loss of generality, Remark 3.
assume that h; H and vi vi E y for all [1,7], compare
We proceed to prove the first claim of this theorem. To this end, note that re-write the of output differences (i.e., the left-hand side of (33)) in terms of the subsequences i E S using the definition of the overlap o from (5) according to
T
N
Σ (=1++-1
1l 1111=
Σ I jmm+1
lB - 12
IIV - 111+
(36)
Regarding the first sum on the right-hand side of (36), each vī coincides with the output of the first training subsequence in the sense that y = y, (0, 00 xd) = (0, (0,0", x = Vf11 for all j E N. Using the Cauchy-Schin Schwarz and Young's 's inequality together with (35), follows that
|| v 1| 1Ive+ 2ز 25 2 =
< 2||y - 2/ 2 + 211v2 - V3
for all j N and E S. Consequently, (36) can be bounded
N
S
Σ
uP11 <
Σ
Σ
11
t=
j=
1=2 j=o,+1
+
2
| 11
(37)
i=2
We first focus on the last sum on the right-hand side in (37). Here, each summand corresponds to the difference of network outputs resulting from different initial hidden states, but using the same
21

---

Published as a conference paper at ICLR 2026
parameter 0* which satisfies Assumption 1 by (16e). In particular, for each i E 8 with 2 2. from (12) we obtain
11 0ūjel1 llus.-j-1 x²) 1,(0, 0r xx!]11
11u; (h 0- xx)) (0, 0" )?
for all j E NV, where we have used the definition of ei == Ca maxMeN |||| (noting that ci is well-defined due to compactness of H and hī E H for all E T', see above). By summing over j and applying the geometric series, it follows that
8
N
N Σ
X²
Σ
l1, - 33 ==0++1
Σ
C1A2S
= Σ
(1- 22()-oa))
x2
i=2
1=0++1
i=2
S C2 (S (8-1)220mmn 20man
(38)
where 02 == :
The first line of the right-hand side in (37) can be bounded using the fact that the burn-in phase length m satisfies m S Omin I o, for all i 2 2 together with the turpikk property provided by Corollary 1, which leads to
5
Σ #=
1b5n W 1 +2
Σ
Σ
2 >
Σ
|R 3/
1=0--1
1=1 j=m+1
S 2 4KKS
(39)
The combination (39) then yields
T
Σ
y111 < (0 - 1)A 20mia 8SKA". llvi
t=++1
Using the definition E₁ : 2- max(c2, 4 and dividing both sides by (T (7-m) > 0, we obtain (33).
To prove the second statement, note that the performance criterion in (9) is Lipschitz on compact sets (recall out observations in the beginning of this proof, the arguments that lead to (23), and the definition of the stage cost (14)); consequently, the left-hand side of (34) can be bounded according to
P"O,0"; D) - P(hb. Db; D) = 0P(0,00; - P(h 0³; :D]l
T
1
Σ
v211 - -1112
121
T t=m+11
5 T
Σ
llat
(40)
t=m+1
Using concavity of the square root and Jensen's 's inequality, it holds that
T
T
T
1/2
Σ ll - w
Σ (=
Vili - 1/1 1
T m
Σ
1l - 3/112
(41)
VT \==+1
(=m+1
From (40), (41), and (33) proven above, we obtain
P(0,0": D) P(hi. D) < m VT mm
L T m
(E) ((S 1)A 20min + ()²) 1/2
(S 20mi + SA
S L
T
This establishes (34) and hence finishes this proof.
22

---

Published as a conference paper at ICLR 2026
-1
20
40
60
80
100
20
40
60
80
100
20
40
60
80
100
1
20
40
60
80
100
3
40
100
40
80
100
time t
time t
Figure 5: Supplementary results for the synthetic data experiment performed in Section 5.1. The left column refers the training data, the right column to the test data. The corresponding input-output time series are shown in the top row. In the middle row, we compare the outputs of the RNN trained with m (without burn-in)
terms of the outputs yJj; (belonging to subsequences of the training data X; i S) and the outputs y/ (generated by processing the entire input sequence X with the benchmark output 3d and the ground truth data y/" The bottom row shows the respective results when training the RNNs using the burn-in phase length N
20.
B ADDITIONAL MATERIAL FOR THE EXPERIMENTS
B. SYNTHETIC DATA
In this section, provide additional illustrative figures and interpretations of the results of the synthetic experiment performed Section 5.1.
Specifically, the left column in Figure compares the generated outputs of the trained RNN, the benchmark RNN, and the ground truth data. Without using a burn-in phase (i.e., m 0, middle row), the transient resulting from the zero initialization causes the output 3/ ye (e xd)) to be relatively far away from the benchmark output y/ Yr (hb, 0b, x (which employs optimal initialization, compare Section 4.1) and the ground truth data. This is because the optimizer tries to minimize the influence of this transient by minimizing the convergence rate, which comes at the expense of long- oon-tterm accuracy. In contrast, when using burn-in phase m > 12, this transient is effectively excluded from the cost function, letting the outputs Vÿバ yı (0, 0*, X converge closely to the benchmark (and the ground truth data).
According to the right column Figure 5, these observations directly carry over to performance of the trained RNN on test data. Consequently, tuning the burn-in phase during training enables the RNN to achieve benchmark performance on both the training and the test data.
B.2 DATA SETS FROM SYSTEM IDENTIFICATION AND TIME SERIES FORECASTING
This section contains supplementary material on the experiments conducted in Section 5.2. Parameters depending the respective data sets reported in Table 2.
23

---

Published as a conference paper at ICLR 2026
Table 2: Training parameters for the experiments performed in Section 5.2.
data set
T
Tty
<x
epochs
Tiest
Tvat
Silver-Box
20,000
100..1000
10
100
400
10,000
10,000
W-H
20,000
100..1000
10
100
400
10,000
10,000
RLC
6667
100..1000
3-10
100
1000
3333
3333
Electricity (5)
10,000
48, 96, 192
96
5.10 5
24
400
3360
3360
Traffic (9)
10,000
48, 96, 192
96
5-10-s
24
400
3360
3360
Solar (2)
20,000
48, 96, 192
96
2-10
24
300
5000
5000
The variables correspond to the training sequence length T, the TBPTT subsequence (lookback window) length N, the RNN output dimension Ny, the step size a used in the Adam algorithm, the batch size the test sequence length Tiest, and the validation sequence length Tval- The values in parentheses indicate the considered variate of the respective time series data (selected randomly).
For all experiments, we divide the training sequence length T into S = N + subsequences of length N, leading to the one-step overlap = 1 for all i E S with i 2 2, compare Section 3.1. This allows us leverage the maximum number subsequences during TBPTT training and ultimately ensure a fair comparison between the considered methods and parameter choices. A thorough comparison various options for the number of subsequences and how they might overlap is out of the scope of this paper, as this strongly depends availability of additional prior knowledge and expertise about the data sets in order to decide which sequences actually contain information and are important for training, and which do not. However, this might indeed be relevant in certain scenarios (especially for very long training sequences), where proper selection of the subsequences used for training (and hence the value of S) can potentially significantly accelerate the training procedure without sacrificing overall performance.
We trained the RNNs in standard Python environment using PyTorch and the Adam opti- mizer (Kingma & Ba, 2015). Computations are performed on a standard laptop (Intel Core Ultra 7 265H, 32 GB RAM, NVIDIA RTX PRO 500 Blackwell GPU).
For the training of RNNs for time series forecasting applied to the Electricity. Traffic, and Solar- Energy data sets, we first split the available time series into training segment of length and test segment of length Trest- Then, we pre-train the RNNs on a subsequence of the training data of length T- using different values of N and m. The remaining subsequence of length Tval is used for evaluating the pre-training performance and ultimately selecting m' based on the MSE achieved during pre-training and validation. Then, for every N and m, we re-train the RNN using the full training sequence of length T, and evaluate the forecast performance (and hence choice of m on the test segment length T Trest. Stateful TEPTT and BPTT training is performed directly using the full training sequence of length T.
In Figure 6, we compare the MSE achieved by each RNN after TBPTT training for various combina- tions of the window length N and burn in phase m. The top corresponds to system identification datasets, and the bottom row to time series forecasting data sets. For the latter, we show the training and test MSE achieved after re-training the RNNs on the full training sequence of length T using the previously selected values of N and m.
In all experiments, we find that the performance achieved both the training and test data is influenced by the choices of N and m. as suggested by our theoretical results. Moreover, proper tuning of N leads to a reduction of the MSE achieved on the training and test data of more than 60% in some cases, see Table 1.
24

---

Published as a conference paper at ICLR 2026
Wiemer-Hammerstein
RLC
10
2
660
10
10
1-10
5-000
N-100
5500
220
=11000
N=200
N =10000
10
10a
10
10
10
10
burn-in phase
burn-in phase
Traffic
Solar
N 448
N 48
0.35
N 96
0.3g
N 96
=1192
X
X
0.3
0.2
X
x
0.25
0.1
0000
o
0.2
50
100
150
200
50
100
150
200
burn-in phase m
burn-in phase m
Figure 6: MSE after (re- (re)ntaan the respective RNNs using TBPTT and various combinations the window length N and burn-in phase m. -markrrs correspond to the MSE achieved the training data, 'x mmaakerss to the MSE on the test data. Gray lines refer to BPTT-trained RNN evaluated on the training (solid) and test data (dashed), This figure complements Figure 4.
25


Reviews:
[{"review_id": "0dxjidq0c0", "reviewer": "Reviewer_5iFc", "summary": "This paper propose a method called burn in phase in training RNN using truncated backpropagation through time (TBPTT) that can improve the performance of existing TBPTT method. The back propagation through time training method for RNN is suffered from the recurrent  unrolling. The TBPTT tries to improve the training speed as a cost of accuracy. It divides long sequences into shorter subsequence and increase the number of mini batch which can be parallelized. However, this method will lower the accuracy as each subsequence has a different initial hidden state compare the original sequence. The author propose a method that can reduce this error. The rationale is that for RNN with decaying memory, it will gradually forget the initial condition, thus by using a loss that discard the first m time indices. The performance can be improved. The author provides theoretical proofs for their proposed method. As well as experimental verifications.", "questions": "1. Line 307-308 says that \"Comparing the corresponding solutions therefore allows us to investigate and quantify the effects of using zero initialization in TBPTT compared to using the optimal initializations.\"\n\nWhy is the optimal initialization as defined in equation (11) is minimized over {h_t}.  $h_0$ is the initial condition, which is a free variable that can be optimized. But how does $h_1$ to $h_T$ be optimized? does that depend on the inputs and $\\theta$. I do not understand the state in line 305-306 which says that the additional constraint forces the underlying hidden state sequences to be a subsequence of the same hidden state sequence.\n\n\nIt appears to me that the most optimal setting is each subsequence $i \\in \\{1,...S\\}$ has its own initialization,so the minimization is taken over $S$ number of initial hidden state for the $S$ different subsequences. \n\n\n \n2. It appears to me that Assumption 1 is saying RNN has exponential memory decay. This means that the dependence on initial condition decays exponentially which is characterized by $\\lambda$. Is this equivalent to a more simple assumption which assume the largest eigenvalue of $W_{hh}$ is smaller than 1. And does the $\\lambda$ has any relation with the largest eigenvalue of $W_{hh}$?\n\n\n3. The notation in Theorem 1 means that $\\theta^\\*$ and $\\theta^b$ are any parameters from the parameter space. But the superscript $*$ and $b$ are used before to indicate model using zero initial and optimal initial. The notation here need to be clarified.\n\n4. Is \"batch-averaged output differences\" a standard term? If not, please explain the physical neaning of this quantity.\n\n5. \"As this constant does not depend on N , they must become smaller if N is increased\" What is \"they\" refers here, this argument need to be more \nrigorous \n\n6. The RHS of Theorem 1.2 has a term $\\frac{\\lambda^m}{N-m}$, this function first decays and then increases on interval $(0,N)$. And the best possible $m$ depends on $\\lambda$ and can be explicitly solved. Why line 356 to line 358 only discuss this optimal $m$ vaguely? Is it true that the this optimal $m$ results from the upper bound is not the actually optimal choice of $m$?\n\n7. What is the trade-off discussed in line 380 to line 384. Currently what I can infer from the discussion here is that we want larger $m$, which requires larger $o_{min}$, which results in larger $S$.  This means that larger $m$ results in larger $S$. \n    Firstly, my question is do we always want a larger $m$. \n    Secondly, is the argument here rigorous. Actually when $m$ is large, the upper bound exploded, however, this does not necessarily means that LHS also explode. To be rigorous, I think a lower bound is needed. If not, at least there should be some numerical evidence.\n\n\n8. Regarding the synthetic experiment:\n\n    a. For the synthetic data experiment. How does it verifies your theories. Is the exponential decay rate follows that $\\lambda$ of the unknown system. \n    b. Is it even possible to know the $\\lambda$ of the unknown system.\n    c. You talked about best choice of $m$ depends on $\\lambda$, but from this experiment, can only see larger $m$ the better. What if when $m$ is very large? Does it make the model behave bad as predicated by the upper bound, if not, this means that the bound is too loose.\n    d. Is there any experiment to show the trade-off mention in line 380 to 384.\n\n    e. Why do you conduct experiments on a unknown system, which cannot fully verify your results. Is it more reasonable to conduct experiments on a known system which you can control the $\\lambda$, to fully support your theoretical results.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The proposed method is easy yet effective. The author provides theoretical proofs for their proposed method. \n\nTheorem 1.1 An error bound is developed between the model with different initial hidden states, shows that the difference converges if the model has decaying memory (lambda  small)\nTheorem 1.2 A bound for training regret is developed. This is the difference between ideal case approximation error where an ideal initial state is known and the approximation of using subsequence with burn in.\n\n\nTheorem 2.1 Shows 1 on full sequence \nTheorem 2.2 Shows 2 on full sequence", "weaknesses": "Some of the concept and terms does not have enough motivation, and the intuition and explanation is not very clear for the Theorems, which can cause non expert hard to follow the results.\n\nThere is no recall of notations in Theorem part. Need frequently refer back to find notations.\n\nCertain part of the experiments can be improved. \n\nDetails please see questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3DPgN3d9oE", "reviewer": "Reviewer_4peq", "summary": "This work studies the Truncated Back-propagation Through Time (TBPTT) algorithm for training recurrent neural networks (RNNs) on time series task. They analyze the theoretical bounds on accuracy and performance loss when optimizing over subsequences compared to the full data sequence. Their bounds reveal that the burn-in period m serves as an important hyper-parameter in guaranteeing better performance. RNN hidden states are 0 initialized and up to time m, the network does not generate any output due to negatively affected internal transients. This paper provides theoretical bounds (on regret/performance loss) showing how the burn-in period m influences performance, and shows empirical evidence across synthetic and real-world time-series/system-identification datasets that tuning the burn-in phase improves results.", "questions": "- Other than the contractive dynamics, are there well known dynamics that would satisfy the assumption 1?\n- While the impact of the burn-in period is relatively simple to observe in synthetic data setup, but it comes harder as we evaluate on harder and more challenging real-world time series forecasting problem. Does this imply that the analysis breaks down under this setup?\n- In the case of real-world time series and sysid problems, the experimental setup uses LSTMs as the RNN choice, is there any clue or intuition as to whether the assumption 1 is valid in this setup ? \n- Since typically in TBPTT, the hidden state is carried over from the previous segment, how far off are the empirical numbers compared to this setup? \n- Have you evaluated the effect of N and S in the Tab.1 ? \n- Although, this paper analyzes the regret w.r.t. TBPTT performance with the optimal hidden state initialization and the correct dynamics, how would the analysis change if we were to compare V* with a V_b which measures the performance of the BPTT algorithm (i.e., to predict the future, one can look at the entire past)? \n- Would it be much harder to transfer the theoretical analysis to other state space models?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "- This work analyzes TBPTT with zero initialized hidden states from the lens of regret w.r.t. training benchmark and shows that the burn-in period (m) affects the performance guarantees, both training regret and performance regret.\n- Experimental evaluations on the synthetic data survey as a good validation of the proposed bounds. It shows that larger values of the burn-in phase (m) allow for longer transient phases, and enables RNN to exploit the exploration to reach closer to the benchmark output.", "weaknesses": "- Although the theoretical bounds are valid in the synthetic data setup, the theoretical analysis starts to dwindle down as we reach towards harder real-world setup\n- It is unclear how the assumption 1 holds true in realistic settings with more complex choice of RNNs (such as LSTMs, GRUs, etc.)\n- Although theoretically the burn-in period helps in simplifying and analyzing the regret, it is unclear if the burn-in period is really the bottleneck in more sophisticated TBPTT (such as non-zero hidden state initialization, hidden state propagation between sequences, normalization techniques, etc.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "EWVKq7P9UL", "reviewer": "Reviewer_xtPM", "summary": "This paper assesses the effects of training recurrent neural networks using truncated back-propagation through time (TBPTT). In TBPTT the full training sequence is split into a number of overlapping sub-sequences to which learning is applied. Each sub-sequence includes an initial part (the burn-in phase) that is excluded from the loss. This paper seeks to access the loss in performance produced by approximating BPTT with TBPTT, and the effects of the burn-in period on performance.", "questions": "For all the empirical experiments, why not provide results for BPTT to provide a baseline to compare to?\n\nFor the three real-world data-sets why not perform experiments for a range of different combinations of m and N?\n\nCan you quantify how much TBPTT reduces computational costs? e.g. What is the training time for each method that you compare?\n\nGiven that different hyper-parameters are required to obtain good performance with TBPTT on different tasks, how could hyper-parameter tuning be performed efficiently? Is the time taken to perform hyper-parameter tuning smaller than the time it would take to perform BPTT?\n\nHyper-parameter tuning is a particular concern as the optimal hyper-parameters found for the training data do not seem to produce the best results on the testing data. It would therefore be particularly informative to compare the test-data performance of BPTT with that of TBPTT using the hyper-parameters optimised on the training data.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "The paper combines both theoretical and empirical research.\n\nIt well structured and generally clearly written.", "weaknesses": "The theoretical and practical consequences of this work are unclear.\n\nThe most detailed empirical results are provided for a synthetic task that is so simple it is unlikely to have any practical relevance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)