Long-term Forecasting using Higher Order Tensor RNNs

Introduction

One of the central questions in science is forecasting: given the past history, how well can we predict the future? In many domains with complex multi-variate correlation structures and nonlinear dynamics, forecasting is highly challenging since the system has long-term temporal dependencies and higher-order dynamics. Examples of such systems abound in science and engineering, from biological neural network activity, fluid turbulence, to climate and traffic systems (see, e.g., Figure 1). Since current forecasting systems are unable to faithfully represent the higher-order dynamics, they have limited ability for accurate long-term forecasting.

Therefore, a fundamental challenge is accurately modeling nonlinear dynamics and obtaining stable long-term predictions, given a dataset of realizations of the dynamics. Here, the forecasting problem can be stated as follows: how can we efficiently learn a model that, given only a few initial states, can predict a sequence of future states over a long horizon of $T$ time-steps accurately and reliably?

Common approaches to forecasting include classic linear time series models such as auto-regressive moving average (ARMA), state space models such as hidden Markov model (HMM), and deep neural networks. See a survey on time series forecasting by (Box et al., 2015) and the references therein. A recurrent neural network (RNN), as well as its memory-based extensions such as the LSTM, is a class of models that have achieved state of the art performance on sequence prediction tasks from demand forecasting (Flunkert et al., 2017) to speech recognition (Soltau et al., 2016) and video analysis (LeCun et al., 2015). But most of these methods focus on short-term predictions, and often fail to generalize to nonlinear dynamics and forecast over long time horizons.

In this work, we propose HOT-RNN, a model class that is more expressive and empirically generalizes better than standard RNNs, for the same model capacity. HOT-RNN explicitly models the 1) higher-order dynamics, by incorporating a longer history and higher-order state interactions of previous hidden states; and 2) using tensor trains decomposition that greatly reduces the number of model parameters, while mostly preserving the correlation structure of the full-rank model. We prove that HOT-RNN is exponentially more expressive than standard RNNs for functions that satisfy certain regularity conditions. Our contributions can be summarized as follows:

We propose a novel family of RNNs HOT-RNNs to encode non-Markovian dynamics and higher-order state interactions. To address the memory issue, we propose a tensor-train decomposition that makes learning tractable.

We provide theoretical guarantees for the expressiveness of HOT-RNNs for nonlinear dynamics, and characterize the target dynamics and its HOT-RNN representation. In contrast, no such theoretical results are known for standard recurrent networks.

We show that HOT-RNNs can forecast more accurately for significantly longer time horizons compared to standard RNNs and LSTMs on simulated data and real-world environments with nonlinear dynamics.

Related Work

Time series forecasting is at the core of many dynamics modeling tasks. In statistics, classic work such as the ARMA or ARIMA model (Box et al., 2015) model a stochastic process with assumptions of linear dynamics. In the control and dynamical system community, estimating dynamics models from measurement data is also known as system identification (Ljung, 2001). System identification often requires strong parametric assumptions which are often challenging to find from first principles. Moreover, finding (approximate) solutions of complex nonlinear differential equations demands high computational cost. In our work, we instead take a “mode-free” approach to learn a powerful approximate nonlinear dynamics model.

Recurrent Neural Networks

Using neural networks to model time series data has a long history (Schmidhuber, 2015). Recent developments in deep leaning and RNNs has led to non-linear forecasting models such as deep AutoRegressive (Flunkert et al., 2017), Predictive State Representation (Downey et al., 2017), Deep State Space model (Rangapuram et al., 2018). However, these works usually study short-term forecasting and use RNNs that contain only the most recent state. Our method contrasts with this by explicitly modeling higher-order dyanmics to capture long-term dependencies.

There are several classic work on higher-order RNNs. For example, (Giles et al., 1989) proposes a higher-order RNN to simulate a deterministic finite state machine and recognize regular grammars. The model considers a multiplicative structure of inputs and the most recent hidden state, but is limited to two-way interactions. (Sutskever et al., 2011) also studies tensor RNNs that allow a different hidden-to-hidden weight matrix for every input dimension. Soltani and Jiang (2016) proposes a higher-order RNN that concatenates a sequence of past hidden states, but the underlying state interactions are still linear. Moreover, hierarchical RNNs (Zheng et al., 2016) have been used to model sequential data at multiple temporal resolutions. Our method generalizes all these works to capture higher-order interactions using a hidden-to-hidden tensor.

Tensor methods

Tensor methods have tight connections with neural networks. For example, (Novikov et al., 2015; Stoudenmire and Schwab, 2016) employ tensor-train to compress the weights in neural networks. (Yang et al., 2017) extends this idea to RNNs by reshaping the inputs into a tensor and factorizes the input-hidden weight tensor. However, the purpose of these works is model compression in the input space whereas our method learns the dynamics in the hidden state space. Theoretically, (Cohen et al., 2016) shows convolutional neural networks and hierarchical tensor factorizations are equivalent. (Khrulkov et al., 2017) provides expressiveness analysis for shallow networks using tensor train.

Tensor methods have also been used for sequence modeling. For example, one can apply tensor decomposition as method of moments estimators for latent variable models such as Hidden Markov Models (HMMs) (Anandkumar et al., 2012). Tensor methods have also shown promises in reducing the model dimensionality of multivariate spatiotemporal learning problems (Yu and Liu, 2016), as well as nonlinear system identification (Decuyper et al., 2019). Most recently, Schlag and Schmidhuber (2018) combine tensor product of relational information and recurrent neural networks for natural language reasoning tasks. This work however, to the best of our knowledge, is the first to consider tensor networks within RNNs for sequence learning in environments with nonlinear dynamics.

Higher-Order Tensor RNNs

where $\xi^{i}$ can be an arbitrary (smooth) function of the state ${\mathbf{x}}_{t}$ and its derivatives. Continuous time dynamics are usually described by differential equations while difference equations are employed for discrete time. In continuous time, a classic example is the first-order Lorenz attractor, whose realizations showcase the “butterfly-effect”, a characteristic set of double-spiral orbits. In discrete-time, a non-trivial example is the 1-dimensional Genz dynamics, whose difference equation is:

where $x_{t}$ denotes the system state at time $t$ and $c,w$ are the parameters. Due to the nonlinear nature of the dynamics, such systems exhibit higher-order correlations, long-term dependencies and sensitivity to error propagation, and thus form a challenging setting for forecasting.

Given a sequence of initial states ${\mathbf{x}}_{0}\ldots{\mathbf{x}}_{t}$, the forecasting problem aims to learn a dynamics model $F$ that outputs a sequence of future states ${\mathbf{x}}_{t+1}\ldots{\mathbf{x}}_{T}$.

The system is governed by some unknown dynamics. Hence, accurately approximating the dynamics is critical to learning a good forecasting model and making predictions for long time horizons.

First-order Markovian Models

In deep learning, popular approaches such as recurrent neural networks (RNNs) employ first-order hidden-state models to approximate the dynamics. An RNN with a single cell recursively computes a hidden state ${\mathbf{h}}_{t}$ using the most recent hidden state ${\mathbf{h}}_{t-1}$, generating the output ${\mathbf{y}}_{t}$ from the hidden state ${\mathbf{h}}_{t}$ :

where $f$ is the state transition function, $g$ is the output function and $\{\theta_{f},\theta_{g}\}$ are the corresponding model parameters. A common parametrization scheme for (4) applies a nonlinear activation function such as sigmoid $\sigma$ to a linear map of ${\mathbf{x}}_{t}$ and ${\mathbf{h}}_{t-1}$ as:

where $W^{hx},W^{xh}$ and $W^{hh}$ are the transition weight matrices and ${\mathbf{b}}^{h},{\mathbf{b}}^{x}$ are the biases.

RNNs have many different variations, including LSTMs (Hochreiter and Schmidhuber, 1997) and GRUs (Chung et al., 2014). Although a RNN can approximate any function in theory, its hidden state ${\mathbf{h}}_{t}$ only depends on the previous state ${\mathbf{h}}_{t-1}$ and the input ${\mathbf{x}}_{t}$. Such models do not explicitly capture higher-order dynamics and only implicitly encode long-term dependencies between all historical states ${\mathbf{h}}_{0}\ldots{\mathbf{h}}_{t}$. This limits the representation power of RNNs, especially for forecasting in environments with nonlinear dynamics. Hence, instead of using a wide RNN with many hidden units, we exploit the recurrent cell to design higher-order tensor RNNs that can approximate complex non-linear governing equations.

1 Higher-Order Non-Markovian Models

To effectively learn nonlinear dynamics with higher-order temporal dependency, we propose a family of models that generalizes standard RNNs: higher-order recurrent neural networks, or HOT-RNN. We design HOT-RNNs with two goals in mind: explicitly modeling 1) $L$-order Markov processes with $L$ steps of temporal memory and 2) polynomial interactions between the hidden states ${\mathbf{h}}_{\cdot}$ and ${\mathbf{x}}_{t}$.

First, we consider longer “history”: we keep length $L$ historic states: ${\mathbf{h}}_{t},\cdots,{\mathbf{h}}_{t-L}$:

where $f$ represents the state transition function. In principle, early work (Giles et al., 1989) has shown that with a large enough hidden state size, such recurrent structures are capable of approximating any dynamical system.

which concatenates $L$ previous hidden states. To compute ${\mathbf{h}}_{t}$, we construct a $P$-dimensional transition weight tensor to model degree-$P$ polynomial interactions between hidden states:

where $\alpha$ indices the hidden dimension, $i_{\cdot}$ indices the higher-order terms and $P$ is the total polynomial order. We included the bias unit $1$ in ${\mathbf{s}}$ to account for the first order term, so that ${\mathbf{s}}_{i_{1}}\otimes\cdots\otimes{\mathbf{s}}_{i_{p}}=[1,{\mathbf{h}}_{t},{\mathbf{h}}_{t}{\mathbf{h}}_{t-1},\cdots]$ can include all polynomial expansions of hidden states up to order $P$.

The HOT-RNN with LSTM cell, or “HOT-LSTM”, is defined analogously as:

where $\circ$ denotes the Hadamard product. Note that the bias units are again included.

HOT-RNN is a basic unit that can be incorporated in most of the existing recurrent neural architectures such as convolutional RNN (Xingjian et al., 2015) and hierarchical RNN (Chung et al., 2016). In this work, we use HOT-RNN as a module for sequence-to-sequence (seq2seq) framework (Sutskever et al., 2014) in order to perform long-term forecasting.

As shown in Figure 3, seq2seq models consist of an encoder-decoder pair. The encoder takes an input sequence and learns a hidden representation. The decoder initializes with this hidden representation and generates an output sequence. Both the encoder and the decoder contain multiple layers of higher-order tensor recurrent cells (red). The augmented state ${\mathbf{s}}_{t-1}$ (grey) concatenates the past $L$ hidden states; the HOT-RNN cell takes ${\mathbf{s}}_{t-1}$ and outputs the next hidden state. The encoder encodes the initial states $x_{0},\ldots,x_{t}$ and the decoder predicts $x_{t+1},\ldots,x_{T}$. For each time step $t$, the decoder uses its previous prediction ${\mathbf{y}}_{t}$ as an input.

2 Dimension Reduction with Tensor-Train

Unfortunately, due to the “curse of dimensionality”, the number of parameters in ${\mathcal{W}}_{\alpha}$ with hidden size $H$ grows exponentially as $O(HL^{P})$, which makes the higher-order model prohibitively large to train. To overcome this difficulty, we utilize tensor networks to approximate the weight tensor. Such networks encode a structural decomposition of tensors into low-dimensional components and have been shown to provide the most general approximation to smooth tensors (Orús, 2014). The most commonly used tensor networks are linear tensor networks (LTN), also known as tensor-trains in numerical analysis or matrix-product states in quantum physics (Oseledets, 2011).

with $\alpha_{0}=\alpha_{P}=1$, as depicted in Figure (3). When $r_{0}=r_{P}=1$ the $\{r_{p}\}$ are called the tensor-train rank. With tensor-train decomposition, we can reduce the number of parameters of HOT-RNN from $(HL+1)^{P}$ to $(HL+1)R^{2}P$, with $R=\max_{p}{r_{p}}$ as the upper bound on the tensor-train rank. Thus, a major benefit of tensor-train is that they do not suffer from the curse of dimensionality, which is in sharp contrast to many classical tensor decomposition models, such as the Tucker decomposition.

Approximation Theorem for HOT-RNNs

A significant benefit of using HOT-RNN is that we can theoretically characterize its expressiveness for approximating the underlying dynamics. The main idea is to analyze a class of functions that satisfies certain regularity conditions. For such functions, tensor-train representations preserve the weak differentiability and yield a compact representation.

The following theorem characterizes the representation power of HOT-RNN, viewed as a one-layer hidden neural network, in terms of 1) the regularity of the target function $f$, 2) the dimension of the input space, 3) the tensor train rank and 4) the order of the tensor:

Let the target function $f\in\mathcal{H}^{k}_{\mu}$ be a Hölder continuous function defined on a input domain $\mathcal{I}=I_{1}\times\cdots\times I_{d}$, with bounded derivatives up to order $k$ and finite Fourier magnitude distribution $C_{f}$. A single layer HOT-RNN with $h$ hidden units, $\hat{f}$ can approximate $f$ with approximation error $\epsilon$ at most:

where $C_{f}=\int|\omega|_{1}|\hat{f}(\omega)d\omega|$, $d$ is the dimension of the function, i.e., the size of the state space, $r$ is the tensor-train rank, $p$ is the degree of the higher-order polynomials i.e., the order of the tensor, and $C(k)$ is the coefficient of the spectral expansion of $f$.

Remarks: The result above shows that the number of weights required to approximate the target function $f$ is dictated by its regularity (i.e., its Hölder-continuity order $k$). The expressiveness of HOT-RNN is driven by the selection of the rank $r$ and the polynomial degree $p$; moreover, it improves for functions with increasing regularity. Compared with “first-order” regular RNNs, HOT-RNNs are exponentially more powerful for large rank: if the order $p$ increases, we require fewer hidden units $h$.

Proof sketch: For the full proof, see the Appendix. We design HOT-RNN to approximate the underlying system dynamics. The target function $f({\mathbf{x}})$ represents the state transition function, as in (3.1). We first show that if $f$ preserves weak derivatives, then it has a compact tensor-train representation. Formally, let us assume that $f$ is a Sobolev function: $f\in\mathcal{H}^{k}_{\mu}$, defined on the input space ${\mathcal{I}}=I_{1}\times I_{2}\times\cdots I_{d}$, where each $I_{i}$ is a set of vectors. The space $\mathcal{H}^{k}_{\mu}$ is defined as the functions that have bounded derivatives up to some order $k$ and are $L_{\mu}$-integrable.

where $D^{(i)}f$ is the $i$-th weak derivative of $f$ and $\mu\geq 0$.A weak derivative generalizes the derivative concept for (non)-differentiable functions and is implicitly defined as: e.g. $v\in L^{1}([a,b])$ is a weak derivative of $u\in L^{1}([a,b])$ if for all smooth $\varphi$ with $\varphi(a)=\varphi(b)=0$: $\int_{a}^{b}u(t)\varphi^{\prime}(t)=-\int_{a}^{b}v(t)\varphi(t)$. It is known that any Sobolev function admits a Schmidt decomposition: $f(\cdot)=\sum_{i=0}^{\infty}\sqrt{\lambda_{i}}\gamma(\cdot)_{i}\otimes\phi(\cdot)_{i}$, where $\{\lambda\}$ are the eigenvalues and $\{\gamma\},\{\phi\}$ are the associated eigenfunctions. Hence, for ${\mathbf{x}}\in\mathcal{I}$, we can represent the target function $f({\mathbf{x}})$ as an infinite summation of products of a set of basis functions:

where $\{{\mathcal{A}}^{j}(x_{j})_{\alpha_{j-1}\alpha_{j}}\}$ are basis functions over each input dimension. These basis functions satisfy $\langle{\mathcal{A}}^{j}(\cdot)_{im},{\mathcal{A}}^{j}(\cdot)_{in}\rangle=\delta_{mn}$ for all $j$. If we truncate (11) to a low dimensional subspace (${\mathbf{r}}<\infty$), we obtain a functional approximation of the state transition function $f({\mathbf{x}})$. This approximation is also known as the functional tensor-train (FTT):

In practice, HOT-RNN implements a polynomial expansion of the states using $[{\mathbf{s}},{\mathbf{s}}^{\otimes 2},\cdots,{\mathbf{s}}^{\otimes P}]$, where $P$ is the degree of the polynomial. The final function represented by HOT-RNN is a polynomial approximation of the functional tensor-train function $f_{FTT}$.

Given a target function $f({\mathbf{x}})=f({\mathbf{s}}\otimes\dots\otimes{\mathbf{s}})$, we can express it using FTT and the polynomial expansion of the states ${\mathbf{s}}$. This allows us to characterize HOT-RNN using a family of functions that it can represent. Combined with the classic neural network approximation theory Barron (1993), we can bound the approximation error for HOT-RNN with one hidden layer. The above results applies to the full family of HOT-RNNs, including those using vanilla RNN or LSTM as the recurrent cell.

Variance Bound for HOT-RNN

Assuming the time series is governed by a system whose order of dynamics is at most $P$, represented as a joint probabilistic distribution $P(X_{1},\cdots,X_{P})$. The variance for the HOT-RNN estimator $\hat{C}$ and the true variance $C$ of the population statistics is upper bounded by:

Proof : The tensor-train distribution forms a Gibbs field, thus based on Hammersley-Clifford Theorem, a Gibbs field satisfies global Markov property, therefore, it must be a Conditional Random Field (CRF) as the conditional probability distribution factorizes. According to the duality between tensor network and graphical model Robeva and Seigal (2017), tensor train is the dual graph of CRF. Figure 4 visualizes such dual relationship in the graphical model template. After establishing the relationship between conditional random field and functional tensor-train, we can bound the variance of the HOT-RNN estimator.

We first state the well-celebrated Hammersley-Clifford theorem that gives the sufficient and necessary conditions of which a probability distribution is a Markov Random field.

A graphical model $G$ is a Markov Random Field if and only if the probability distribution $P(X)$ on $G$ is a Gibbs distribution:

where $Z=\sum_{x}\prod_{c\in C_{G}}\psi(X_{c})$ is the normalization constant, $\psi$ are functions defined on maximal cliques, and $C_{G}$ is a set of all maximal cliques in the graph.

When the underlying distribution is a Markov random field and belongs to the exponential family, we can generalize Theorem 3 to kernel functions and obtain the following results.

Altun et al. (2012) Given a Markov random field $X$ with respect to a graphical model $G$, if the sufficient statistics $\Phi(X)=(\Phi(X_{c_{1}}),\cdots,\Phi(X_{c_{i}}))$, then the kernels $k(X,X^{\prime})=\langle\Phi(X),\Phi(X^{\prime})\rangle$ satisfy

where $\Phi(X_{c})$ are the sufficient statistics defined on maximal cliques, $k_{c}(X,X^{\prime})=\langle\Phi(X_{C}),\Phi(X^{\prime}_{C})\rangle$

which is the functional tensor-train model.

Experiments

We conducted exhaustive experiments to examine the behavior of the proposed HOT-RNN model on both synthetic and real-world time series data. The source code is available at https://github.com/yuqirose/tensor_train_RNN.

We validated the accuracy and efficiency of HOT-RNN on the following three datasets.

Genz functions are often used as basis for evaluating high-dimensional function approximation. In particular, they have been used to analyze tensor-train decompositions (Bigoni et al., 2016). There are in total $7$ different Genz functions. (1) $g_{1}(x)=\cos(2\pi w+cx)$, (2) $g_{2}(x)=(c^{-2}+(x+w)^{-2})^{-1}$, (3) $g_{3}(x)=(1+cx)^{-2}$, (4) $e^{-c^{2}\pi(x-w)^{2}}$ (5) $e^{-c^{2}\pi|x-w|}$ (6) $g_{6}(x)=\begin{cases}0\quad x>w\\ e^{cx}\quad else\end{cases}$. For each function, we generated a dataset with $10,000$ samples using (2) with $w=0.5$ and $c=1.0$ and random initial points draw from a range of $[-0.1,0.1]$.

Traffic

We use the traffic data of Los Angeles County highway network collected from California department of transportation http://pems.dot.ca.gov/. The dataset consists of $4$ month speed readings aggregated every $5$ minutes . Due to large number of missing values ($\sim 30\%$) in the raw data, we impute the missing values using the average values of non-missing entries from other sensors at the same time. In total, after processing, the dataset covers $35\,136,$ time-series. We treat each sequence as daily traffic of $288$ time stamps. We up-sample the dataset every $20$ minutes, which results in a dataset of $8\,784$ sequences of daily measurements. We select $15$ sensors as a joint forecasting tasks.

Climate

We use the daily maximum temperature data from the U.S. Historical Climatology Network (USHCN) daily http://cdiac.ornl.gov/ftp/ushcn_daily/. The dataset contains daily measurements for $5$ climate variables for approximately $124$ years. The records were collected across more than $1\,200$ locations and span over $45\,384$ days. We analyze the area in California which contains $54$ stations. We removed the first $10$ years of day, most of which has no observations. We treat the temperature reading per year as one sequence and impute the missing observations using other non-missing entries from other stations across years. We augment the datasets by rotating the sequence every $7$ days, which results in a data set of $5\,928$ sequences.

Figure 5 visualizes the time series from Genz dynamics, traffic and climate systems, respectively. To test the stationarity of the time series, we also perform a Dickey–Fuller test on the real-world traffic and climate data. Dickey–Fuller test is a commonly used statistical test procedure to determine whether a time series is stationary. Its null hypothesis is that a unit root is present in an autoregressive model, hence the time series is not stationary. The test statistics of the traffic and climate data is shown in Table 1, which demonstrate the non-stationarity of the time series.

2 Training Details

We use a seq2seq architecture with HOT-RNN using LSTM as recurrent cells (HOT-LSTM). For all experiments, we use the length-$T$ sequence regression loss $L(y,\hat{y})=\sum_{t=1}^{T}||\hat{y}_{t}-y_{t}||^{2}_{2},$ where $y_{t}=x_{t+1},\hat{y}_{t}$ are the ground truth and model prediction respectively. For all datasets, we used a $80\%-10\%-10\%$ train-validation-test split and train for a maximum of $1e^{4}$ steps. We compute the moving average of the validation loss as an early stopping criteria. We also did not include scheduled sampling Bengio et al. (2015), as we found training with scheduled sampling became highly unstable under a range of annealing schedules.

Hyperparameter Search

All models are trained using RMS-prop with a learning rate decay of $0.8$. We performed an exhaustive search over the hyper-parameters for validation. Table 2 reports the search range of different hyper-parameters used in this work.

Baselines

We compared HOT-RNN against 2 sets of natural baselines: 1st-order RNN (vanilla RNN, LSTM), and matrix RNNs (vanilla MRNN, MLSTM), which use matrix products of multiple hidden states without factorization (Soltani and Jiang, 2016). We observed that HOT-RNN with RNN cells outperforms vanilla RNN and MRNN, but using LSTM cells performs best in all experiments. We also evaluated the classic ARIMA time series model with AR lags of $1\sim 5$, and MA lags of $1\sim 3$. We observed that it consistently performs $\sim 5\%$ worse than LSTM.

3 Long-term Forecasting Accuracy

We evaluate the long-term forecasting accuracy of the proposed method and the baselines. For traffic, we forecast up to $18$ hours ahead with $5$ hours as inputs. For climate, we forecast up to $300$ days ahead given $60$ days of observations. For Genz dynamics, we forecast for $80$ steps given $5$ initial steps. We report the forecasting results averaged over $3$ runs.

Figure 6 shows the test prediction error (in RMSE) for varying forecasting horizons for different datasets. We can see that HOT-LSTM notably outperforms all baselines on all datasets in this setting. In particular, HOT-LSTM is more robust to long-term error propagation. We observe two salient benefits of using HOT-RNNs over the unfactorized models. First, MRNN and MLSTM can suffer from overfitting as the number of weights increases. Second, on traffic, unfactorized models also show considerable instability in their long-term predictions. These results suggest that HOT-RNNs learn more stable representations that generalize better for long-term horizons. To compare the performance on high-dimensional time series, we also evaluated on the unsupervised video prediction task for Moving MNIST. We forecast $20$ and $40$ frames ahead given $10$ initial frames. The per-pixel forecasting RMSE results are shown in Table 3. We observe a small gain ($\sim 2-5\%$) of HOT-LSTM over the baselines. This is likely due to the fact that the underlying circular dynamics are still pretty simple. Moreover, the high-dimensional inputs have spatial structure that are hard to learn by RNNs alone (note we do not use convolutional features). We expect HOT-LSTM to improve over baselines even more with more complicated dynamics and using convolutional features.

4 Visualization of Predictions

To get intuition for the learned models, we visualize predictions from the best performing HOT-LSTM and baselines. Figure 7 shows the predictions for the Genz function “corner-peak” as the state-transition function from three realizations of Genz dynamics. We can see that HOT-LSTM can almost perfectly recover the original function, while LSTM and MLSTM only correctly predict the mean. These baselines cannot capture the dynamics fully, often predicting an incorrect range and phase for the dynamics.

Figure 8 shows predictions for the traffic and climate datasets. This work uses deterministic models, hence the predictions correspond to the trend. We can see that the HOT-LSTM aligns significantly better with ground truth in long-term forecasting. As the ground truth time series is highly nonlinear and noisy, LSTM often deviates from the general trend. While both MLSTM and HOT-LSTM can correctly learn the trend, HOT-LSTM captures more detailed curvatures due to higher-order structure.

5 Model Capacity

The number of parameters for HOT-RNN is $\mathcal{O}(HL+1)R^{2}P$ with hidden size $H$, lag $L$, rank $R$ and order $P$. This gives us more flexibility to decide the model capacity. Fewer parameters may have limited representation power, while more parameters would cause overfitting.

Note that the memory complexity only grows quadratically with the rank $R$ while Theorem 1 shows the expressiveness of HOT-RNN improves exponentially. In practice, we used cross-validation to select the values for these hyper-parameters. The best models on real-world climate and traffic data are listed in Table 4. We can see that the number of parameters of HOT-LSTM model is comparable with that of MLSTM and LSTM.

6 Speed Performance Trade-off

We now investigate potential trade-offs between accuracy and computation. Figure 9 displays the validation loss with respect to the number of steps, for the best performing models on long-term forecasting. We see that HOT-RNNs converge faster than other models, and achieve lower validation-loss. This suggests that HOT-RNN has a more efficient representation of the nonlinear dynamics, and can learn much faster as a result.

7 Sensitivity Analysis

The HOT-LSTM model has several hyperparameters, such as tensor-train rank and lag $L$. We study the sensitivity of HOT-LSTM to these hyperparameters; Table 5 shows the results. In the top row, we report the prediction RMSE for the largest forecasting horizon w.r.t tensor ranks for all the datasets with lag $3$. When the rank is too low, the model does not have enough capacity to capture non-linear dynamics. When the rank is too high, the model starts to overfit. In the bottom row, we report the effect of changing lag $L$. For each setting, the best $r$ is determined by cross-validation. Note that the best lag $L$ also varies for different forecasting horizons.

8 Chaotic Nonlinear Dynamics

Chaotic dynamics such as Lorenz attractor is notoriously different to lean in non-linear dynamics. In such systems, the dynamics are highly sensitive to perturbations in the input state: two close points can move exponentially far apart under the dynamics. We also evaluated tensor-train neural networks on long-term forecasting for Lorenz attractor and report the results.

The Lorenz attractor system describes a two-dimensional flow of fluids:

This system has chaotic solutions (for certain parameter values) that revolve around the so-called Lorenz attractor. We simulated $10\,000$ trajectories with the discretized time interval length $0.01$. We sample from each trajectory every $10$ units in Euclidean distance.

As shown in Figure 10, the blue trajectory represents the discretized dynamics and red circles are sampled observations. The dynamics is generated using $\sigma=10$ $\rho=28$, $\beta=2.667$. The initial condition of each trajectory is sampled uniformly random from the interval of $[-0.1,0.1]$.

Figure 11 shows $45$ steps ahead predictions for all models. HORNN is the full tensor HOT-RNN using vanilla RNN unit without the tensor-train decomposition. We can see all the tensor models perform better than vanilla RNN or MRNN. HOT-RNN shows slight improvement at the beginning state.

We have also evaluated HOT-RNN on long-term forecasting for chaotic dynamics, such as the Lorenz dynamics. Such dynamics are highly sensitive to input perturbations: two close points can move exponentially far apart under the dynamics. This makes long-term forecasting highly challenging, as small errors can lead to catastrophic long-term errors. Figure 12 shows that HOT-RNN can predict up to $T=40$ steps into the future, but diverges quickly beyond that. We have found no state-of-the-art prediction model is stable beyond $40$ time step in this setting.

Discussion

In this paper, We studied long-term forecasting under nonlinear dynamics. We proposed a novel class of RNNs – HOT-RNN that directly learns the nonlinear dynamics using higher-order structures. We provided the first approximation guarantees for its representation power. We demonstrated the benefits of HOT-RNN to forecast accurately for significantly longer time horizon in both synthetic and real-world multivariate time series data.

In terms of future work, forecasting chaotic dynamics, still presents a significant challenge to any sequential prediction model. Hence, it would be worthwhile to study how to learn robust models for chaotic dynamics. For other sequence modeling tasks, such as language, there does not (or is not known to) exist a succinct analytical description of the data-generating process. It would also be interesting to go beyond forecasting and further investigate the effectiveness of HOT-RNNs in such domains as well.

We would like to acknowledge support for this project from the National Science Foundation (NSF grant IIS-9988642) and the Multidisciplinary Research Program of the Department of Defense (MURI N00014-00-1-0637).

References

Appendix A.

We provide theoretical guarantees for the proposed HOT-RNN model by analyzing a class of functions that satisfy some regularity condition. For such functions, tensor-train decomposition preserve weak differentiability and yield a compact representation. We combine this property with neural network theory to bound the approximation error for HOT-RNN with one hidden layer, in terms of: 1) the regularity of the target function $f$, 2) the dimension of the input, and 3) the tensor train rank.

In the context of HOT-RNN, the target function $f({\mathbf{x}})$ with ${\mathbf{x}}={\mathbf{s}}\otimes\ldots\otimes{\mathbf{s}}$, is the system dynamics that describes state transitions. Let us assume that $f({\mathbf{x}})$ is a Sobolev function: $f\in\mathcal{H}^{k}_{\mu}$, defined on the input space ${\mathcal{I}}=I_{1}\times I_{2}\times\cdots I_{d}$, where each $I_{i}$ is a set of vectors. The space $\mathcal{H}^{k}_{\mu}$ is defined as the set of functions that have bounded derivatives up to some order $k$ and are $L_{\mu}$-integrable:

where $D^{(i)}f$ is the $i$-th weak derivative of $f$ and $\mu\geq 0$.A weak derivative generalizes the derivative concept for (non)-differentiable functions and is implicitly defined as: e.g. $v\in L^{1}([a,b])$ is a weak derivative of $u\in L^{1}([a,b])$ if for all smooth $\varphi$ with $\varphi(a)=\varphi(b)=0$: $\int_{a}^{b}u(t)\varphi^{\prime}(t)=-\int_{a}^{b}v(t)\varphi(t)$.

Any Sobolev function admits a Schmidt decomposition: $f(\cdot)=\sum_{i=0}^{\infty}\sqrt{\lambda_{i}}\gamma(\cdot)_{i}\otimes\phi(\cdot)_{i}$, where $\{\lambda\}$ are the eigenvalues and $\{\gamma\},\{\phi\}$ are the associated eigenfunctions. Applying the Schmidt decomposition along $x_{1}$, we have

We can apply similar Schmidt decomposition along $x_{2}$

Recursively performing such operation, and let $\gamma(x_{d})_{\alpha_{d-1},\alpha_{d}}=\sqrt{\lambda_{\alpha_{d-1}}}\phi(x_{d})_{\alpha_{d-1}}$ and ${\mathcal{A}}^{j}(x_{j})_{\alpha_{j-1},\alpha_{j}}=\gamma(x_{j})_{\alpha_{j-1},\alpha_{j}}$, the target function $f\in\mathcal{H}^{k}_{\mu}$ can be decomposed as:

where $\{{\mathcal{A}}^{j}(\cdot)_{\alpha_{j-1}\alpha_{j}}\}$ are basis functions, satisfying $\langle{\mathcal{A}}^{j}(\cdot)_{im},{\mathcal{A}}^{j}(\cdot)_{in}\rangle=\delta_{mn}$. We can truncate Eqn 17 to a low dimensional subspace (${\mathbf{r}}<\infty$), and obtain the functional tensor-train (FTT) approximation of the target function $f$:

FTT approximation in Eqn 17 projects the target function to a subspace with finite basis. And the approximation error can be bounded using the following Lemma:

Let $f\in\mathcal{H}^{k}_{\mu}$ for $k>0$. Let $P$ be the approximating polynomial with degree $p$, Then

So far, we have obtained the tensor-train approximation error with the regularity of the target function $f$. Next we will connect the tensor-train approximation and the approximation error of neural networks with one layer hidden units. Given a neural network with one hidden layer and sigmoid activation function, following Lemma describes the classic result of describes the error between a target function $f$ and the single hidden-layer neural network that approximates it best:

Given a function $f$ with finite Fourier magnitude distribution $C_{f}$, there exists a neural network with $n$ hidden units $f_{n}$, such that

where $C_{f}=\int|\omega|_{1}|\hat{f}(\omega)|d\omega$ with Fourier representation $f(x)=\int e^{i\omega x}\hat{f}(\omega)d\omega$.

We can now generalize Barron’s approximation lemma 7 to HOT-RNN. The target time series is $f({\mathbf{x}})=f({\mathbf{s}}\otimes\dots\otimes{\mathbf{s}})$. We can express the function using FTT, followed by the polynomial expansion of the states concatenation $P_{TT}$. The approximation error of HOT-RNN, viewed as one layer hidden

Where $p$ is the order of tensor and $r$ is the tensor-train rank. As the rank of the tensor-train and the polynomial order increase, the required size of the hidden units become smaller, up to a constant that depends on the regularity of the underlying dynamics $f$.

.2 Additional Experiments

Genz functions are often used as basis for evaluating high-dimensional function approximation. Figure 14 visualizes different Genz functions, realizations of dynamics and predictions from HOT-LSTM and baselines. We can see for “oscillatory”, “product peak” and “Gaussian ”, HOT-LSTM can better capture the complex dynamics, leading to more accurate predictions.

Moving MNIST

Moving MNIST Srivastava et al. (2015) generates around $50,000$ video sequences of length $100$ on the fly. The video is generated by moving the digits in the MNIST image dataset along a given trajectory within a canvas of size $48\times 48$. The trajectory reflects the dynamics of the movement. In this experiment, we used $cos$ and $sin$ velocity.