On the approximation properties of recur- rent encoder-decoder architectures
Year: 2022

Abstract:
Encoder-decoder architectures have recently gained popularity in sequence to sequence modelling, featuring in state-of-the-art models such as transformers. However, a mathematical understanding of their working principles still remains limited. In this paper, we study the approximation properties of recurrent encoder-decoder architectures. Prior work established theoretical results for RNNs in the linear setting, where approximation capabilities can be related to smoothness and memory of target temporal relationships. Here, we uncover that the encoder and decoder together form a particular “temporal product structure” which determines the approximation efficiency. Moreover, the encoder-decoder architecture generalises RNNs with the capability to learn time-inhomogeneous relationships. Our results provide the theoretical understanding of approximation properties of the recurrent encoder-decoder architecture, which precisely characterises, in the considered setting, the types of temporal relationships that can be efficiently learned.

heading: 1 Introduction
text: Encoder-decoder is an increasingly popular architecture for sequence to sequence modelling problems (Sutskever et al., 2014; Chiu et al., 2018; Venugopalan et al., 2015). The core of this architecture is to first encode the input sequence into a vector using the encoder and then map the vector into the output sequence through the decoder. In particular, such architecture forms the main component in the transformer network (Vaswani et al., 2017), which has become a powerful method for modelling sequence to sequence relationships (Parmar et al., 2018; Beltagy et al., 2020; Li et al., 2019). The encoder-decoder family of structures differ significantly from direct application of recurrent neural networks (RNNs, Elman (1990)) and its generalisations (Hochreiter & Schmidhuber, 1997; Cho et al., 2014b) for processing sequences. However, both architectures can be considered as modelling mappings between sequences, albeit with different underlying structures. Hence, a natural but unresolved question is: how are these approaches fundamentally different? Answering this question is not only of theoretical importance but also of practical interest. Currently, architectural selection for different time series modelling tasks is predominantly empirical. Thus, it is desirable to develop a concrete mathematical framework to understand the key differences between separate architectures in order to guide practitioners in a principled way. ∗Equal contribution †Corresponding author
In this paper, we investigate the approximation properties of encoder-decoder architectures. Approximation is one of the most basic and important problems for supervised learning. It considers to what extent the model can fit a target. In particular, we prove a general approximation result in the linear setting, which characterises the types of temporal input-output relationships that can be efficiently approximated by encoder-decoder architectures. These results reveal that such architectures essentially generalise RNNs by lifting the requirement of time-homogeneity (see Remark 3.2) in the target relationships. Hence, it can be used to tackle a broader class of sequence to sequence problems. Furthermore, of particular interest is the identification of a “temporal product structure” — a precise property of the target temporal relationship that highlights another intrinsic difference between recurrent encoder-decoders and RNNs. Our main contributions can be summarised as follows.
1. We prove a universal approximation result for recurrent encoder-decoder architectures in the linear setting, including the approximation rates.
2. We show that in the considered setting, the recurrent encoder-decoder generalises the RNNs and can approximate time-inhomogeneous relationships, which further adapt to additional temporal product structures in the target relationship. This answers precisely how encoder-decoders are different from RNNs, at least in the considered setting.
Organisation. In Section 2, we review the related work on encoder-decoder architectures and general approximation theories of sequence modelling. The approximation problem is formulated in Section 3. Our main results, their consequences and numerical illustrations are presented in Section 4. All the proofs and numerical details are included in appendices.
Notations. For consistency, we adhere to the following notations. Boldfaced letters are reserved for sequences or paths, which can be understood as functions of time. Lower case letters can mean vectors or scalars. Matrices are denoted by capital letters. For α ∈ N, Cα denotes the space of functions with continuous derivatives up to order-α.

heading: 2 Related work
text: We first review some previous works on sequence to sequence modelling. The encoderdecoder architecture first appeared in Kalchbrenner & Blunsom (2013), where they map the input sequence into a vector using convolutional neural networks (CNNs), and then using a recurrent structure to map the vector to the output sequence. With the flexibility of manipulating the underlying structure of encoder and decoder, numerous models based on this architecture have come out thereafter. For instance, Cho et al. (2014b) used gated RNNs as both the encoder and decoder, while in the later work (Cho et al., 2014a), they proposed a CNN-based decoder. In Sutskever et al. (2014), they proposed a deep LSTM for both the encoder and decoder. Bahdanau et al. (2015) first introduced the attention mechanism, which was further developed in the well-known transformer networks (Vaswani et al., 2017). However, most of the research on encoder-decoder architectures focused on applications. A theoretical understanding is helpful for its further improvement and development. From the theoretical point of view, Ye & Sung (2019) studied several theoretical properties of CNN encoder-decoders, including expressiveness, generalisation capability and optimisation landscape. Of particular relevance to the current work is expressiveness, which considers the relationships that can be generated from the architecture. However, this is not approximation. Yun et al. (2020) proved the universal approximation property of transformers for certain classes of functions, for example, permutation equivariant functions, but they did not consider the actual dynamical properties of target relationships that affect approximation. Dynamical proprieties such as memory, smoothness and low rank structures are essential, because they can precisely characterise different temporal relationships and affect the approximation capabilities of models. Assuming the target generated from a hidden dynamical system is one approach, which is widely applied (Maass et al., 2007; Schäfer & Zimmermann, 2007; Doya, 1993; Funahashi & Nakamura, 1993). In contrast, a functional-based approach is
recently introduced, where the target temporal relationships are generated from functionals satisfying specific properties such as linearity, continuity, regularity and time-homogeneity (Li et al., 2021). In Li et al. (2021), the approximation properties of linear RNN models are studied, and the results therein show that the approximation efficiency is related to the memory structure. In Jiang et al. (2021), similar formulations are applied to investigate convolutional architectures, where the results suggest that targets with certain spectrum regularity can be well approximated by dilated CNNs. Under this framework, the target temporal relationship that can be efficiently approximated is characterised by properties such as memory, smoothness and sparsity. This enables us to make precise mathematical comparisons between different architectures. Our results in this work reveal that the encoder-decoders have a special temporal product structure which is intrinsically different from other sequence modelling architectures.

heading: 3 Problem formulation
text: In this section, we precisely define the input space, output space, concept space and hypothesis space, respectively.
Functional formulation of temporal modelling. First, we define the input and output space precisely. A temporal sequence can be viewed as a function of time t. The input space is defined by X = C0((−∞, 0],Rd). This is the space of continuous functions from (−∞, 0] to Rd vanishing at infinity, where d ∈ N+ is the dimension. Denote the element in X by x = {xt ∈ Rd : t ≤ 0}, we equip X with the supremum norm ‖x‖X := supt≤0 ‖xt‖∞. We take the outputs space as Y = Cb([0,∞),R), the space of bounded continuous functions from [0,∞) to R. We consider real-valued outputs, since each dimension can be handled individually for vector-valued outputs. The mapping between input and output sequences can be formulated as a sequence of functionals, i.e. yt = Ht(x), t ≥ 0. The output yt at the time step t depends on the input sequence x. The ground truth relation between inputs and outputs is formulated by the sequence of functionals H := {Ht : t ≥ 0}. We provide an example to illustrate the above formulation. Given an input x, the output y is a smoothed version of x, resulting from convolving x with the Gaussian kernel g(s) =
1√ 2π exp(−
s2 2 ). This relation can be formulated as yt = Ht(x) = ∫∞ 0 g(t+ s)x−sds.
The RNN encoder-decoder model. For the supervised learning problem, our goal is to use a model to learn the target relationship H. First, we define the model. Among all different variants of the encoder-decoder architectures, the RNN encoder-decoder introduced in Cho et al. (2014b) can be considered as the most simple and representative model, where the encoder and decoder are both RNNs. We study this particular model as we try to eliminate other factors and only focus on the encoder-decoder architecture itself. Under our setting, the simplified model of Cho et al. (2014b) with RNNs as both encoder and decoder can be formulated as
hs = σE(WEhs−1 + UExs + bE), v = hτ , gt = σD(WDgt−1 + bD), g0 = v, ot = WOgt + bO,
(1)
where ht, gt are hidden states of the encoder and decoder respectively. Recurrent activation functions are denoted by σE and σD. Here, τ denotes the terminating time step of the encoder, and v is the summary of the input sequence, which is called as the coding vector. The model prediction is denoted as ot ∈ R. All the other notations are model parameters. Equation (1) describes the following model dynamics. First, the encoder reads the entire input x, and then summarises the input into a fixed size coding vector v, which is also the last hidden state of the encoder. Next, the coding vector is passed into the decoder as the initial state, and then the decoder produces an output at each time step. Note that the encoder has a terminating time, and the decoder has a starting time. This is the reason why we take the input and output as semi-infinite sequences.
We study a linear, residual and continuous-time idealisation of the model dynamics (1):
ḣs = Whs + Uxs, v = Qh0, s ≤ 0 ġt = V gt, g0 = Pv, ot = c>gt, t ≥ 0,
(2)
where W ∈ RmE×mE , U ∈ RmE×d, Q ∈ RN×mE , V ∈ RmD×mD , P ∈ RmD×N and c ∈ RmD are parameters. mE and mD denote the width of encoder and decoder, respectively. The coding vector v has dimension N , where we apply linear transformations to control it. We assume h−∞ = 0, which is the usual choice for the initial condition of RNN hidden states. Since our goal is to investigate approximation problems over large time horizons, we are supposed to consider the stable RNN encoder-decoders, where
W ∈ WmE := {W ∈ RmE×mE : eigenvalues of W have negative real parts}, (3) V ∈ VmD := {V ∈ RmD×mD : eigenvalues of V have negative real parts}. (4)
The hypothesis space of RNN encoder-decoder models with arbitrary widths and coding vector dimension is defined as Ĥ := ⋃ mE ,mD,N∈N+ ĤmE ,mD,N , where
ĤmE ,mD,N := { Ĥ := {Ĥt : t ≥ 0} : Ĥt(x) = c>eV tP ∫ ∞ 0 QeWsUx−sds, with
(W,U,Q, V, P, c) ∈ WmE × RmE×d × RN×mE × VmD × RmD×N × RmD } . (5)
The widths mE ,mD and the coding vector dimension N together control the capacity/complexity of the hypothesis space. Note that the assumptions on eigenvalues of W and V ensure that the parameterized linear functionals are continuous. Due to the mathematical form (5), not all functionals can be represented by RNN encoderdecoders. To achieve a good approximation, the target functionals must possess certain structures. We introduce the following definitions to clarify these structures. Definition 3.1. Let H = {Ht : t ≥ 0} be a sequence of functionals.
1. For any t ≥ 0, the functional Ht is linear and continuous if for any λ1, λ2 ∈ R and x1,x2 ∈ X , we have Ht(λ1x1 + λ2x2) = λ1Ht(x1) + λ2Ht(x2), and ‖Ht‖ := supx∈X ,‖x‖X≤1 |Ht(x)| <∞, where ‖Ht‖ denotes the induced functional norm.
2. For any t ≥ 0, the functional Ht is regular if for any sequence {x(n)}∞n=1 ⊂ X such that limn→∞ x(n)s = 0 for almost every s ≤ 0 (Lebesgue measure), we have limn→∞Ht(x(n)) = 0.
For a sequence of functionals H, we define its norm by ‖H‖ := ∫ ∞
0 ‖Ht‖dt.
Remark 3.1. The definitions of linear and continuous functionals are standard. One can view regular functionals as those not determined by inputs on arbitrarily small time intervals, e.g. an infinitely thin spike (i.e. δ-functions).
Given the above definitions, we immediately have the following observation. Proposition 3.1. Let Ĥ ∈ Ĥ be a sequence of functionals in the RNN encoder-decoder hypothesis space (see (5)). Then for any t ≥ 0, Ĥt ∈ Ĥ is a linear, continuous and regular functional. Furthermore, ‖Ĥt‖ decays exponentially as a function of t.
The proof is found in Appendix A. This proposition characterises properties of the encoderdecoder hypothesis space. In particular, it is different from the RNN hypothesis space discussed in Li et al. (2021), since the encoder-decoder is not necessarily time-homogeneous. Remark 3.2. A sequence of functionals H is time-homogeneous if for any t, τ ≥ 0, Ht(x) = Ht+τ (x(τ)), with x(τ)s = xs−τ for all s ∈ R. That is, if the input is shifted to the right by
τ , the output is also shifted by τ . Temporal convolution is an example of time-homogeneous operation (recall the Gaussian convolution discussed in Section 3. An example of timeinhomogeneous relationship is video captioning: shifts in the sequence of input video frames do not necessarily lead to corresponding shifts in the caption text sequence.
Relation with RNNs. Here, we emphasise the differences between the encoder-decoder hypothesis space and the RNN hypothesis space discussed in Li et al. (2021), where Ĥ (RNN) t (x) = ∫∞ 0 c
>eW (t+s)Ux−sds. A key difference is that the encoder-decoder has a structure involving two temporal parameters t and s, while the RNN only has one depending on t + s, due to the time-homogeneity. Owing to this difference and the fact that Ĥ(RNN) ⊂ Ĥ, the encoder-decoder hypothesis space (5) is more general, with the extra capability to learn time-inhomogeneous relationships. Furthermore, eV t and eWs adapt to a temporal product structure, which is an intrinsic difference between encoder-decoders and other architectures. We will discuss this in detail in the next section.

heading: 4 Approximation results
text: One of the most fundamental problems for supervised learning is the approximation problem. It basically concerns the capacity of the hypothesis space to fit the concept space. In general, there are two levels of approximation problems that can be discussed. The first is known as the universal approximation, which considers the density of the hypothesis space in the concept space. The second is the approximation rate, which aims to characterise quantitatively the approximation accuracy concerning the capacity/complexity of the hypothesis space (e.g. the number of trainable parameters). In this section, both of them are developed for RNN encoder-decoders.

heading: 4.1 Universal approximation
text: We first present the most basic density result, which states that any linear, continuous, and regular temporal relationship can be approximated by RNN encoder-decoders up to arbitrary accuracy. The proof is found in Appendix B. Theorem 4.1. Let H be a sequence of linear, continuous and regular functionals defined on X , and satisfy ‖H‖ <∞. Then for any > 0, there exists Ĥ ∈ Ĥ such that
‖H − Ĥ‖ ≡ ∫ ∞
0 ‖Ht − Ĥt‖dt < . (6)
Here, we highlight two important observations while deriving Theorem 4.1. First, one can show that each sequence of functionals H ∈ H can be associated with a unique twoparameter “representation” ρ(t, s), such that Ht(x) = ∫∞ 0 x > −sρ(t, s)ds. Recall the model
form Ĥt(x) = ∫∞ 0 x > −sρ̂(t, s)ds, where ρ̂(t, s) := [c>eV tPQeWsU ]> denotes the corresponding representation. The functional approximation is then reduced to function approximation in the sense of representations, i.e. ‖H−Ĥ‖ ≤ ‖ρ− ρ̂‖L1([0,∞)2). It turns out that ρ directly affects the rate of approximation and gives rise to intrinsic properties. We will discuss this in detail in Section 4.3. In addition, we again emphasise the differences between the present work and Li et al. (2021). In Li et al. (2021), the target relationships are assumed to be time-homogeneous with the representation Ht(x) = ∫∞ 0 ρ(t+ s)x−sds, which only depends on t+s. However, the setting here does not assume time-homogeneity, hence implies a more general representation ρ depending on the two temporal directions t and s simultaneously.

heading: 4.2 General approximation rates
text: While the density result (Theorem 4.1) ensures the universal approximation property of the RNN encoder-decoder, it does not identify targets that can be efficiently approximated. To achieve this, we focus on approximation rates next. We characterise the temporal structure of a target relationship by observing its responses to “constant” input signals. Here, we consider the approximation rates for the model with “large size” coding vector, where the
dimension N ≥ m̄ := min{mE ,mD}. This is the scenario where we fix the widths but take an oversized coding vector. Theorem 4.2. Let H be a sequence of linear, continuous and regular functionals defined on X , and satisfy ‖H‖ < ∞. Consider the output of piece-wise constant signals yci (t, s) = Ht(ei1(−∞,−s]), t, s ≥ 0, i = 1, 2, . . . , d, where {ei}di=1 denotes the standard basis of Rd. Assume that there exist α ∈ N+, β > 0 such that for any i = 1, 2, . . . , d,
yci ∈ Cα+1([0,∞)2), (7)
eβ(t+s) ∂k+l
∂tk∂sl yci (t, s) = o(1) as ‖(t, s)‖ → ∞, (k, l) ∈ N× N+, k + l ≤ α+ 1. (8)
Then for any mE ,mD, N ∈ N+, there exists Ĥ ∈ ĤmE ,mD,N such that
‖H − Ĥ‖ ≤ C(α)γd β2 ( 1 mαE + 1 mαD ) , (9)
where C(α), γ > 0 are both universal constants with dependence only on α and (α, β), respectively, and γ := max
i∈N+, i≤d max k,l∈N, k+l≤α+1 sup t,s≥0
β−(k+l)eβ(t+s) ∣∣∣ ∂k+l∂tk∂sl yci (t, s)∣∣∣ <∞. Here,
the number of trainable parameters is dN(mE +mD) with N ≥ m̄.
The proof is found in Appendix C. First, note that the error bound does not depend on the coding vector size N , as long as N ≥ m̄. This is because further increasing N beyond m̄ only increases the number of trainable parameters, but does not increase the model capacity (see Remark C.2). Only the model widths mE ,mD affect the approximation capabilities. Next, we focus on the classes of target relationships that can be well approximated. Here, α characterises the smoothness of H, and β characterises the temporal decay rates of the output responding to a constant signal under H. This is a notion of memory in the target relationship. The error bound (9) indicates that a sequence of target functionals can be efficiently approximated by the encoder-decoder if it is smooth (large α), and has fast decayed memory (large β). The characterisation in smoothness and memory decay also appears in the approximation results of RNNs (Li et al., 2021), where the upper bound is C(α)γdβmα . However, our results for encoder-decoders suggest extra structures, where the bound involves two (instead of one) temporal parameters together with smoothness and decay memories in both. The two-parameter temporal dependence allows the encoder-decoder to approximate timeinhomogeneous relationships, which generalises the RNN. This two-parameter structure further leads to adaptation to a specific low rank type of target relationships, resulting in finer approximation rates as we discuss next.

heading: 4.3 Approximation rates via temporal product structure
text: Motivation of temporal product structure. In contrast with Theorem 4.2, we next consider the model with N < m̄ = min{mE ,mD}. In this situation, the model has fewer parameters, and we aim to characterise the target relationships by further exploiting the structure of the two-parameter representation ρ(t, s). This leads to a finer approximation rate by considering mE ,mD, N together. We first motivate how the “temporal product structure” arises, and how it relates to the approximation. Detailed discussions and proofs are found in Appendix D. For the illustration purpose, we set the input dimension d = 1. Recall Q ∈ RN×mE , P ∈ RmD×N , then the representation ρ̂ of the encoder-decoder functional can be rewritten as
ρ̂(t, s) = c>eV tP ·QeWsu = N∑ n=1 ( mD∑ i,j=1 ciPjn [ eV t ] ij )( mE∑ i,j=1 uiQnj [ eWs ] ji )
= N∑ n=1 ϕ̂n(t)φ̂n(s). (10)
This is a tensor product structure over the (t, s) time domain (determined by the encoder {φ̂n} and decoder {ϕ̂n} successively). We call it the temporal product structure. As is shown later, this structure significantly affects approximation rates. When {φ̂n} and {ϕ̂n} are selected as the “bases” along s, t direction, respectively, N is considered as the rank of the temporal product. We also define N as the rank of the model, which is understood as the maximum rank of temporal products that the encoder-decoder model can represent.
The rank concept of temporal relationships. Recall that the given number of trainable parameters is dN(mE + mD). Hence, a low rank model may achieve fewer trainable parameters. When investigating relationships that can be well approximated by low rank models, a natural conjecture would be “low rank” targets. What is the meaning of “low rank” for a temporal relationship? It is well-known that in linear algebra, an operator is low rank means that its range space is low-dimensional. This idea can be also applied to temporal relationships. For a “low rank” temporal relationship, the output sequence is more “regular”, meaning that the output sequences (viewed as functions) are in a low-dimensional function space. We provide an intuitive numerical illustration for better understanding.
Figure 1 shows the outputs of a high rank (a) and a low rank (b) target relationship on the same set of random input sequences. Different colours refer to different instances of inputs. In the first case (high rank), the temporal structure of the outputs is very complex and depends sensitively on the inputs. However, in the second case (low rank), the output sequences are much more regular, and only macroscopic structures (e.g. scale/offset) appear. Remark 4.1. In the research of approximation theories for temporal sequences, prior works also related a notion of rank to approximation properties of the dilated convolutional structure (Jiang et al., 2021). Here, we emphasise that the notion of rank considered in our work is very different from that in Jiang et al. (2021), which mainly concerns the tensorisation of a discrete-time sequence according to the width of convolution filters.
POD as an analogue of SVD. Now, we characterise low rank and high rank temporal relationships in a mathematical way. We will introduce the concepts informally, and rigorous definitions and arguments can be found in Appendix D. For a matrix, we can assess its rank by performing the singular value decomposition (SVD). This method can be extended to the temporal relationships using proper orthogonal decomposition (POD; Liang et al. (2002), Berkooz et al. (1993), Chatterjee (2000)). The basic insight is that the function ρ can be decomposed into the following form: ρ(t, s) = ∑N0 n=1 σnϕn(t)φn(s), where N0 ≤ ∞, {ϕn} and {φn} are orthonormal bases, and σ1 ≥ σ2 ≥ · · · ≥ 0 denote the singular values. This procedure can be viewed as applying SVD to an infinite-dimensional space (when N0 = ∞). An analogue of Eckart–Young theorem (Eckart & Young, 1936), which characterises the best low rank approximation, also holds for POD. It roughly states that infrank(ρ̂)=N‖ρ − ρ̂‖2L2 = ∑N0 n=N+1 σ 2 n. That is,
any target ρ has a rank-N best approximation, with error equalling to the tail sum of the squared singular values. In other words, a target with fast decayed σn (low “effective rank”) has smaller approximation errors. This forms the basis of our next result, which states that if the target relationship possesses an effective low rank structure in terms of the decay of singular values, then one can achieve an efficient approximation using encoder-decoder structures by limiting the size of coding vectors. Detailed definitions for {σn} and proofs are found in Appendix D. Theorem 4.3. Assume the same conditions as in Theorem 4.2. Then for any mE ,mD, N ∈ N+ with N ≤ m̄, there exists Ĥ ∈ ĤmE ,mD,N such that
‖H − Ĥ‖ . C(α)γd β2
{( 1 + √ m̄−N ) · (
1 mαE + 1 mαD
) + (
m̄∑ n=N+1 σ2n
)1/2
+ (
m̄∑ n=N+1 σn
)1/2 · ( 1
m α/2 E + 1 m α/2 D
)} , (11)
where . hides universal positive constants, and m̄ = min{mE ,mD}. Here, the number of trainable parameters is dN(mE +mD) with N ≤ m̄.
This is a finer approximation rate compared to Theorem 4.2, where both the widths mE ,mD and the coding vector size N affect the model capacity for approximation. Besides the smoothness and memory decay, we have the additional rank structure of the target relationship, which is characterised by its singular values {σn}. We again focus on the class of functionals that can be well approximated. Smoothness α and decay rate β is the same as Theorem 4.2. The difference lies in the rank structure indicated by {σn}: the error bound is small if {σn} has a small tail ∑m̄ n=N+1 σ 2 n. It suggests that a target with fast decayed {σn} or low “effective rank” can be well approximated by the RNN encoder-decoder with fewer parameters. Due to the Eckart–Young-like low rank approximation, we can appropriately select N based on the decay rate of singular values. Here, we emphasise that the temporal product is an intrinsic structure arising from the encoder-decoder architecture. Recall the dynamics of the encoder-decoder: it first encodes the input sequence into a coding vector, and then decodes an entire output sequence from it. In this sense, the coding vector is the only interaction between the input and output. Thus, the coding vector size N is an essential measure of the model capacity concerning the dependence of outputs on inputs. Here, we show that this concept can be formalised as a notion of rank, which can pinpoint the precise types of input-output relationships that encoder-decoder architectures are well adapted to.
Numerical illustrations. Here, we utilise numerical examples to illustrate the above discussions. We observe how the decay rate of singular values, the rank N0 of the target relationships, and the model rank N affect the approximation error ‖H − Ĥ‖. However, it is not always possible to construct the best approximation. Instead, we perform some training steps to achieve an upper bound of the approximation error, which is consistent with our theoretical results. In Figure 2, we train linear encoder-decoder models to learn three relationships of different ranks determined by various decay patterns of singular values, given in (a), (b) and (c). Different colours denote targets with different ranks. From Figure 2, we have the following observations consistent with previous discussions. First, observe that increasing the model rank N makes approximation errors smaller, as expected. Moreover, note that when increasing N , the speeds of error decrements are different. If the singular values decay fast, the approximation errors also decay fast. This implies that a target with fast decayed singular values can be approximated efficiently with fewer parameters (smaller N). In addition, for each experiment, we are able to achieve low approximation errors by choosing N m. The errors will remain unchanged or decrease much more slowly when further increasing N . This suggests that in practice, one can choose N such that it covers the major singular values of the target in order to improve the approximation efficiency.
In Figure 3, we perform experiments on the forced Lorentz 96 system (Lorenz, 1996), which parameterises a high-dimensional and nonlinear relationship between input forcing and model states. The parameters K,J in the Lorenz 96 system control the overall complexity of the target (see Appendix E.3 for details). We use the RNN encoder-decoder with tanh activations to learn this target. Although our theories are developed in the linear regime, the low rank approximation phenomenon also appears in this nonlinear setting. The error decrements saturate when increasing the coding vector size N beyond a threshold, suggesting the existence of some implicit notion of “rank” of the target nonlinear functional. This “rank” increases with the target complexity (mainly K).

heading: 5 Conclusion
text: We theoretically study the approximation properties of the RNN encoder-decoder in a linear setting. We prove a universal approximation result for linear temporal relationships utilising encoder-decoder architectures, and show that they generalise RNNs to the timeinhomogeneous setting. Moreover, we discover an important temporal product structure that characterises the types of input-output relationships especially suited for the efficient approximation using encoder-decoders. This elucidates the key differences between these novel architectures and classical methods for temporal modelling, and forms a basic step towards understanding the intricacies of modern deep learning.
Reproducibility statements. Detailed proofs for theoretical results, and complete settings of numerical examples are found in the appendix. The source code for numerical tests can be made available upon request. Here is a quick reference:
Proposition 3.1 Properties of RNN encoder-decoder functionals Appendix A Theorem 4.1 Universal approximation theorem Appendix B Theorem 4.2 General approximation rates Appendix C Theorem 4.3 Approximation rates concerning temporal product structure Appendix D
Figure 1 Illustration of high/low rank temporal relationships Appendix E.1 Figure 2 Numerical examples on singular values Appendix E.2 Figure 3 Numerical examples on Lorenz 96 systems Appendix E.3

heading: Acknowledgements
text: ZL is supported by Peking University under BICMR mathematical scholarship. HJ is supported by National University of Singapore under PGF scholarship. QL is supported by the National Research Foundation, Singapore, under the NRF fellowship (NRF-NRFF13-20210005).

heading: A Properties of model functionals
text: In this section, we prove observations of the hypothesis space reported in Proposition 3.1.
Proof of Proposition 3.1. Recall that Ĥt(x; θ) = ∫ ∞
0 c>eV tMeWsUx−sds. (12)
Fix any θ = (W,V,U,M, c). The linearity is obvious. Since both W ∈ WmE and V ∈ VmD have eigenvalues with negative real parts, there exist c1, c2, c′1, c′2 > 0, such that ‖eV t‖∞ ≤ c1e −c2t and ‖eWs‖∞ ≤ c′1e−c ′ 2s for any t, s ≥ 0, hence
|c>eV tMeWsUx−s| ≤ ‖c‖1‖eV tMeWsUx−s‖∞ ≤ ‖c‖1‖eV t‖∞‖M‖∞‖eWs‖∞‖U‖∞‖x−s‖∞ . e−c2te−c ′ 2s‖x‖X , (13)
where . hides universal positive constants depending on parameters θ. Therefore |̂Ht(x; θ)| . e−c2t‖x‖X ⇒ ‖Ĥt(·; θ)‖ . e−c2t. (14)
That is, the functional Ĥt(·; θ) is bounded (i.e. continuous) with an exponentially-decayed norm (as a function of t). Finally, by (13) and Lebesgue’s dominated convergence theorem, the functional Ĥt(·; θ) is also regular. The proof is completed.

heading: B Universal approximation
text: In this section, we provide the proof of Theorem 4.1, i.e. the universal approximation property of RNN encoder-decoders. As is stated in the main text, the key step is to utilise the classical representation theorem, which helps us to reduce the approximation problem of functionals to functions.

heading: B.1 Preliminaries
text: First, we list the background definitions and notations used in the following theorems. Let (X,A) be a measure space (with A as the σ-algebra of subsets of X).
• The space X is called locally compact if for any x ∈ X, x has a compact neighbourhood.
• X is a Hausdorff space if all distinct points in X are pair-wisely separable by neighbourhoods. That is, for any x, y ∈ X, there exists a neighbourhood ∆x of x and a neighbourhood ∆y of y, such that ∆x ∩∆y = ∅.
• The measure µ is called a finite measure, if it satisfies µ(X) < ∞. The measure ν is called a σ-finite measure, if X can be covered with at most countably many measurable sets with finite measure. That is, there are measurable sets {An}∞n=1 ⊂ A with ν(An) < ∞ for all n ∈ N+, such that ⋃∞ n=1An = X. Obviously, a finite
measure is also σ-finite. • LetX = (−∞, 0]. For a measure µ on the measure space ((−∞, 0],A), µ is absolutely
continuous with respect to the Lebesgue measure ν, if for every measurable set A, ν(A) = 0 implies µ(A) = 0, which is written as µ ν.
Denote by C0(X) the linear space of continuous functions defined on X vanishing at infinity. We have the following classical representation theorem. Theorem B.1 (Riesz-Markov-Kakutani representation theorem). Let X be a locally compact Hausdorff space. For any continuous linear functional ψ on C0(X), there is a unique, regular, countably additive and signed measure µ on X, such that
ψ(f) = ∫ X f(x)dµ(x), ∀f ∈ C0(X), (15)
with ‖ψ‖ = |µ|(X). Here, |µ|(X) denotes the total variation of (the signed measure) µ, which is defined as |µ|(X) := supP ∑k i=1 |µ(Ai)|, where P : X = ⋃k i=1Ai is a partition over

heading: X with Ai ∈ A for all i = 1, 2, · · · , k.
text: Proof. Well-known, see e.g. Bogachev (2007) (CH 7.10.4).
Remark B.1. It is straightforward to verify that |µ|(X) = supA∈A(|µ(A)|+ |µ(Ac)|). Furthermore, if µ has a density dµ/dν with respect to a countably additive, nonnegative measure ν, then we have |µ|(X) = ‖dµ/dν‖L1(ν).
To handle signed measures, the following Jordan decomposition theorem (Bogachev, 2007) is necessary. Theorem B.2. Let µ be a signed measure on the measure space (X,A). Then, there are two mutually singular (non-negative) measures µ+ and µ− on (X,A), such that µ = µ+ − µ−. Moreover, such a pair (µ+, µ−) is unique.
Based on this, we have the following proposition to characterise absolutely continuous signed measures. Proposition B.1. If µ and ν are signed measures, then we have µ ν ⇔ µ+ ν and µ− ν.
We also need the following Radon-Nikodym theorem (Bogachev, 2007). Theorem B.3. Let (X,A, ν) be a σ-finite measure space, and let µ be a σ-finite signed measure, such that µ ν. Then there exists a unique measurable function f , such that µ(A) = ∫ A fdν for every measurable set A.

heading: B.2 Proofs
text: Before we prove the universal approximation theorem (Theorem 4.1), we need some lemmas. Lemma B.1. Let {Ht : t ≥ 0} be a family of linear, continuous and regular functionals defined on X . Then there exists a integrable function ρ : [0,∞)2 → Rd, i.e.
‖ρ‖L1([0,∞)2) := d∑ i=1 ‖ρi‖L1([0,∞)2) <∞, (16)
such that
Ht(x) = ∫ ∞
0 x>−sρ(t, s)ds, ∀x ∈ X . (17)
In particular, we have ‖H‖ = ∫∞
0 ‖Ht‖dt = ‖ρ‖L1([0,∞)2).
Proof. Obviously, (−∞, 0] is a locally compact Hausdorff space. For any t ≥ 0, since Ht is linear continuous, according to the Riesz-Markov-Kakutani representation theorem (Theorem B.1), there exists a unique, regular, countably additive and signed measure µt, such that
Ht(x) = ∫ 0 −∞ x>s dµt(s), ∀x ∈ X , (18)
with ∑d i=1 |µt,i|((−∞, 0]) = ‖Ht‖. We show that for any t ≥ 0, i = 1, 2, · · · , d, µt,i is absolutely continuous with respect to ν (the Lebesgue measure), i.e. µt,i ν. According to Theorem B.2 and Proposition B.1, one can assume µt,i to be non-negative without loss of generality. Take a measurable set A ⊂ (−∞, 0] with ν(A) = 0, the aim is to show µt,i(A) = 0. Let A′ = (−∞, 0] \ A. Since both A and A′ are measurable, there exist Kn ⊂ A, K ′n ⊂ A′ with Kn,K ′n closed, such that µt,i(A \Kn) ≤ 1/n, µt,i(A′ \K ′n) ≤ 1/n
and ν(A′ \K ′n) ≤ 1/n for any n ∈ N+. Fix any i ∈ {1, 2, . . . , d}, we construct the sequence of input signals {x(n)}∞n=1 as
x (n) s,j =  0, s ≤ 0, j 6= i, 0, s ∈ K ′n, j = i, j = 1, 2, . . . , d, 1, s ∈ Kn, j = i,
(19)
which can then be continuously extended to (−∞, 0] by defining x(n)s,i := d(s,K′n)
d(s,Kn)+d(s,K′n) ∈
[0, 1]. 1
We deduce that limn→∞ x(n)s,i = 0 for ν-a.e. s ≤ 0. In fact, let S := {s ≤ 0 : limn→∞ x (n) s,i = 0}, we have K ′n ⊂ S since for any s ∈ K ′n, x (n) s,i = 0. Hence, (−∞, 0] \ S ⊂ A ∪ (A′ \K ′n), which gives ν((−∞, 0]\S) ≤ ν(A)+ν(A′ \K ′n) ≤ 1/n→ 0 as n→∞. Due to the regularity of Ht, we get limn→∞Ht(x(n)) = 0. By (18) and (19), we have
Ht(x(n)) = d∑ j=1 ∫ 0 −∞ x (n) s,j dµt,j(s) = ∫ 0 −∞ x (n) s,i dµt,i(s)
= ∫ Kn + ∫ A\Kn + ∫ A′\K′n x (n) s,i dµt,i(s) = µt,i(Kn) + I1,n + I2,n, (20)
where µt,i(Kn) = µt,i(A) − µt,i(A \ Kn) ∈ [µt,i(A) − 1/n, µt,i(A)], and |I1,n| + |I2,n| ≤ ∫ A\Kn + ∫ A′\K′n
1dµt,i(s) = µt,i(A \ Kn) + µt,i(A′ \ K ′n) ≤ 2/n, which gives limn→∞Ht(x(n)) = µt,i(A). Therefore, µt,i(A) = 0. Notice that |µt,i((−∞, 0])| ≤ |µt,i|((−∞, 0]) ≤ ‖Ht‖ < ∞ for a.e. t ≥ 0 (since ‖H‖ =∫∞
0 ‖Ht‖dt <∞), we get that ((−∞, 0],A, µt,i) is a finite measure space, and hence σ-finite. Obviously, ((−∞, 0],A, ν) is a σ-finite measure space. According to the Radon-Nikodym theorem (Theorem B.3), there exists a unique measurable function ρt,i : (−∞, 0]→ R, such that µt,i(A) = ∫ A ρt,i(s)dν(s) for every measurable set A. Hence, we have
Ht(x) = ∫ 0 −∞ x>s ρt(s)ds = ∫ ∞ 0 x>−sρ(t, s)ds, ∀x ∈ X , (21)
with ρ(t, s) := ρt(−s). In addition, by Remark B.1, we have |µt,i|((−∞, 0]) =∫ 0 −∞ |ρt,i(s)|ds, which gives
‖H‖ = ∫ ∞
0 ‖Ht‖dt = d∑ i=1 ∫ ∞ 0 |µt,i|((−∞, 0])dt = d∑ i=1 ‖ρi‖L1([0,∞)2) = ‖ρ‖L1([0,∞)2). (22)
The proof is completed.
Based on this representation theorem, the problem of functional approximation is reduced as function approximation. That is,
‖H − Ĥ‖ = ∫ ∞
0 ‖Ht − Ĥt‖dt = ∫ ∞ 0 sup ‖x‖X≤1 ∣∣Ht(x)− Ĥt(x; θ)∣∣dt = ∫ ∞
0 sup ‖x‖X≤1 ∣∣∣∣∫ ∞ 0 x>−s(ρ(t, s)− ρ̂(t, s))ds ∣∣∣∣dt
≤ ∫ ∞
0 sup ‖x‖X≤1 ∫ ∞ 0 ‖x−s‖∞‖ρ(t, s)− ρ̂(t, s)‖1dsdt
≤ d∑ i=1 ∫ ∞ 0 ∫ ∞ 0 |ρi(t, s)− ρ̂i(t, s)|dsdt, (23)
1Here, d(s,B) := inf{|s− a| : a ∈ B} is the distance between a point s and a set B.
i.e.
‖H − Ĥ‖ ≤ ‖ρ− ρ̂‖L1([0,∞)2) := d∑ i=1 ‖ρi − ρ̂i‖L1([0,∞)2). (24)
Lemma B.2. Let ρ(t, s) : [0,∞)2 → R with ‖ρ‖L1([0,∞)2) < ∞. Then for any > 0, there exists a polynomial p(u, v) = ∑m j=1,k=1 cjku
jvk, such that∫ ∞ 0 ∫ ∞ 0 |ρ(t, s)− p(e−t, e−s)|dtds < . (25)
Proof. Fix any > 0. Consider the following transformation R(u, v) = { 1 uvρ(− ln u,− ln v), u, v ∈ (0, 1], 0, uv = 0. (26)
This transformation preserves the norm with ‖ρ‖L1([0,∞)2) = ‖R‖L1([0,1]2).
First, according to the density of continuous functions in Lp space (Rudin, 1987, Theorem 3.14), there exists R̃ ∈ C([0, 1]2), such that ‖R− R̃‖L1([0,1]2) < /2. Next, by the density of polynomials in the space of continuous functions (Lorentz, 2005, Theorem 6), there exists a polynomial q(u, v) = ∑m j,k=0 cjku
jvk, such that ‖R̃ − q‖L∞([0,1]2) < /2. Finally, let p(u, v) = uvq(u, v), we have∫ ∞
0 ∫ ∞ 0 |ρ(t, s)− p(e−t, e−s)|dtds = ∫ 1 0 ∫ 1 0 ∣∣∣∣R(u, v)− 1uv p(u, v) ∣∣∣∣dudv
= ‖R− q‖L1([0,1]2) ≤ ‖R− R̃‖L1([0,1]2) + ‖R̃− q‖L∞([0,1]2) < /2 + /2 = , (27)
which completes the proof.
Now we are ready to prove the universal approximation theorem.
Proof of Theorem 4.1. According to the representation theorem (Lemma B.1), we have Ht(x) = ∫ ∞
0 x>−sρ(t, s)ds, ∀x ∈ X , (28)
where ‖ρ‖L1([0,∞)2) = ‖H‖ < ∞. Therefore, by Lemma B.2, there exists pi(u, v) =∑m j,k=1 c (i) jku
jvk, i = 1, 2, . . . , d, where m is the maximal degree of {pi}di=1, such that d∑ i=1 ∫ ∞ 0 ∫ ∞ 0 |ρi(t, s)− pi(e−t, e−s)|dtds < . (29)
Let c = u = 1m, V = W̃ = −diag (1, 2, . . . ,m), W = diag(W̃ , W̃ , · · · , W̃ ) ∈ Rdm×dm, U = diag(u, u, · · · , u) ∈ Rdm×d, (30)
M = PQ = (M1,M2, · · · ,Md) ∈ Rm×dm, [Mi]jk = c(i)jk , we get
ρ̂(t, s)> = c>eV tPQeWsU = c>eV tMeWsU = c>eV t · (M1,M2, · · · ,Md) · diag(eWs, eWs, · · · , eWs) · diag(u, u, · · · , u) = (c>eV tM1, c>eV tM2, · · · , c>eV tMd) · diag(eWsu, eWsu, · · · , eWsu) = (c>eV tM1eWsu, c>eV tM2eWsu, · · · , c>eV tMdeWsu), (31)
with ρ̂i(t, s) = c>eV tMieWsu = pi(e−t, e−s), i = 1, 2, . . . , d. (32)
Therefore, by (24) and (29), we have
‖H − Ĥ‖ ≤ d∑ i=1 ‖ρi − ρ̂i‖L1([0,∞)2)
= d∑ i=1 ∫ ∞ 0 ∫ ∞ 0 ∣∣ρi(t, s)− pi(e−t, e−s)∣∣dsdt < , (33) which completes the proof.

heading: C General approximation rates
text: In this section, the proof of Theorem 4.2 is given. Again, by (24), the aim now is to investigate the function approximation ‖ρ− ρ̂‖. Since one can handle each spatial dimension separately (similarly with (30) and (31)), we firstly derive the estimates by assuming d = 1, and then extend the obtained results to the case of multi-dimensional inputs (for general d ∈ N+).
Conditions on representation. To characterise the accuracy of using the model c>eV tMeWsu (with M := PQ) to approximate the target ρ(t, s), the first stuff is to translate the conditions on the output (of piece-wise constant signals) to the representation. Recall that yc(t, s) = Ht(1(−∞,−s]), t, s ≥ 0, we get ρ(t, s) = − ddsHt(1(−∞,−s]). Hence, the assumptions on yc in Theorem 4.2 is equivalent to the following smoothness and exponential decay conditions on ρ. That is, there exist α ∈ N+, β > 0 such that
ρ ∈ Cα([0,+∞)2), (34)
eβ(t+s) ∂k+l
∂tk∂sl ρ(t, s) = o(1) as ‖(t, s)‖ → ∞, k, l ∈ N, k + l ≤ α. (35)
Note that the last decay condition implies
sup t,s≥0 β−(k+l)eβ(t+s) ∣∣∣∣ ∂k+l∂tk∂sl ρ(t, s) ∣∣∣∣ ≤ γ, k, l ∈ N, k + l ≤ α (36) for some γ > 0.

heading: C.1 Basics
text: Let Ω ∈ Rd be a bounded set. Define the spaces Cα(Ω) := {f ∈ C(Ω̄) : Dif ∈ C(Ω̄) for all |i| ≤ α}, α ∈ N, (37)
Cα,µ(Ω) := {f ∈ Cα(Ω) : |Dif(x)−Dif(y)| ≤ K‖x− y‖µ2 for some K > 0, for all x, y ∈ Ω and |i| = α}, (38)
and the “norm”
|f |α,µ,Ω := sup |i|=α sup x,y∈Ω |Dif(x)−Dif(y)| ‖x− y‖µ2 , ∀f ∈ Cα,µ(Ω), (39)
with the shorthand ‖ · ‖Ω := | · |0,0,Ω. Theorem C.1 (Multivariate Jackson’s theorem (Schultz, 1969), Theorem 4.10). Let Ω ∈ Rd be a regular, 2 bounded and open set, and f ∈ Cα,µ(Ω̄) for some α ∈ N, µ ∈ [0, 1]. Then for any n ∈ N+, we have
inf p∈Pdn ‖f − p‖Ω̄ ≤ C(α, µ) nα+µ |f |α,µ,Ω̄, (40)
2It is proved that every bounded, open and convex set is regular. See Morrey (1966) (Lemma 3.4.1).
where Pdn denotes the set of all polynomials with the degree of no more than n in each variable, C(α, µ) > 0 is a universal constant only depending on α, µ and Ω.
A commonly used case is when µ = 0. That is, for f ∈ Cα(Ω̄), we get
|f |α,0,Ω̄ ≤ 2 max|i|=α ‖D if‖L∞(Ω̄) <∞. (41)
For any x, x0 ∈ Ω̄, let p̃(x) := p(x) + f(x0)− p(x0), and ẽ(x) := f(x)− p̃(x). Then p̃ ∈ Pdn, and
|f(x)− p̃(x)| ≤ |ẽ(x0)|+ |ẽ(x)− ẽ(x0)| ≤ sup x,y∈Ω̄ |ẽ(x)− ẽ(y)| = |ẽ|0,0,Ω̄ = ‖f − p̃‖Ω̄, ∀x ∈ Ω̄. (42)
This gives the following convenient corollary. Corollary C.1. Let Ω ∈ Rd be a regular, bounded and open set, and f ∈ Cα(Ω̄) for some α ∈ N. Then for any n ∈ N+, there exists p ∈ Pdn such that
‖f − p‖L∞(Ω̄) ≤ Cα nα max |i|=α ‖Dif‖L∞(Ω̄), (43)
where Pdn denotes the set of all polynomials with the degree of no more than n in each variable, Cα > 0 is a universal constant only depending on α and Ω.

heading: C.2 Proofs
text: Now we are ready to present the proof.
Proof of Theorem 4.2. Step 1: domain transform. Consider the transform from the infinite domain [0,∞)2 to the compact one [0, 1]2:
R(u, v) = { 1 uvρ(−c0 ln u,−c0 ln v), u, v ∈ (0, 1], 0, uv = 0, (44)
where c0 := (α+ 1)/β > 0 is a fixed constant. A straightforward computation by induction shows that, for any k, l ∈ N, k + l ≤ α, and any u, v ∈ (0, 1],
∂k+l ∂uk∂vl R(u, v) = (−1) k+l uk+1vl+1 k∑ i=0 l∑ j=0 C(k, i)C ′(l, j)ci+j0 ∂i+j ∂ti∂sj ρ(−c0 ln u,−c0 ln v), (45)
where C(k, i), C ′(l, j) are some integer constants, and (t, s) = (−c0 ln u,−c0 ln v) is a oneto-one mapping between (0, 1]2 and [0,∞)2. By (36), we get∣∣∣∣ ∂k+l∂uk∂vlR(u, v) ∣∣∣∣ ≤ 1uk+1vl+1 k∑ i=0 l∑ j=0 |C(k, i)||C ′(l, j)|ci+j0
∣∣∣∣ ∂i+j∂ti∂sj ρ(−c0 ln u,−c0 ln v) ∣∣∣∣ ,∣∣∣∣ ∂k+l∂uk∂vlR(e− tc0 , e− sc0 ) ∣∣∣∣ ≤ e (k+1)c0 te (l+1)c0 s k∑ i=0 l∑ j=0 |C(k, i)||C ′(l, j)|ci+j0 ∣∣∣∣ ∂i+j∂ti∂sj ρ(t, s) ∣∣∣∣
≤ k∑ i=0 l∑ j=0 |C(k, i)||C ′(l, j)|(α+ 1)i+j · β−(i+j)eβ(t+s) ∣∣∣∣ ∂i+j∂ti∂sj ρ(t, s) ∣∣∣∣ ≤ k∑ i=0 l∑ j=0 |C(k, i)||C ′(l, j)|(α+ 1)i+jγ ≤ C(α)γ, (46)
where C(α) > 0 is a universal constant only depending on α. According to the decay condition (35), we get
lim (u,v)→(0,0)
∂i+j
∂ti∂sj ρ(−c0 ln u,−c0 ln v) uk+1vl+1 = lim (t,s)→(+∞,+∞) e (k+1) c0 te (l+1) c0 s ∂ i+j ∂ti∂sj ρ(t, s)
= lim (t,s)→(+∞,+∞)
e (k−α) α+1 βte (l−α) α+1 βs · eβ(t+s) ∂
i+j
∂ti∂sj ρ(t, s)
= 0, (47) and similarly for u0, v0 ∈ (0, 1],
lim (u,v)→(u0,0)
∂i+j
∂ti∂sj ρ(−c0 ln u,−c0 ln v) uk+1vl+1 = lim (t,s)→(−c0 lnu0,+∞) e (k+1) c0 te (l+1) c0 s ∂ i+j ∂ti∂sj ρ(t, s)
= lim ‖(t,s)‖→+∞
e (k−α) α+1 βte (l−α) α+1 βs · eβ(t+s) ∂
i+j
∂ti∂sj ρ(t, s)
= 0, (48)
lim (u,v)→(0,v0)
∂i+j
∂ti∂sj ρ(−c0 ln u,−c0 ln v) uk+1vl+1 = 0. (49)
This gives ∂k+l
∂uk∂vl R(u, v) = 0, (u, v) ∈ [0, 1]× {0} ∪ {0} × [0, 1], (50)
and
M0 := max k,l∈N, k+l≤α max (u,v)∈[0,1]2 ∣∣∣∣ ∂k+l∂uk∂vlR(u, v) ∣∣∣∣ ≤ C(α)γ (51)
by (46) and (50). Hence, R(u, v) ∈ Cα([0, 1]2) with bounded derivatives. Step 2: polynomial approximation. According to Corollary C.1 and (51), we obtain that there exists R̃n ∈ P2n, such that
‖R− R̃n‖L∞([0,1]2) ≤ Cα nα max k,l∈N, k+l=α ∥∥∥∥ ∂k+l∂uk∂vlR(u, v) ∥∥∥∥ L∞([0,1]2)
≤ C(α)γ nα , ∀n ∈ N+, (52)
where C(α) > 0 is a universal constant only related to α. Furthermore, let R̂n(u, v) := R̃n(u, v)− R̃n(u, 0)− R̃n(0, v) + R̃n(0, 0), we get R̂n ∈ P2n with R̂n(u, 0) = R̂n(0, v) = 0 for any u, v ∈ [0, 1]. By (50), we have R(u, v) = 0 for any (u, v) ∈ [0, 1]×{0}∪{0}× [0, 1], then
‖R− R̂n‖L∞([0,1]2) ≤ ‖R− R̃n‖L∞([0,1]2) + ‖R̃n − R̂n‖L∞([0,1]2) ≤ ‖R− R̃n‖L∞([0,1]2) + ‖R̃n(u, 0)−R(u, 0)‖L∞([0,1]2)
+ ‖R̃n(0, v)−R(0, v)‖L∞([0,1]2) + ‖R̃n(0, 0)−R(0, 0)‖L∞([0,1]2) . ‖R− R̃n‖L∞([0,1]2), (53)
i.e. we can further require the approximator satisfying the zero half-boundary condition (i.e. vanishing on (u, v) ∈ [0, 1] × {0} ∪ {0} × [0, 1]) without effecting the approximation accuracy. Let m̄ := min{mE ,mD}, we get
‖R− R̃‖L∞([0,1]2) ≤ C(α)γ m̄α
≤ C(α)γ (
1 mαE + 1 mαD
) , (54)
where
R̃ := R̃m̄ ∈ P2m̄, R̃(u, v) := m̄∑ i=1 m̄∑ j=1 r̃iju ivj . (55)
Let
c = 1m̄, u = 1m̄, M = [r̃ij ] ∈ Rm̄×m̄, (56) V = −diag(2, 3, · · · , m̄+ 1)/c0, W = −diag(2, 3, · · · , m̄+ 1)/c0, (57)
then by (24) and (54), we have
‖H − Ĥ‖ ≤ ‖ρ− ρ̂‖L1([0,∞)2) (58) = ∥∥ρ(t, s)− c>eV tMeWsu∥∥
L1([0,∞)2) = ∥∥∥ρ(t, s)− e− tc0 e− sc0 R̃(e− tc0 , e− sc0 )∥∥∥
L1([0,∞)2) = c20 ∥∥R− R̃∥∥ L1([0,1]2) (59)
≤ c20 ∥∥R− R̃∥∥
L∞([0,1]2) ≤ C(α)γ β2m̄α ≤ C(α)γ β2 ( 1 mαE + 1 mαD ) . (60)
That is, one can achieve an approximation accuracy scaling like (1/m̄)α with m̄2 parameters. The proof is completed.
Remark C.1. The extension to multi-dimensional inputs (general d ∈ N+) is found in the last paragraph of Appendix D.2. Remark C.2. Recall M = PQ ∈ RmD×mE with P ∈ RmD×N , Q ∈ RN×mE , we only need to investigate the case of N ≤ min{mE ,mD} = m̄, since rank(M) ≤ m̄.

heading: D Approximation rates via temporal product structure
text: In this section, the proof of Theorem 4.3 is provided.

heading: D.1 Proper Orthogonal Decomposition
text: Proper orthogonal decomposition (POD; (Liang et al., 2002), (Berkooz et al., 1993), (Chatterjee, 2000)) is a method for model reduction, which is commonly applied to numerical PDEs and fluids mechanics. It can be viewed as an extension of singular value decomposition (SVD) and principal component analysis (PCA) to infinite-dimensional spaces. Fix any R ∈ L∞([0, 1]2). 3 Define the POD operator
K : φ(v) 7→ ∫ 1
0 ∫ 1 0 R(u, v)φ(v)dv ·R(u, v)du, φ(v) ∈ L2[0, 1]. (61)
Proposition D.1. The operator K is linear, bounded, compact, self-adjoint and nonnegative.
Proof. (i) The linearity is obvious. (ii) Let
R : φ(v) 7→ ∫ 1
0 R(u, v)φ(v)dv, φ(v) ∈ L2[0, 1], (62)
then
(Kφ)(v) = ∫ 1
0 R(u, v)(Rφ)(u)du, φ(v) ∈ L2[0, 1]. (63)
By the Cauchy-Schwartz inequality, we get
‖Rφ‖L2[0,1] ≤ ‖R‖L2([0,1]2)‖φ‖L2[0,1], (64) 3Here, we use the same notation as (44), since the function defined there is also bounded.
which gives ‖Kφ‖L2[0,1] ≤ ‖R‖L2([0,1]2)‖Rφ‖L2[0,1] ≤ ‖R‖2L2([0,1]2)‖φ‖L2[0,1]. (65)
Note that R(u, v) ∈ L∞([0, 1]2) ⊂ L2([0, 1]2), hence both K and R are bounded operator from L2[0, 1] to itself. (iii) It is well-known that the Hilbert–Schmidt integral operator
(Cψ)(v) = ∫ 1
0 R(u, v)ψ(u)du, ψ(u) ∈ L2[0, 1] (66)
is a compact operator from L2[0, 1] to itself. Therefore
Kφ = CRφ, φ ∈ L2[0, 1], (67) which gives K = CR. Since R is bounded and C is compact, we get K is also compact. (iv) By Fubini’s theorem, it is straightforward to verify that
〈Kφ, ψ〉L2[0,1] = ∫ 1
0 ∫ 1 0 R(u,w) (∫ 1 0 R(u, v)φ(v)dv ) du · ψ(w)dw
= ∫ 1
0 ∫ 1 0 ∫ 1 0 R(u,w)R(u, v)φ(v)ψ(w)dvdudw
= ∫ 1
0 ∫ 1 0 R(u, v) (∫ 1 0 R(u,w)ψ(w)dw ) du · φ(v)dv
= 〈φ,Kψ〉L2[0,1], φ, ψ ∈ L2[0, 1]. (68)
(v) By Fubini’s theorem, it is straightforward to verify that 〈Kφ, φ〉L2[0,1] = ∫ 1
0 ∫ 1 0 (Rφ)(u)R(u,w)du · φ(w)dw
= ∫ 1
0 ∫ 1 0 (Rφ)(u)R(u,w)φ(w)dudw
= ∫ 1
0 (Rφ)(u) (∫ 1 0 R(u,w)φ(w)dw ) du
= ∫ 1
0 (Rφ)2(u)du ≥ 0, φ ∈ L2[0, 1]. (69)
The proof is completed.
Combining (i)—(iv) and applying Hilbert–Schmidt’s expansion theorem, we obtain that L2[0, 1] has an orthonormal basis {φn}n∈N ⋃ {ψξ}ξ∈Ξ, such that
• Kφn = λnφn, λn 6= 0 for n ∈ N , and Kψξ = 0 for ξ ∈ Ξ, where N is a finite or countable set. If N is not finite, we have limn→∞ λn = 0;
• For any ψ ∈ L2[0, 1], we have ψ = ∑ n∈N 〈ψ, φn〉L2[0,1]φn + ∑ ξ∈Ξ 〈ψ,ψξ〉L2[0,1]ψξ, (70)
where the second summation has at most countable non-zero terms, and Kψ = ∑ n∈N λn〈ψ, φn〉L2[0,1]φn. (71)
Here, all the series converge under the norm ‖ · ‖L2[0,1]. Without loss of generality, N = {1, 2, · · · , N0} for N0 ∈ N+ or N0 = +∞ (i.e. N = N+). By (69), we get
0 ≤ 〈Kφn, φn〉L2[0,1] = λn〈φn, φn〉2L2[0,1] = λn, ∀n ∈ N , (72)
i.e. all the eigenvalues of K are non-negative, and λn 6= 0 for n ∈ N implies λn > 0, ∀n ∈ N . In addition, limn→∞ λn = 0 implies that one can index all the eigenvalues in a non-increasing sequence: λ1 ≥ λ2 ≥ · · · ≥ λn ≥ · · · ≥ 0. Then, we can present the POD estimate. Theorem D.1. For any R ∈ L∞([0, 1]2), we have∫ 1
0 ∥∥∥∥∥R(u, v)− N∑ n=1 〈R(u, v), φn(v)〉L2[0,1]φn(v) ∥∥∥∥∥ 2
L2[0,1]
du = N0∑
n=N+1 λn, ∀N ∈ N. (73)
Proof. Combining (69) and (72) gives λn = 〈Kφn, φn〉L2[0,1] = ∫ 1
0 (Rφn)2(u)du, ∀n ∈ N . (74)
Similarly,∫ 1 0 (Rψξ)2(u)du = 〈Kψξ, ψξ〉L2[0,1] = ∑ n∈N λn〈ψξ, φn〉2L2[0,1] = 0, ∀ξ ∈ Ξ, (75)
which gives
(Rψξ)(u) = ∫ 1
0 R(u, v)ψξ(v)dv = 0, a.e. u ∈ [0, 1]. (76)
Notice that Ru(v) := R(u, v) ∈ Cα[0, 1] ⊂ L2[0, 1] for any u ∈ [0, 1]. By (70) and (76), we get
‖Ru‖2L2[0,1] = 〈 ∑ n∈N 〈Ru, φn〉L2[0,1]φn + ∑ ξ∈Ξ 〈Ru, ψξ〉L2[0,1]ψξ,
∑ n∈N 〈Ru, φn〉L2[0,1]φn + ∑ ξ∈Ξ 〈Ru, ψξ〉L2[0,1]ψξ 〉 L2[0,1]
= ∑ n∈N 〈Ru, φn〉2L2[0,1] + ∑ ξ∈Ξ 〈Ru, ψξ〉2L2[0,1]
= ∑ n∈N (Rφn)2(u), a.e. u ∈ [0, 1]. (77)
Hence for any N ∈ N, we have∥∥∥∥∥Ru − N∑ n=1 〈Ru, φn〉L2[0,1]φn ∥∥∥∥∥ 2
L2[0,1]
= ‖Ru‖2L2[0,1] − N∑ n=1 〈Ru, φn〉2L2[0,1] (78) = ∑ n∈N (Rφn)2(u)− N∑ n=1 (Rφn)2(u)
= N0∑
n=N+1 (Rφn)2(u), (79)
where the summation is zero by convention if the subscript is larger than the superscript. This by (74) implies∫ 1
0 ∥∥∥∥∥R(u, v)− N∑ n=1 〈R(u, v), φn(v)〉L2[0,1]φn(v) ∥∥∥∥∥ 2
L2[0,1]
du = N0∑
n=N+1 λn. (80)
Here, the equality holds as a consequence of Beppo Levi’s monotone convergence lemma and Lebesgue’s dominated convergence theorem, and one has ∑ n∈N λn < +∞. In fact, for
N0 = +∞, let Sn = ∑n k=1 λk, we get Sn increasing (since λk ≥ 0 for all k ∈ N ). By (74) and (78), we have
Sn = n∑ k=1 ∫ 1 0 (Rφk)2(u)du = ∫ 1 0 n∑ k=1 〈Ru, φk〉2L2[0,1]du ≤ ∫ 1 0 ‖Ru‖2L2[0,1]du = ‖R‖2L2([0,1]2),
(81) which gives that Sn converges as n→∞. The proof is completed.
Remark D.1. Let ϕn(u) := 〈R(u, v), φn(v)〉L2[0,1], then we have the POD estimate R(u, v) ≈ ∑ n ϕn(u)φn(v), where the error is characterised by the tail sum of eigenvalues of the POD operator.
Recall the POD operator defined in (61). We similarly define K̃ : φ(v) 7→ ∫ 1
0 ∫ 1 0 R̃(u, v)φ(v)dv · R̃(u, v)du, φ(v) ∈ L2[0, 1], (82)
where R̃ is defined as (55), i.e. the approximator constructed in the general approximation theorem before. Obviously, as a polynomial, R̃ ∈ L∞([0, 1]2). Hence, by Proposition D.1, K̃ is also linear, bounded, compact, self-adjoint and non-negative. In addition, according to Theorem D.1, we have the following POD estimate∫ 1
0 ∥∥∥∥∥R̃(u, v)− N∑ n=1 〈 R̃(u, v), φ̃n(v) 〉 L2[0,1] φ̃n(v) ∥∥∥∥∥ 2
L2[0,1]
du = Ñ0∑
n=N+1 λ̃n, (83)
where {λ̃n}Ñ0n=1 are eigenvalues of K̃ satisfying λ̃1 ≥ λ̃2 ≥ · · · ≥ λ̃n ≥ · · · ≥ 0, and {φ̃n} Ñ0 n=1 ⊂ L2[0, 1] are the corresponding orthonormal eigenfunctions, i.e. K̃φ̃n = λ̃nφ̃n, λ̃n > 0 for n ∈ {1, 2, · · · , Ñ0}. Lemma D.1. K̃ is a finite-rank operator. That is, Ñ0 ≤ m̄ = min{mE ,mD} <∞. Proof. Let ϕ̃n(u) := 〈 R̃(u, v), φ̃n(v) 〉 L2[0,1]. We first show that both ϕ̃n(u) and φ̃n(v) are
polynomials. In fact, since R̃(u, v) = ∑m̄ i=1 ∑m̄ j=1 r̃iju
ivj , we have∣∣∣∣ ∂k∂uk R̃(u, v) ∣∣∣∣ = ∣∣∣∣∣∣ m̄∑ i=k m̄∑ j=1 r̃ij i! (i− k)!u i−kvj ∣∣∣∣∣∣ ≤
m̄∑ i=k m̄∑ j=1 |r̃ij | i! (i− k)! , C1(k, m̄), k = 1, 2, · · · , (84)
with the convention that the summation is zero if the subscript is larger than the superscript, i.e. C1(k, m̄) = 0 for any k > m̄. Let C1(m̄) := max{‖R̃‖L∞([0,1]2),max1≤k≤m̄ C1(k, m̄)}, then we have∣∣∣∣ ∂k∂uk R̃(u, v)φ̃n(v)
∣∣∣∣ ≤ C1(m̄)|φ̃n(v)| ∈ L2[0, 1] ⊂ L1[0, 1], k = 0, 1, · · · . (85) According to Lebesgue’s dominated convergence theorem, we get by induction that
dk
duk ϕ̃n(u) =
dk
duk ∫ 1 0 R̃(u, v)φ̃n(v)dv = ∫ 1 0 ∂k ∂uk R̃(u, v)φ̃n(v)dv, k = 0, 1, · · · . (86)
Similarly, we get∣∣∣∣ ∂k∂vk R̃(u, v) ∣∣∣∣ = ∣∣∣∣∣∣ m̄∑ i=1 m̄∑ j=k r̃iju i j! (j − k)!v j−k ∣∣∣∣∣∣ ≤
m̄∑ i=1 m̄∑ j=k |r̃ij | j! (j − k)! , C2(k, m̄), k = 1, 2, · · · , (87)
and C2(k, m̄) = 0 for any k > m̄. Let C2(m̄) := max{‖R̃‖L∞([0,1]2),max1≤k≤m̄ C2(k, m̄)}, then for k = 0, 1, · · · , we have∣∣∣∣ ∂k∂vk R̃(u, v) ∫ 1 0 R̃(u, v)φ̃n(v)dv ∣∣∣∣ ≤ C2(m̄)‖R̃‖L∞([0,1]2)‖φ̃n‖L1[0,1] ≤ C22 (m̄)‖φ̃n‖L2[0,1] = C22 (m̄) ⊂ L1[0, 1]. (88)
According to Lebesgue’s dominated convergence theorem, we get by induction that
dk
dvk φ̃n(v) = 1 λ̃n ∫ 1 0 ∂k ∂vk R̃(u, v) ∫ 1 0 R̃(u, v)φ̃n(v)dvdu, k = 0, 1, · · · . (89)
That is, ϕ̃n, φ̃n ∈ C∞[0, 1] for any n = 1, 2, · · · , Ñ0. Since R̃(u, 0) = R̃(0, v) = 0 for u, v ∈ [0, 1], we get ϕ̃n(0) = φ̃n(0) = 0. Furthermore, we have d k duk ϕ̃n(u) = d k dvk φ̃n(v) = 0 for k > m̄, hence ϕ̃n, φ̃n ∈ Pm̄. Since {φ̃n}Ñ0n=1 ⊂ L2[0, 1] are orthonormal, we must have Ñ0 ≤ m̄ <∞. 4 The proof is completed.

heading: D.2 Approximation rates
text: Perturbation of eigenvalues. First, we need to bound the gap between the eigenvalues {λ̃n}Ñ0n=1 and {λn} N0 n=1 (corresponding to the function R defined in (44)). The following theorem is necessary. Theorem D.2 (Courant–Fischer–Weyl min-max principle; Lax (2002) (Chapter 28, Theorem 4)). Let B be a compact, self-adjoint operator on a Hilbert space H, whose positive eigenvalues are listed in a decreasing order µ1 ≥ µ2 ≥ · · · ≥ µk ≥ · · · > 0. Then
max Sk min x∈Sk, ‖x‖H=1 〈Bx, x〉H = µk, (90)
where Sk ⊂ H is any k-dimensional linear subspace.
Based on it, we have the following lemma to characterise the perturbation of singular values. Lemma D.2. For any R1, R2 ∈ L∞([0, 1]2), we have the estimate∣∣∣∣√λR1k −√λR2k ∣∣∣∣ ≤ ‖R1 −R2‖L2([0,1]2). (91) Proof. According to Theorem D.2 and by (69), we have
λRk = maxSk min φ∈Sk, ‖φ‖L2[0,1]=1 〈KRφ, φ〉L2[0,1] = maxSk min φ∈Sk, ‖φ‖L2[0,1]=1 ‖RRφ‖2L2[0,1]. (92)
Note that RR1 −RR2 = RR1−R2 and by (64), we have√ λR1k = maxSk min φ∈Sk, ‖φ‖L2[0,1]=1 ‖RR1φ‖L2[0,1]
≤ max Sk min φ∈Sk, ‖φ‖L2[0,1]=1
( ‖(RR1 −RR2)φ‖L2[0,1] + ‖RR2φ‖L2[0,1] ) ≤ max Sk min φ∈Sk, ‖φ‖L2[0,1]=1 ( ‖RR1−R2‖+ ‖RR2φ‖L2[0,1]
) ≤ max Sk min φ∈Sk, ‖φ‖L2[0,1]=1 ‖RR2φ‖L2[0,1] + ‖R1 −R2‖L2([0,1]2)
= √ λR2k + ‖R1 −R2‖L2([0,1]2), (93)
4In fact, for any p ∈ Pm̄ with p(0) = 0, we have p ∈ span{v, v2, · · · , vm̄}. Through a standard Schmidt-orthogonalization, we can get ek(v) ∈ Pk, k = 1, 2, · · · , m̄, such that 〈ei, ej〉L2[0,1] = δij , and p ∈ span{e1(v), e2(v), · · · , em̄(v)}. That is, if 〈p, q〉L2[0,1] = 0 for some p, q ∈ Pm̄ with p(0) = 0, q(0) = 0, then their coordinates under the basis {ek}m̄k=1 are orthogonal. Hence, the Ñ0 orthogonal m̄-dimensional coordinates here leads to at most m̄ non-zeros.
and similarly,√ λR2k = maxSk min φ∈Sk, ‖φ‖L2[0,1]=1 ‖RR2φ‖L2[0,1]
≤ max Sk min φ∈Sk, ‖φ‖L2[0,1]=1
( ‖(RR2 −RR1)φ‖L2[0,1] + ‖RR1φ‖L2[0,1] ) ≤ max Sk min φ∈Sk, ‖φ‖L2[0,1]=1 ( ‖RR2−R1‖+ ‖RR1φ‖L2[0,1]
) ≤ max Sk min φ∈Sk, ‖φ‖L2[0,1]=1 ‖RR1φ‖L2[0,1] + ‖R2 −R1‖L2([0,1]2)
= √ λR1k + ‖R1 −R2‖L2([0,1]2), (94)
which completes the proof.
Proofs. Now we are ready to derive the final estimate.
Proof of Theorem 4.3. By Lemma D.2 and (54), we get∣∣∣∣√λk −√λ̃k∣∣∣∣ ≤ ‖R− R̃‖L2([0,1]2) ≤ ‖R− R̃‖L∞([0,1]2) ≤ C(α)γm̄α . (95) Combining (54), (83) and (95) gives that
1 c20 ∥∥∥∥∥ρ(t, s)− N∑ n=1 e− t c0 ϕ̃n(e− t c0 ) · e− s c0 φ̃n(e− s c0 ) ∥∥∥∥∥ L1([0,∞)2)
= ∥∥∥∥∥R(u, v)− N∑ n=1 ϕ̃n(u)φ̃n(v) ∥∥∥∥∥ L1([0,1]2) ≤ ∥∥R(u, v)− R̃(u, v)∥∥ L1([0,1]2) + ∥∥∥∥∥R̃(u, v)− N∑ n=1 ϕ̃n(u)φ̃n(v) ∥∥∥∥∥ L1([0,1]2) ≤ ∥∥R(u, v)− R̃(u, v)∥∥ L∞([0,1]2) + ∥∥∥∥∥R̃(u, v)− N∑ n=1 ϕ̃n(u)φ̃n(v) ∥∥∥∥∥ L2([0,1]2)
≤ C(α)γ m̄α + √√√√ Ñ0∑ n=N+1 λ̃n ≤ C(α)γ m̄α + √√√√ Ñ0∑ n=N+1 |λ̃n − λn|+ Ñ0∑ n=N+1 λn
≤ C(α)γ m̄α + √√√√ Ñ0∑ n=N+1 λn + √√√√ Ñ0∑ n=N+1 ∣∣∣∣√λ̃n −√λn∣∣∣∣ ∣∣∣∣√λ̃n +√λn∣∣∣∣ ≤ C(α)γ
m̄α + √√√√ Ñ0∑ n=N+1 λn + √√√√ Ñ0∑ n=N+1 ∣∣∣∣√λ̃n −√λn∣∣∣∣2 +√2 √√√√ Ñ0∑ n=N+1 √ λn ∣∣∣∣√λ̃n −√λn∣∣∣∣ .C(α)γ  ( 1 + √ Ñ0 −N ) · 1 m̄α + √√√√ Ñ0∑ n=N+1 λn + √√√√ Ñ0∑ n=N+1 √ λn · 1 m̄α/2
 , (96) where . hides universal positive constants.
For the corresponding parameters, recall that ϕ̃n, φ̃n ∈ Pm̄ with ϕ̃n(0) = φ̃n(0) = 0, we can write
ϕ̃n(u) = m̄∑ i=1 Pinu i, φ̃n(v) = m̄∑ j=1 Qnjv j (97)
for any n = 1, 2, · · · , Ñ0. Let c = 1m̄, u = 1m̄, (98) V = −diag(2, 3, · · · , m̄+ 1)/c0, W = −diag(2, 3, · · · , m̄+ 1)/c0, (99) M = PQ with P = [Pin] ∈ Rm̄×N , Q = [Qnj ] ∈ RN×m̄, (100)
then we have
c>eV tMeWsu = c>eV tP ·QeWsu = N∑ n=1 m̄∑ i=1 e− i+1 c0 tPin · m̄∑ j=1 Qnje − j+1c0 s
= N∑ n=1 e− t c0 ϕ̃n(e− t c0 ) · e− s c0 φ̃n(e− s c0 ). (101)
Plugging this into (96) gives
‖H − Ĥ‖ ≤ ‖ρ− ρ̂‖L1([0,∞)2) (102) = ∥∥ρ(t, s)− c>eV tMeWsu∥∥
L1([0,∞)2)
. C(α)γ β2  ( 1 + √ Ñ0 −N ) · 1 m̄α + √√√√ Ñ0∑ n=N+1 λn + √√√√ Ñ0∑ n=N+1 √ λn · 1 m̄α/2  , (103)
and the number of trainable parameters is 2Nm̄. Together with Lemma D.1, the proof is completed.
Extension to multi-dimensional inputs. The above results can be naturally extended to the general case where a d-dimensional input is given (∀d ∈ N+). In fact, let∥∥ρi(t, s)− c>eV tMieWsu∥∥L1([0,∞)2) . , i = 1, 2, · · · , d, (104) for some 0 < 1, then we take
M = (M1,M2, · · · ,Md) ∈ Rm̄×dm̄, (105) W = diag(W,W, · · · ,W ) ∈ Rdm̄×dm̄, U = diag(u, u, · · · , u) ∈ Rdm̄×d, (106)
and have c>eV tMeWsU = c>eV t · (M1,M2, · · · ,Md) · diag(eWs, eWs, · · · , eWs) · diag(u, u, · · · , u)
= (c>eV tM1, c>eV tM2, · · · , c>eV tMd) · diag(eWsu, eWsu, · · · , eWsu) = (c>eV tM1eWsu, c>eV tM2eWsu, · · · , c>eV tMdeWsu). (107)
If the form PQ(= M) is required, it is sufficient to take Mi = PiQi, i = 1, 2, · · · , d, and
P = (P1, P2, · · · , Pd) ∈ Rm̄×dN , Q = diag(Q1, Q2, · · · , Qd) ∈ RdN×dm̄. (108) Therefore, we obtain
d∑ i=1 ∥∥ρi(t, s)− [c>eV tMeWsU]i∥∥L1([0,∞)2) . d , (109) with the number of parameters increased by d-times compared to the corresponding onedimensional setting.

heading: D.3 Case analysis
text: Bounds under different cases. Now we make the comparison between (102) and (58). Recall that {λn}N0n=1 is a positive decreased sequence (with limn→∞ λn = 0 and ∑∞ n=1 λn ≤ ‖R‖2L2([0,1]2) by (81), if N0 = +∞), and Ñ0 ≤ m̄, we have the following cases.
• if Ñ0 = o(m̄), we set N = Ñ0 in (102) and get the same bound as (58), but the number of parameters is only O(Ñ0m̄) = o(m̄2); • if Ñ0 = O(m̄), then (102) implies that∥∥ρ(t, s)− c>eV tMeWsu∥∥ L1([0,∞)2)
. C(α)γ β2 (1 +√m̄−N) · 1m̄α + √√√√ m̄∑ n=N+1 λn + √√√√ m̄∑ n=N+1 √ λn · 1 m̄α/2  . C(α)γ β2  1m̄α− 12 + √√√√ m̄∑ n=N+1 λn + √√√√ m̄∑ n=N+1 √ λn · 1 m̄α/2
 . (110) We are supposed to require that
1 m̄α− 1 2 & max  √√√√ m̄∑ n=N+1 λn, √√√√ m̄∑ n=N+1 √ λn · 1 m̄α/2  ⇔
m̄∑ n=N+1 λn . 1 m̄2α−1 , m̄∑ n=N+1 √ λn . 1 m̄α−1 . (111)
We give the following typical examples to illustrate sufficient conditions to guarantee (111).
– If λn = O(n−r) with r > 2α+ 1 ≥ 3 (α ∈ N+), since m̄∑
n=N+1 n−r ≤ ∫ m̄ N x−rdx = 1 r − 1 ( 1 Nr−1 − 1 m̄r−1 ) , ∀r > 1, (112)
we get by (111) that m̄∑
n=N+1 λn . m̄∑ n=N+1 n−r . 1 Nr−1 . 1 m̄2α−1 ⇔ N & m̄ 2α−1 r−1 , (113)
m̄∑ n=N+1 √ λn . m̄∑ n=N+1 n− r 2 . 1 N r 2−1 . 1 m̄α−1 ⇔ N & m̄ α−1 r 2−1 . (114)
Assume that N ∼ m̄δ with δ ∈ [0, 1), then by (113) and (114), we require δ ≥ max{ 2α−1r−1 , α−1 r 2−1 }, i.e.
r ≥ max {
2α− 1 δ + 1, 2 ( α− 1 δ + 1 )} = 2α− 1 δ + 1. (115)
Meanwhile, the POD-estimate (110) achieves an accuracy scaling like (1/m̄)α− 12 with O(m̄1+δ) parameters, while under the same capacity, the accuracy of (58) scales like (1/m̄)
α(1+δ) 2 . The former beats the latter if α − 12 > α(1+δ) 2 ,
i.e. δ < 1 − 1/α (α ≥ 2). When δ = 1 − 1/α, (115) becomes r ≥ max{ 2α2−1α−1 , 2(α + 1)} = 2α2−1 α−1 . That is to say, r∗ := 2α2−1 α−1 can be viewed as an upper bound of the critical point where the two estimates are comparable. When r > r∗, the POD-estimate outperforms the non-POD-estimate, and this effect gets more notable with r increasing. In fact, when r > r∗, we take N = m̄ 2α−1 r−1 > m̄ α−1 r 2−1 (hence satisfying (113) and (114)), which gives an O((1/m̄)α− 12 ) approximation error with O(m̄1+ 2α−1 r−1 ) parameters for the PODestimate, while O(m̄2− 1α ) trainable parameters are needed to achieve the same accuracy using the non-POD-estimate. Since 2α−1r−1 < 1 − 1 α , we get that the POD-estimate outperforms the non-POD-estimate. When r → +∞, the number of trainable parameters for the POD-estimate is O(m̄), much better than the non-POD-estimate.
Remark D.2. Let K(v, w) := ∫ 1
0 R(u, v)R(u,w)du, v, w ∈ [0, 1], (116)
we get
(Kφ)(w) = ∫ 1
0 K(v, w)φ(v)dv, φ ∈ L2[0, 1]. (117)
Recall that R ∈ Cα([0, 1]2) ⊂ L∞([0, 1]2), we get K ∈ L∞([0, 1]2) ⊂ L2([0, 1]2), and hence K can be also viewed as a Hilbert-Schmidt integral operator with the kernel K, where K is obviously symmetric as K(v, w) = K(w, v), and positive definite since 0 ≤ 〈Kφ, φ〉L2[0,1] = ∫ 1 0 ∫ 1 0 K(v, w)φ(v)φ(w)dvdw, φ ∈ L 2[0, 1] by (69). Since the derivatives of R (up to α-order) are continuous on [0, 1]2 (hence bounded and integrable), we get K ∈ Cα,α([0, 1]2) (α-differentiable for both arguments). According to Chang & Ha (1999) (Theorem 1), we have
∞∑ n=N+2α+1 λn . N −α, ∀N ∈ N+ ⇒ λn = O(n−(α+1)), n→∞. (118)
That is to say, this general setting (only assume smoothness of the kernel) can not guarantee the sufficient condition provided here (λn = O(n−r) with r > 2α
2−1 α−1 ), i.e. the point that the POD-estimate outperforms the non-POD-
estimate requires a faster decay of the eigenvalues. – If λn = O(e−ωn) with ω > 0, we get by (111) that
m̄∑ n=N+1 λn . ∞∑ n=N+1 e−ωn = e −ω(N+1) 1− e−ω . 1 m̄2α−1 ⇔ N & 1 ω ln ( m̄2α−1 1− e−ω ) − 1,
(119) m̄∑
n=N+1
√ λn . ∞∑ n=N+1 e− ω 2 n = e −ω2 (N+1) 1− e−ω2 . 1 m̄α−1 ⇔ N & 2 ω ln ( m̄α−1 1− e−ω2 ) − 1.
(120)
That is to say, for any α ∈ N+, ω > 0, we have N ∼ (2α − 1) ln m̄. This implies an O((1/m̄)α− 12 ) approximation error with O(m̄ ln m̄) parameters for the POD-estimate, while O(m̄2− 1α ) parameters are needed to achieve the same accuracy using the non-POD-estimate.
– If N0 < +∞ (i.e. K is a finite rank operator by (71)), one can just take N = N0 and get by (110) anO((1/m̄)α− 12 ) approximation error. For m̄ ∈ N+ sufficiently large such that m̄ ∼ Nκ0 for some κ 1, the number of trainable parameters for the POD-estimate is O(Nm̄) = O(Nκ+10 ), while O ( N κ(2− 1α ) 0 ) parameters are
needed to achieve the same accuracy using the non-POD-estimate. Obviously, if α ≥ 2, we get κ(2− 1α ) ≥ 3 2κ κ+ 1.
Eigenvalue approximation. One can apply (91) in Lemma D.2 to estimate the eigenvalues {λn}N0n=1, where R1 = R and R2 is taken as some approximator of R, say R̂. Let λ̂ := λR̂, we recall (91) as ∣∣∣∣√λk −√λ̂k∣∣∣∣ ≤ ‖R− R̂‖L2([0,1]2). (121) We provide a naive method here. That is, one can take R̂ as a piece-wise constant approximation of R. Fix any n ∈ N+. Let I0 := {0}, Ii := ( i−1n , i n ] for i = 1, 2, · · · , n, then
{Ii×Ij}ni,j=0 is the uniform partition over [0, 1]2. Let R̂(u, v) := ∑n i,j=0R( in , j n )1Ii×Ij (u, v).
Then for any φ ∈ L2[0, 1], if w ∈ Ik, k = 0, 1, · · · , n, we have (KR̂φ)(w) = ∫ 1
0 ∫ 1 0 R̂(u, v)φ(v)dv · R̂(u,w)du
= n∑ i=0 n∑ j=0 n∑ i′=0 n∑ j′=0 R ( i n , j n ) R ( i′ n , j′ n )∫ 1 0 ∫ 1 0 1Ii×Ij (u, v)φ(v)dv · 1Ii′×Ij′ (u,w)du
= n∑ i=0 n∑ j=0 R ( i n , j n ) R ( i n , k n )∫ Ii ∫ Ij 1Ii×Ij (u, v)φ(v)dv · 1Ii×Ik(u,w)du
= n∑ i=0 n∑ j=0 R ( i n , j n ) R ( i n , k n ) · 1 n ∫ Ij φ(v)dv, (122)
which is a constant only related to k. That is, KR̂φ is a piece-wise constant function, i.e. KR̂φ ∈ span {1Ik} n k=0, or range(KR̂) ⊂ span {1Ik} n k=0. Obviously, {1Ik} n k=0 is an orthogonal set, which implies that the operator KR̂ is of finite rank at most n + 1. Let R̂ij := 1nR( i−1 n , j−1 n ), i, j = 1, 2, · · · , n + 1, {σk} n+1 k=1 be the singular value of R̂ := [R̂ij ] ∈ R(n+1)×(n+1) and V ∈ R(n+1)×(n+1) with columns as the corresponding right singular vectors. Set [e1, e2, · · · , en+1] := [1I0 ,1I1 , · · · ,1In ]V , then by (122), we get for l = 1, 2, · · · , n+ 1,
(KR̂el)(w) = n∑ i=0 n∑ j=0 R ( i n , j n ) R ( i n , k n ) · 1 n2 Vj+1,l
= n+1∑ i=1 n+1∑ j=1 R̂ijR̂i,k+1Vjl = [ R̂>R̂V:,l ] k+1 , (123)
which gives
KR̂el = n∑ k=0 [ R̂>R̂V:,l ] k+1 1Ik = σ2l n∑ k=0 Vk+1,l1Ik = σ2l el, (124)
i.e. the set {σ2k} n+1 k=1 collects the all the eigenvalues of KR̂, which can be obtained by the SVD of R̂. For the error estimate, it is straightforward to have∥∥∥R(u, v)− R̂(u, v)∥∥∥2
L2([0,1]2) = n∑ i=1 n∑ j=1 ∫ i n i−1 n ∫ j n j−1 n ∣∣∣∣R(u, v)−R( in , jn )∣∣∣∣2 dudv ≤ n∑ i=1 n∑ j=1 ∫ i n i−1 n ∫ j n j−1 n max (u,v)∈[0,1]2 ‖∇R(u, v)‖22 · ∥∥∥∥(u− in , v − jn )∥∥∥∥2 2 dudv ≤ n∑ i=1 n∑ j=1 ∫ i n i−1 n ∫ j n j−1 n 2C2(α)γ2 · 2 n2 dudv = 4C 2(α)γ2 n2 , (125)
hence by (121), ∣∣√λk − σk∣∣ ≤ 2C(α)γn for any n ∈ N+.

heading: E Numerical settings
text: According to Lemma B.1, the target input-output temporal relationship has a Riesz representation form. Under the discrete-time regime, we are supposed to set
Ht(x) = T∑ s=1 ρ(t, s)xs, (126)
where T ∈ N+ is the path length.

heading: E.1 Settings of Figure 1
text: For the input, we generate 6 sequences using Gaussian i.i.d. random variables with the path length T = 30. Hence, the output (target) is Ht(x) = ∑30 s=1 ρ(t, s)xs.
The high rank target has the representation ρhigh(t, s) = {
cos(t), t = s, 0, t 6= s, (127)
while the low rank target has the representation
ρlow(t, s) = 99∑ n=0 1 n+ 1 cos(nπt) cos(nπs). (128)

heading: E.2 Settings of Figure 2
text: Consider the target with the following representation
ρ(t, s) = e− t c0 e− s c0 R(e− t c0 , e− s c0 ), (129)
where
R(u, v) = ∞∑ n=1 σnϕn(u)φn(v), (130)
with both {ϕk} and {φk} as orthonormal bases. In this way, we construct a target with singular values {σk}. Under the discrete-time setting, we are supposed to use the following linear, width-m RNN encoder-decoder
Ĥt(x) = τ∑ s=1 c>V tPQW s−1Ux(τ − s), (131)
where P ∈ Rm×N , Q ∈ RN×m with m = mD = mE . Recall that the approximation error is derived as
‖H − Ĥ‖ . ‖R− R̂‖L1([0,1]2) ≤ ‖R− R̂‖L2([0,1]2), (132) where R̂(u, v) := ∑N0 i=1 ϕ̃i(u)φ̃i(v) ∈ P2m with ϕ̃n, φ̃n ∈ Pm. In the numerical experiments, we first construct an R(u, v), and then fit it with the polynomial R̂(u, v) using the method of least squares. The norm ‖R− R̂‖L2([0,1]2) is used to evaluate the approximation error. In the particular example reported in Figure 2, we set m = 128, ϕn(u) = √ 2 sin(nπu) and φn(v) = √ 2 sin(nπv). The singular values are taken as: (a) σn = { n− 1 8 , n ≤ N0
0, n > N0 ; (b) σn = { n−1, n ≤ N0 0, n > N0 ; (c) σn = n−2, with N0 = 2, 4, 6, 8. As an infinite sum, R is constructed under a finite truncation with the first 50 terms.

heading: E.3 Settings of Figure 3
text: We perform experiments on nonlinear targets to show that the insight of low rank approximation also holds in the nonlinear case.
Nonlinear target. We consider the forced Lorenz 96 system (Lorenz, 1996), which is an important example of reduced order modelling for convection dynamics, with applications in weather forecasting. Mathematically, the system has K output variables {yk} and JK hidden variables {zj,k} with k = 1, 2, . . . ,K and j = 1, 2, . . . , J . The parameters K,J control the number of
variables in the system, and can be viewed as a complexity measure. The input {xk} is an external temporal forcing. The system satisfies the following dynamics
dyk dt = −yk−1(yk−2 − yk+1)− yk + xk − 1 J J∑ j=1 zj,k, (133)
dzj,k dt = −zj+1,k(zj+2,k − zj−1,k)− zj,k + yk, (134)
with cyclic indices yk+K = yk, zj,k+K = zj,k and zj+J,k = zj,k. Here, we take {xk} as randomly generated input sequences with the path length 64. We have tested for several cases with different parameters: i) J = 6 with K = 1, 5, 10, 20; ii) K = 5 with J = 5, 15, 25, 100. Note that the forced Lorenz 96 system parameterizes a highly nonlinear functional.
Nonlinear model. We learn the above system using RNN encoder-decoders with nonlinear activations, i.e.
hs = σ(WEhs−1 + UExs + bE), v = σ(Qhτ + b1), gt = σ(WDgt−1 + bD), g0 = σ(Pv + b2), ot = WOgt + bO,
(135)
where σ is the element-wise tanh activation. Let m = 128 be the hidden dimension, N = 1, 2, . . . , 32 be the size of the coding vector v, we have xs, ot, bO ∈ RK , hs, bE , bD, b2 ∈ Rm, WE ,WD ∈ Rm×m, b1 ∈ RN , WO ∈ Rm×K , Q ∈ Rm×N , P ∈ RN×m. Note that we construct the model with a fixed hidden dimension m but different N , thus only sizes of Q,P, b1, b2 are varying, while sizes of other parameters remain unchanged.
Training and initialisation. We denote the model with the coding vector size N as EncDec(N). We utilise the Adam optimiser and train from EncDec(1) to EncDec(32). For EncDec(1), we use a normal random initialisation, and train for 3000 epoches until a stable error. For EncDec(N) with N > 1, we use the parameters trained from EncDec(N−1) as the initialisation. For the parameters Q,P, b1, b2, we pad them to match the size of EncDec(N) with normal distributions as initialisations. It is shown that the low rank approximation phenomena discovered in the linear setting also appears in this nonlinear case.

Review summaries:
This paper provides theoretical studies for why encoder-decoder can be seen as a generalization of RNNs in time-inhomogenous sequence modeling. The authors put an impressive amount of effort into mathematically defining approximation properties of RNN encoder-decoders beginning from a universal approximation result.   The authors first show the universal approximation property of RNN encoder-decoders, and subsequently, they show approximation rates of targets for RNN encoder-decoders. They introduce a notion of temporal products that can characterize the temporal relationships in the input/output pair. 

This paper provides theoretical insight for approximation properties of RNN encoder-decoder architecture in linear setting. More specifically, they study supervised learning problem of temporal modelling where a first RNN encodes a given sequence into a coding vector and a second RNN is responsible to decode said vector to a target output sequence. Linear setting here refers to a linear and continuous-time idealization in eq 6.  Their analysis summarized as following: 1.  Universal approximation: any linear, continuous, and regular temporal relation can be approximated by RNN encoder-decoder up to arbitrary accuracy. 2.  Approximation rate for large size coding vector: beside width of encoder and decoder, this error bound depends on $\alpha$ smoothness of $\mathbf{H}$ and $\beta$ temporal decay rates of the output of constant signal under $\mathbf{H}$ where $\mathbf{H}$ is a sequence of linear, continuous, and regular functionals on inputs. Each $\mathbf{H}$ has a unique two-parameter representation $\rho(t,s)$. Target functionals that are smooth and have fast decaying memory are identified as good target of this approximation. 3.  Approximation rate for small size coding vector: the architecture and the following assumptions in the paper give rise to an intrinsic structure to this type of RNN encoder-decoders called temporal product structure which can be deconstructed into encoder and decoder parts. By relating this structure to rank of temporal relationships they show approximation rate (in small size coding vector) is additionally a function of rank structure of target relationship. This analysis enables us to see the coding vector size as a knob to control number of parameters vs approximation error. 4.  Experiments: they show that the theoretical analysis holds in their experiments.