HiPPO: Recurrent Memory with Optimal Polynomial Projections

Introduction

Modeling and learning from sequential data is a fundamental problem in modern machine learning, underlying tasks such as language modeling, speech recognition, video processing, and reinforcement learning. A core aspect of modeling long-term and complex temporal dependencies is memory, or storing and incorporating information from previous time steps. The challenge is learning a representation of the entire cumulative history using bounded storage, which must be updated online as more data is received.

One established approach is to model a state that evolves over time as it incorporates more information. The deep learning instantiation of this approach is the recurrent neural network (RNN), which is known to suffer from a limited memory horizon (e.g., the “vanishing gradients” problem). Although various heuristics have been proposed to overcome this, such as gates in the successful LSTM and GRU , or higher-order frequencies in the recent Fourier Recurrent Unit and Legendre Memory Unit (LMU) , a unified understanding of memory remains a challenge. Furthermore, existing methods generally require priors on the sequence length or timescale and are ineffective outside this range ; this can be problematic in settings with distribution shift (e.g. arising from different instrument sampling rates in medical data ). Finally, many of them lack theoretical guarantees on how well they capture long-term dependencies, such as gradient bounds. To design a better memory representation, we would ideally (i) have a unified view of these existing methods, (ii) be able to address dependencies of any length without priors on the timescale, and (iii) have a rigorous theoretical understanding of their memory mechanism.

By posing a formal optimization problem underlying recurrent sequence models, the HiPPO framework (Section 2) generalizes and explains previous methods, unlocks new methods appropriate for sequential data at different timescales, and comes with several theoretical guarantees. (i) For example, with a short derivation we exactly recover as a special case the LMU (Section 2.3), which proposes an update rule that projects onto fixed-length sliding windows through time.The LMU was originally motivated by spiking neural networks in modeling biological nervous systems; its derivation is not self-contained but a sketch can be pieced together from . HiPPO also sheds new light on classic techniques such as the gating mechanism of LSTMs and GRUs, which arise in one extreme using only low-order degrees in the approximation (Section 2.5). (ii) By choosing more suitable measures, HiPPO yields a novel mechanism (Scaled Legendre, or LegS) that always takes into account the function’s full history instead of a sliding window. This flexibility removes the need for hyperparameters or priors on the sequence length, allowing LegS to generalize to different input timescales. (iii) The connections to dynamical systems and approximation theory allows us to show several theoretical benefits of HiPPO-LegS: invariance to input timescale, asymptotically more efficient updates, and bounds on gradient flow and approximation error (Section 3).

We integrate the HiPPO memory mechanisms into RNNs, and empirically show that they outperform baselines on standard tasks used to benchmark long-term dependencies. On the permuted MNIST dataset, our hyperparameter-free HiPPO-LegS method achieves a new state-of-the-art accuracy of 98.3%, beating the previous RNN SoTA by over 1 point and even outperforming models with global context such as transformers (Section 4.1). Next, we demonstrate the timescale robustness of HiPPO-LegS on a novel trajectory classification task, where it is able to generalize to unseen timescales and handle missing data whereas RNN and neural ODE baselines fail (Section 4.2). Finally, we validate HiPPO’s theory, including computational efficiency and scalability, allowing fast and accurate online function reconstruction over millions of time steps (Section 4.3). Code for reproducing our experiments is available at https://github.com/HazyResearch/hippo-code.

The HiPPO Framework: High-order Polynomial Projection Operators

We motivate the problem of online function approximation with projections as an approach to learning memory representations (Section 2.1). Section 2.2 describes the general HiPPO framework to derive memory updates, including a precise definition of the technical problem we introduce, and an overview of our approach to solving it. Section 2.3 instantiates the framework to recover the LMU and yield new memory updates (e.g. HiPPO-LagT), demonstrating the generality of the HiPPO framework. Section 2.4 discusses how to convert the main continuous-time results into practical discrete versions. Finally in Section 2.5 we show how gating in RNNs is an instance of HiPPO memory.

Any $N$-dimensional subspace $\mathcal{G}$ of this function space is a suitable candidate for the approximation. The parameter $N$ corresponds to the order of the approximation, or the size of the compression; the projected history can be represented by the $N$ coefficients of its expansion in any basis of $\mathcal{G}$. For the remainder of this paper, we use the polynomials as a natural basis, so that $\mathcal{G}$ is the set of polynomials of degree less than $N$. We note that the polynomial basis is very general; for example, the Fourier basis $\sin(nx),\cos(nx)$ can be seen as polynomials on the unit circle $(e^{2\pi ix})^{n}$ (cf. Section D.4). In Appendix C, we additionally formalize a more general framework that allows different bases other than polynomials by tilting the measure with another function.

Since we care about approximating $f_{\leq t}$ for every time $t$, we also let the measure vary through time. For every $t$, let $\mu^{(t)}$ be a measure supported on $(-\infty,t]$ (since $f_{\leq t}$ is only defined up to time $t$). Overall, we seek some $g^{(t)}\in\mathcal{G}$ that minimizes $\|f_{\leq t}-g^{(t)}\|_{L_{2}(\mu^{(t)})}$. Intuitively, the measure $\mu$ controls the importance of various parts of the input domain, and the basis defines the allowable approximations. The challenge is how to solve the optimization problem in closed form given $\mu^{(t)}$, and how these coefficients can be maintained online as $t\to\infty$.

2 General HiPPO framework

We provide a brief overview of the main ideas behind solving this problem, which provides a surprisingly simple and general strategy for many measure families $\mu^{(t)}$. This framework builds upon a rich history of the well-studied orthogonal polynomials and related transforms in the signal processing literature. Our formal abstraction (Definition 1) departs from prior work on sliding transforms in several ways, which we discuss in detail in Section A.1. For example, our concept of the time-varying measure allows choosing $\mu^{(t)}$ more appropriately, which will lead to solutions with qualitatively different behavior. Appendix C contains the full details and formalisms of our framework.

As mentioned, the approximated function can be represented by the $N$ coefficients of its expansion in any basis; the first key step is to choose a suitable basis $\{g_{n}\}_{n<N}$ of $\mathcal{G}$. Leveraging classic techniques from approximation theory, a natural basis is the set of orthogonal polynomials for the measure $\mu^{(t)}$, which forms an orthogonal basis of the subspace. Then the coefficients of the optimal basis expansion are simply $c^{(t)}_{n}:=\langle f_{\leq t},g_{n}\rangle_{\mu^{(t)}}$.

$\operatorname{proj}_{t}$ takes the function $f$ restricted up to time $t$, $f_{\leq t}\coloneqq f(x)\mid_{x\leq t}$, and maps it to a polynomial $g^{(t)}\in\mathcal{G}$, that minimizes the approximation error $\|f_{\leq t}-g^{(t)}\|_{L_{2}(\mu^{(t)})}$.

Figure 1 illustrates the overall framework when we use uniform measures. Next, we give our main results showing $\operatorname{hippo}$ for several concrete instantiations of the framework.

3 High Order Projection: Measure Families and HiPPO ODEs

Our main theoretical results are instantiations of HiPPO for various measure families $\mu^{(t)}$. We provide two examples of natural sliding window measures and the corresponding projection operators. The unified perspective on memory mechanisms allows us to derive these closed-form solutions with the same strategy, provided in Appendices D.1,D.2. The first explains the core Legendre Memory Unit (LMU) update in a principled way and characterizes its limitations, while the other is novel, demonstrating the generality of the HiPPO framework. Appendix D contrasts the tradeoffs of these measures (Fig. 5), contains proofs of their derivations, and derives additional HiPPO formulas for other bases such as Fourier (recovering the Fourier Recurrent Unit ) and Chebyshev.

The translated Legendre (LegT) measures assign uniform weight to the most recent history $[t-\theta,t]$. There is a hyperparameter $\theta$ representing the length of the sliding window, or the length of history that is being summarized. The translated Laguerre (LagT) measures instead use the exponentially decaying measure, assigning more importance to recent history.

Equation (1) proves the LMU update [71, equation (1)]. Additionally, our derivation (Appendix D.1) shows that outside of the projections, there is another source of approximation. This sliding window update rule requires access to $f(t-\theta)$, which is no longer available; it instead assumes that the current coefficients $c(t)$ are an accurate enough model of the function $f(x)_{x\leq t}$ that $f(t-\theta)$ can be recovered.

4 HiPPO recurrences: from Continuous to Discrete Time with ODE Discretization

Since actual data is inherently discrete (e.g. sequences and time series), we discuss how the HiPPO projection operators can be discretized using standard techniques, so that the continuous-time HiPPO ODEs become discrete-time linear recurrences.

The basic method of discretizating an ODE $\frac{d}{dt}c(t)=u(t,c(t),f(t))$ chooses a step size $\Delta t$ and performs the discrete updates $c(t+\Delta t)=c(t)+\Delta t\cdot u(t,c(t),f(t))$.This is known as the Euler method, used for illustration here; our experiments use the more numerically stable Bilinear and ZOH methods. Section B.3 provides a self-contained overview of our full discretization framework. In general, this process is sensitive to the discretization step size hyperparameter $\Delta t$.

Finally, we note that this provides a way to seamlessly handle timestamped data, even with missing values: the difference between timestamps indicates the (adaptive) $\Delta t$ to use in discretization . Section B.3 contains a full discussion of discretization.

5 Low Order Projection: Memory Mechanisms of Gated RNNs

As a special case, we consider what happens if we do not incorporate higher-order polynomials in the projection problem. Specifically, if $N=1$, then the discretized version of HiPPO-LagT (2) becomes $c(t+\Delta t)=c(t)+\Delta t(-Ac(t)+Bf(t))=(1-\Delta t)c(t)+\Delta tf(t)$, since $A=B=1$. If the inputs $f(t)$ can depend on the hidden state $c(t)$ and the discretization step size $\Delta t$ is chosen adaptively (as a function of input $f(t)$ and state $c(t)$), as in RNNs, then this becomes exactly a gated RNN. For instance, by stacking multiple units in parallel and choosing a specific update function, we obtain the GRU update cell as a special case.The LSTM cell update is similar, with a parameterization known as “tied” gates . In contrast to HiPPO which uses one hidden feature and projects it onto high order polynomials, these models use many hidden features but only project them with degree 1. This view sheds light on these classic techniques by showing how they can be derived from first principles.

HiPPO-LegS: Scaled Measures for Timescale Robustness

Exposing the tight connection between online function approximation and memory allows us to produce memory mechanisms with better theoretical properties, simply by choosing the measure appropriately. Although sliding windows are common in signal processing (Section A.1), a more intuitive approach for memory should scale the window over time to avoid forgetting.

Simply by specifying the desired measure, specializing the HiPPO framework (Sections 2.2, 2.4) yields a new memory mechanism (proof in Section D.3).

The continuous- (3) and discrete- (4) time dynamics for HiPPO-LegS are:

We show that HiPPO-LegS enjoys favorable theoretical properties: it is invariant to input timescale, is fast to compute, and has bounded gradients and approximation error. All proofs are in Appendix E.

As the window size of LegS is adaptive, projection onto this measure is intuitively robust to timescales. Formally, the HiPPO-LegS operator is timescale-equivariant: dilating the input $f$ does not change the approximation coefficients.

For any scalar $\alpha>0$, if $h(t)=f(\alpha t)$, then $\operatorname{hippo}(h)(t)=\operatorname{hippo}(f)(\alpha t)$. In other words, if $\gamma:t\mapsto\alpha t$ is any dilation function, then $\operatorname{hippo}(f\circ\gamma)=\operatorname{hippo}(f)\circ\gamma$.

Informally, this is reflected by HiPPO-LegS having no timescale hyperparameters; in particular, the discrete recurrence (4) is invariant to the discretization step size.(4) uses the Euler method for illustration; HiPPO-LegS is invariant to other discretizations (Section B.3). By contrast, LegT has a hyperparameter $\theta$ for the window size, and both LegT and LagT have a step size hyperparameter $\Delta t$ in the discrete time case. This hyperparameter is important in practice; Section 2.5 showed that $\Delta t$ relates to the gates of RNNs, which are known to be sensitive to their parameterization . We empirically demonstrate the benefits of timescale robustness in Section 4.2.

In order to compute a single step of the discrete HiPPO update, the main operation is multiplication by the (discretized) square matrix $A$. More general discretization specifically requires fast multiplication for any matrix of the form $I+\Delta t\cdot A$ and $(I-\Delta t\cdot A)^{-1}$ for arbitrary step sizes $\Delta t$. Although this is generically a $O(N^{2})$ operation, LegS operators use a fixed $A$ matrix with special structure that turns out to have fast multiplication algorithms for any discretization.It is known that large families of structured matrices related to orthogonal polynomials are efficient .

Under any generalized bilinear transform discretization (cf. Section B.3), each step of the HiPPO-LegS recurrence in equation (4) can be computed in $O(N)$ operations.

Section 4.3 validates the efficiency of HiPPO layers in practice, where unrolling the discretized versions of Theorem 2 is 10x faster than standard matrix multiplication as done in standard RNNs.

Much effort has been spent to alleviate the vanishing gradient problem in RNNs , where backpropagation-based learning is hindered by gradient magnitudes decaying exponentially in time. As LegS is designed for memory, it avoids the vanishing gradient issue.

For any times $t_{0}<t_{1}$, the gradient norm of HiPPO-LegS operator for the output at time $t_{1}$ with respect to input at time $t_{0}$ is $\left\|\frac{\partial c(t_{1})}{\partial f(t_{0})}\right\|=\Theta\left(1/t_{1}\right)$.

The error rate of LegS decreases with the smoothness of the input.

Empirical Validation

The HiPPO dynamics are simple recurrences that can be easily incorporated into various models. We validate three claims that suggest that when incorporated into a simple RNN, these methods–especially HiPPO-LegS–yield a recurrent architecture with improved memory capability. In Section 4.1, the HiPPO-LegS RNN outperforms other RNN approaches in benchmark long-term dependency tasks for RNNs. Section 4.2 shows that HiPPO-LegS RNN is much more robust to timescale shifts compared to other RNN and neural ODE models. Section 4.3 validates the distinct theoretical advantages of the HiPPO-LegS memory mechanism, allowing fast and accurate online function reconstruction over millions of time steps. Experiment details and additional results are described in Appendix F.

We first describe briefly how HiPPO memory updates can be incorporated into a simple neural network architecture, yielding a simple RNN model reminiscent of the classic LSTM. Given inputs $x_{t}$ or features thereof $f_{t}=u(x_{t})$ in any model, the HiPPO framework can be used to memorize the history of features $f_{t}$. Thus, given any RNN update function $h_{t}=\tau(h_{t-1},x_{t})$, we simply replace $h_{t-1}$ with a projected version of the entire history of $h$, as described in Figure 2. The output of each cell is $h_{t}$, which can be passed through any downstream module (e.g. a classification head trained with cross-entropy) to produce predictions.

We map the vector $h_{t-1}$ to 1D with a learned encoding before passing to $\operatorname{hippo}$ (full architecture in App. F.1).

1 Long-range Memory Benchmark Tasks

We consider all of the HiPPO methods (LegT, LagT, and LegS). As we show that many different update dynamics seem to lead to LTI systems that give sensible results (Section 2.3), we additionally consider the Rand baseline that uses random $A$ and $B$ matrices (normalized appropriately) in its updates, to confirm that the precise derived dynamics are important. LegT additionally considers an additional hyperparameter $\theta$, which should be set to the timescale of the data if known a priori; to show the effect of the timescale, we set it to the ideal value as well as values that are too large and small. The MGU is a minimal gated architecture, equivalent to a GRU without the reset gate. The HiPPO architecture we use is simply the MGU with an additional $\operatorname{hippo}$ intermediate layer.

We also compare to several RNN baselines designed for long-term dependencies, including the LSTM , GRU , expRNN , and LMU .In our experiments, LMU refers to the architecture in while LegT uses the one described in Fig. 2.

All methods have the same hidden size in our experiments. In particular, for simplicity and to reduce hyperparameters, HiPPO variants tie the memory size $N$ to the hidden state dimension $d$, so that all methods and baselines have a comparable number of hidden units and parameters. A more detailed comparison of model architectures is in Section F.1.

The permuted MNIST (pMNIST) task feeds inputs to a model pixel-by-pixel in the order of a fixed permutation. The model must process the entire image sequentially – with non-local structure – before outputting a classification label, requiring learning long-term dependencies.

Table 1 shows the validation accuracy on the pMNIST task for the instantiations of our framework and baselines. We highlight that LegS has the best performance of all models. While LegT is close at the optimal hyperparameter $\theta$, its performance can fall off drastically for a mis-specified window length. LagT also performs well at its best hyperparameter $\Delta t$.

Table 1 also compares test accuracy of our methods against reported results from the literature, where the LMU was the state-of-the-art for recurrent models. In addition to RNN-based baselines, other sequence models have been evaluated on this dataset, despite being against the spirit of the task because they have global receptive field instead of being strictly sequential. With a test accuracy of 98.3%, HiPPO-LegS sets a true state-of-the-art accuracy on the permuted MNIST dataset.

This standard RNN task directly tests memorization, where models must regurgitate a sequence of tokens seen at the beginning of the sequence. It is well-known that standard models such as LSTMs struggle to solve this task. Appendix F shows the loss for the Copying task with length $L=200$. Our proposed update LegS solves the task almost perfectly, while LegT is very sensitive to the window length hyperparameter. As expected, most baselines make little progress.

2 Timescale Robustness of HiPPO-LegS

Sequence models generally benefit from priors on the timescale, which take the form of additional hyperparameters in standard models. Examples include the “forget bias” of LSTMs which needs to be modified to address long-term dependencies , or the discretization step size $\Delta t$ of HiPPO-Lag and HiPPO-LegT (Section 2.4). The experiments in Section 4.1 confirm their importance. Fig. 7 (Appendix) and Table 1 ablate these hyperparameters, showing that for example the sliding window length $\theta$ must be set correctly for LegT. Additional ablations for other hyperparameters are in Appendix F.

Recent trends in ML have stressed the importance of understanding robustness under distribution shift, when training and testing distributions are not i.i.d. For time series data, for example, models may be trained on EEG data from one hospital, but deployed at another using instruments with different sampling rates ; or a time series may involve the same trajectory evolving at different speeds. Following Kidger et al. , we consider the Character Trajectories dataset , where the goal is to classify a character from a sequence of pen stroke measurements, collected from one user at a fixed sampling rate. To emulate timescale shift (e.g. testing on another user with slower handwriting), we consider two standard time series generation processes: (1) In the setting of sampling an underlying sequence at a fixed rate, we change the test sampling rate; crucially, the sequences are variable length so the models are unable to detect the sampling rate of the data. (2) In the setting of irregular-sampled (or missing) data with timestamps, we scale the test timestamps.

Recall that the HiPPO framework models the underlying data as a continuous function and interacts with discrete input only through the discretization. Thus, it seamlessly handles missing or irregularly-sampled data by simply evolving according to the given discretization step sizes (details in Section B.3). Combined with LegS timescale invariance (Prop. 3), we expect HiPPO-LegS to work automatically in all these settings. We note that the setting of missing data is a topic of independent interest and we compare against SOTA methods, including the GRU-D which learns a decay between observations, and neural ODE methods which models segments between observations with an ODE.

Table 2 validates that standard models can go catastrophically wrong when tested on sequences at different timescales than expected. Though all methods achieve near-perfect accuracy ($\geq$ 95%) without distribution shift, aside from HiPPO-LegS, no method is able to generalize to unseen timescales.

3 Theoretical Validation and Scalability

We empirically show that HiPPO-LegS can scale to capture dependencies across millions of time steps, and its memory updates are computationally efficient (processing up to 470,000 time steps/s).

We test the ability of different memory mechanisms in approximating an input function, as described in the problem setup in Section 2.1. The model only consists of the memory update (Section 3) and not the additional RNN architecture. We choose random samples from a continuous-time band-limited white noise process, with length $10^{6}$. The model is to traverse the input sequence, and then asked to reconstruct the input, while maintaining no more than 256 units in memory (Fig. 3). This is a difficult task; the LSTM fails with even sequences of length 1000 (MSE $\approx$ 0.25). As shown in Table 3, both the LMU and HiPPO-LegS are able to accurately reconstruct the input function, validating that HiPPO can solve the function approximation problem even for very long sequences. Fig. 3 illustrates the function and its approximations, with HiPPO-LegS almost matching the input function while LSTM unable to do so.

HiPPO-LegS operator is computationally efficient both in theory (Section 3) and in practice. We implement the fast update in C++ with Pytorch binding and show in Table 3 that it can perform 470,000 time step updates per second on a single CPU core, 10x faster than the LSTM and LMU.The LMU is only known to be fast with the simple forward Euler discretization , but not with more sophisticated methods such as bilinear and ZOH that are required to reduce numerical errors for this task.

4 Additional Experiments

We validate that the HiPPO memory updates also perform well on more generic sequence prediction tasks not exclusively focused on memory. Full results and details for these tasks are in Appendix F.

Our RNNs with HiPPO memory updates perform on par with the LSTM, while other long-range memory approaches such as expRNN perform poorly on this more generic task (Section F.6).

This physical simulation task tests the ability to model chaotic dynamical systems. HiPPO-LegS outperforms the LSTM, LMU, and the best hybrid LSTM+LMU model from , reducing normalized MSE by $30\%$ (Section F.7).

Conclusion

We address the fundamental problem of memory in sequential data by proposing a framework (HiPPO) that poses the abstraction of optimal function approximation with respect to time-varying measures. In addition to unifying and explaining existing memory approaches, HiPPO unlocks a new method (HiPPO-LegS) that takes a first step toward timescale robustness and can efficiently handle dependencies across millions of time steps. We anticipate that the study of this core problem will be useful in improving a variety of sequence models, and are excited about future work on integrating our memory mechanisms with other models in addition to RNNs. We hope to realize the benefits of long-range memory on large-scale tasks such as speech recognition, video processing, and reinforcement learning.

We thank Avner May, Mayee Chen, Dan Fu, Aditya Grover, and Daniel Lévy for their helpful feedback. We gratefully acknowledge the support of DARPA under Nos. FA87501720095 (D3M), FA86501827865 (SDH), and FA86501827882 (ASED); NIH under No. U54EB020405 (Mobilize), NSF under Nos. CCF1763315 (Beyond Sparsity), CCF1563078 (Volume to Velocity), and 1937301 (RTML); ONR under No. N000141712266 (Unifying Weak Supervision); the Moore Foundation, NXP, Xilinx, LETI-CEA, Intel, IBM, Microsoft, NEC, Toshiba, TSMC, ARM, Hitachi, BASF, Accenture, Ericsson, Qualcomm, Analog Devices, the Okawa Foundation, American Family Insurance, Google Cloud, Stanford HAI AWS cloud credit, Swiss Re, and members of the Stanford DAWN project: Teradata, Facebook, Google, Ant Financial, NEC, VMWare, and Infosys. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views, policies, or endorsements, either expressed or implied, of DARPA, NIH, ONR, or the U.S. Government. Atri Rudra’s research is supported by NSF grant CCF-1763481.

References

Appendix A Related Work

Our work touches on a variety of topics and related work, which we explore in detail.

The technical contributions in this work build on a rich history of approximation theory in signal processing. Our main framework – orthogonalizing functions with respect to time-varying measures (Section 2) – are related to “online” versions of classical signal processing transforms. In short, these methods compute specific transforms on sliding windows of discrete sequences. Concretely, they calculate $c_{n,k}=\sum_{i=0}^{N-1}f_{k+i}\psi(i,n)$ given signal $(f_{k})$, where $\{\psi(i,n)\}$ is a discrete orthogonal transform. Our technical problem differs in several key aspects:

Examples of sliding transforms considered in the literature include the sliding DFT , sliding DCT , sliding discrete (Walsh-)Hadamard transform , Haar , sliding discrete Hartley transform , and sliding discrete Chebyshev moments . While each of these address a specific transform, we present a general approach (Section 2) that addresses several transforms at once. Furthermore, we are unaware of sliding transform algorithms for the OPs we consider here, in particular the Legendre and Laguerre polynomials. Our derivations in Appendix D cover Legendre, (generalized) Laguerre, Fourier, and Chebyshev continuous sliding transforms.

All mentioned works operate in the sliding window setting, where a fixed-size context window on the discrete signal is taken into account. Our measure-based abstraction for approximation allows considering a new type of scaled measure where the window size increases over time, leading to methods with qualitatively different theoretical (Section 3) and empirical properties (Section 4.2). We are not aware of any previous works addressing this scaled setting.

Even in the fixed-length sliding window case, our solutions to the “translated measure” problems (e.g., HiPPO-LegT Section D.1) solve a continuous-time sliding window problem on an underlying continuous signal, then discretize.

On the other hand, the sliding transform problems calculate transforms directly on a discrete stream. Discrete transforms are equivalent to calculating projection coefficients on a measure (equation (18)) by Gaussian quadrature, which assumes the discrete input is subsampled from a signal at the quadrature nodes . However, since these nodes are non-uniformly spaced in general, the sliding discrete transform is not consistent with a discretization of an underlying continuous signal.

Thus, our main abstraction (Definition 1) has a fundamentally different interpretation than standard transforms, and our approach of first calculating the dynamics of the underlying continuous-time problem (e.g. equation (20)) is correspondingly new.

We remark that our novel scaled measures are fundamentally difficult to address with a standard discrete-time based approach. These discrete sliding methods require a fixed-size context in order to have consistent transform sizes, while the scaled measure would require solving transforms with an increasing number of input points over time.

A.1.2 OPs in ML

More broadly, orthogonal polynomials and orthogonal polynomial transforms have recently found applications in various facets of machine learning. For example, Dao et al. leverage the connection between orthogonal polynomials and quadrature to derive rules for computing kernel features in machine learning. More directly, apply parametrized families of structured matrices directly inspired by orthogonal polynomial transforms () as layers in neural networks. Some particular families of orthogonal polynomials such as the Chebyshev polynomials have desirable approximation properties that find many well-known classical uses in numerical analysis and optimization. More recently, they have been applied to ML models such as graph convolutional neural networks, and generalizations such as Gegenbauer and Jacobi polynomials have been used to analyze optimization dynamics. Generalization of orthogonal polynomials and Fourier transform, expressed as products of butterfly matrices, have found applications in automatic algorithm design , model compression , and replacing hand-crafted preprocessing in speech recognition . Orthogonal polynomials are known to have various efficiency results , and we conjecture that Proposition 4 on the efficiency of HiPPO methods can be extended to arbitrary measures besides the ones considered in this work.

A.2 Memory in Machine Learning

Sequential or temporal data in areas such as language, reinforcement learning, and continual learning can involve increasingly long dependencies. However, direct parametric modeling cannot handle inputs of unknown and potentially unbounded lengths. Many modern solutions such as attention and dilated convolutions , are functions on finite windows, thus sidestepping the need for an explicit memory representation. While this suffices for certain tasks, these approaches can only process a finite context window instead of an entire sequence. Naively increasing the window length poses significant compute and memory challenges. This has spurred various approaches to extend this fixed context window subjected to compute and storage constraints .

We instead focus on the core problem of online processing and memorization of continuous and discrete signals, and anticipate that the study of this foundational problem will be useful in improving a variety of models.

Recurrent neural networks are a natural tool for modeling sequential data online, with the appealing property of having unbounded context; in other words they can summarize history indefinitely. However, due to difficulties in the optimization process (vanishing/exploding gradients ), particular care must be paid to endow them with longer memory. The ubiquitous LSTM and simplifications such as the GRU control the update with gates to smooth the optimization process. With more careful parametrization, the addition of gates alone make RNNs significantly more robust and able to address long-term dependencies . Tallec and Ollivier show that gates are in fact fundamental for recurrent dynamics by allowing time dilations. Many other approaches to endowing RNNs with better memory exist, such as noise injection or non-saturating gates , which can suffer from instability issues. A long line of work controls the spectrum of the recurrent updates with (nearly-) orthogonal matrices to control gradients , but have been found to be less robust across different tasks .

A.3 Directly related methods

The main result of the Legendre Memory Unit is a direct instantiation of our framework using the LegT measure (Section 2.3). The original LMU is motivated by neurobiological advances and approaches the problem from the opposite direction as us: it considers approximating spiking neurons in the frequency domain, while we directly solve an interpretable optimization problem in the time domain. More specifically, they consider time-lagged linear time invariant (LTI) dynamical systems and approximate the dynamics with Padé approximants; Voelker et al. observes that the result also has an interpretation in terms of Legendre polynomials, but not that it is the optimal solution to a natural projection problem. This approach involves heavier machinery, and we were not able to find a complete proof of the update mechanism .

In contrast, our approach directly poses the relevant online signal approximation problem, which ties to orthogonal polynomial families and leads to simple derivations of several related memory mechanisms (Appendix D). Our interpretation in time rather than frequency space, and associated derivation (Section D.1) for the LegT measure, reveals a different set of approximations stemming from the sliding window, which is confirmed empirically (Section F.8).

As the motivations of our work are substantially different from Voelker et al. , yet finds the same memory mechanism in a special case, we highlight the potential connection between these sequence models and biological nervous systems as an area of exploration for future work, such as alternative interpretations of our methods in the frequency domain.

We remark that the term LMU in fact refers to a specific recurrent neural network architecture, which interleaves the projection operator with other specific neural network components. By contrast, we use HiPPO to refer to the projection operator in isolation (Theorem 1), which is a function-to-function or sequence-to-sequence operator independent of model. HiPPO is integrated into an RNN architecture in Section 4, with slight improvements to the LMU architecture, as ablated in Sections F.2 and F.3. As a standalone module, HiPPO can be used as a layer in other types of models.

The Fourier Recurrent Unit (FRU) uses Fourier basis (cosine and sine) to express the input signal, motivated by the discrete Fourier transform. In particular, each recurrent unit computes the discrete Fourier transform of the input signal for a randomly chosen frequency. It is not clear how discrete transform with respect to other bases (e.g., Legendre, Laguerre, Chebyshev) can in turn yield similar memory mechanisms. We show that FRU is also an instantiation of the HiPPO framework (Section D.4), where the Fourier basis can be viewed as orthogonal polynomials $z^{n}$ on the unit circle $\{z\colon\left\lvert z\right\rvert=1\}$.

Zhang et al. prove that if a timescale hyperparameter is chosen appropriately, FRU has bounded gradients, thus avoiding vanishing and exploding gradients. This essentially follows from the fact that $(1-\Delta t)^{T}=\Theta(1)$ if the discretization step size $\Delta t=\Theta(\frac{1}{T})$ is chosen, if the time horizon $T$ is known (cf. Sections B.3 and E). It is easily shown that this property is not intrinsic to the FRU but to sliding window methods, and is shared by all of our translated measure HiPPO methods (all but HiPPO-LegS in Appendix D). We show the stronger property that HiPPO-LegS, which uses scaling rather than sliding windows, also enjoys bounded gradient guarantees, without needing a well-specified timescale hyperparameter (Proposition 5).

HiPPO produces linear ODEs that describe the dynamics of the coefficients. Recent work has also incorporated ODEs into machine learning models. Chen et al. introduce neural ODEs, employing general nonlinear ODEs parameterized by neural networks in the context of normalizing flows and time series modeling. Neural ODEs have shown promising results in modeling irregularly sampled time series , especially when combined with RNNs . Though neural ODEs are expressive , due to their complex parameterization, they often suffer from slow training because of their need for more complicated ODE solvers. On the other hand, HiPPO ODEs are linear and are fast to solve with classical discretization techniques in linear systems, such as Euler method, Bilinear method, and Zero-Order Hold (ZOH) .

Appendix B Technical Preliminaries

We collect here some technical background that will be used in presenting the general HiPPO framework and in deriving specific HiPPO update rules.

Classical OPs families comprise Jacobi (which include Legendre and Chebyshev polynomials as special cases), Laguerre, and Hermite polynomials. The Fourier basis can also be interpreted as OPs on the unit circle in the complex plane.

We will also consider scaling the Legendre polynomials to be orthogonal on the interval $[0,t]$. A change of variables on (5) yields

Therefore, with respect to the measure $\omega_{t}=\mathbf{1}_{[0,t]}/t$ (which is a probability measure for all $t$), the normalized orthogonal polynomials are

In general, the orthonormal basis for any uniform measure consists of $(2n+1)^{\frac{1}{2}}$ times the corresponding linearly shifted version of $P_{n}$.

We note the following recurrence relations on Legendre polynomials ([2, Chapter 12]):

where the sum stops at $P_{0}$ or $P_{1}$.

These will be used in the derivations of the HiPPO-LegT and HiPPO-LegS updates, respectively.

B.1.2 Properties of Laguerre Polynomials

The standard Laguerre polynomials $L_{n}(x)$ are defined to be orthogonal with respect to the weight function $e^{-x}$ supported on $[0,\infty)$, while the generalized Laguerre polynomials (also called associated Laguerre polynomials) $L_{n}^{(\alpha)}$ are defined to be orthogonal with respect to the weight function $x^{\alpha}e^{-x}$ also supported on $[0,\infty)$:

The standard Laguerre polynomials correspond to the case of $\alpha=0$ of generalized Laguerre polynomials.

We note the following recurrence relations on generalized Laguerre polynomials ([2, Chapter 13.2]):

B.1.3 Properties of Chebyshev polynomials

Let $T_{n}$ be the classical Chebyshev polynomials (of the first kind), defined to be orthogonal with respect to the weight function $(1-x^{2})^{1/2}$ supported on $$, and let$p_{n}$be the normalized version of$T_{n}$ (i.e, with norm 1):

We will also consider shifting and scaling the Chebyshev polynomials to be orthogonal on the interval $[t-\theta,t]$ for fixed length $\theta$.

In terms of the original Chebyshev polynomials, these are

B.2 Leibniz Integral Rule

Differentiating through such integrals can be formalized by the Leibniz integral rule, the basic version of which states that

We elide over the formalisms in our derivations (Appendix D) and instead use the following trick. We replace integrand limits with an indicator function; and using the Dirac delta function $\delta$ when differentiating (i.e., using the formalism of distributional derivatives). For example, the above formula can be derived succinctly with this trick:

B.3 ODE Discretization

In our framework, time series inputs will be modeled with a continuous function and then discretized. Here we provide some background on ODE discretization methods, including a new discretization that applies to a specific type of ODE that our new method encounters.

The general formulation of an ODE is $\frac{d}{dt}c(t)=f(t,c(t))$. We will also focus on the linear time-invariant ODE of the form $\frac{d}{dt}c(t)=Ac(t)+Bf(t)$ for some input function $f(t)$, as a special case. The general methodology for discretizing the ODE, for step size $\Delta t$, is to rewrite the ODE as

Many ODE discretization methods corresponds to different ways to approximate the RHS integral:

To approximate the RHS of equation (12), keep the left endpoint $\Delta tf(t,c(t))$. For the linear ODE, we get:

To approximate the RHS of equation (12), keep the right endpoint $\Delta tf(t+\Delta t,c(t+\Delta t))$. For the linear ODE, we get the linear equation and the update:

To approximate the RHS of equation (12), average the endpoints $\Delta t\frac{f(t,c(t))+f(t+\Delta t,c(t+\Delta t))}{2}$. For the linear ODE, again we get a linear equation and the update:

This method approximates the RHS of equation (12) by taking a weighted average of the endpoints $\Delta t[(1-\alpha)f(t,c(t))+\alpha f(t+\Delta t,c(t+\Delta t))]$, for some parameter $\alpha\in$. For the linear ODE, again we get a linear equation and the update:

GBT generalizes the three methods mentioned above: forward Euler corresponds to $\alpha=0$, backward Euler to $\alpha=1$, and bilinear to $\alpha=1/2$.

In the case of HiPPO-LegS, we have a linear ODE of the form $\frac{d}{dt}c(t)=\frac{1}{t}Ac(t)+\frac{1}{t}Bf(t)$. Adapting the GBT discretization (which generalizes forward/backward Euler and bilinear) to this linear ODE, we obtain:

We highlight that this system is invariant to the discretization step size $\Delta t$. Indeed, if $c^{(k)}\coloneqq c(k\Delta t)$ and $f_{k}\coloneqq f(k\Delta t)$ then we have the recurrence

To understand the impact of approximation error in discretization, in Fig. 4, we show the absolute error for the HiPPO-LegS updates in function approximation (Section F.8) for different discretization methods: forward Euler, backward Euler, and bilinear. The bilinear method generally provide sufficiently accurate approximation. We will use bilinear as the discretization method for the LegS updates for the experiments.

Appendix C General HiPPO Framework

We present the general HiPPO framework, as described in Section 2, in more details. We also generalize it to include bases other than polynomials.

The first step refers to the $\operatorname{proj}_{t}$ operator and the second the $\operatorname{coef}_{t}$ operator in Definition 1.

We first describe the parameters of the $\operatorname{hippo}$ operator (a measure and basis) in more detail in Section C.1. We define the projection $\operatorname{proj}_{t}$ and coefficient $\operatorname{coef}_{t}$ operators in Section C.2. Then we give a general strategy to calculate these coefficients $c(t)$, by deriving a differential equation that governs the coefficient dynamics (Section C.3). Finally we discuss how to turn the continuous $\operatorname{hippo}$ operator into a discrete one that can be applied to sequence data (Section C.4).

We describe and motivate the ingredients of HiPPO in more detail here. Recall that the high level goal is online function approximation; this requires both a set of valid approximations and a notion of approximation quality.

We also assume for simplicity that the measures $\mu^{(t)}$ are normalized to be probability measures; arbitrary scaling does not affect the optimal projection.

Note that the $P_{n}^{(t)}$ are not required to be normalized, while the $p_{n}^{(t)}$ are.

Tilted measure and basis Our goal is simply to store a compressed representation of functions, which can use any basis, not necessarily OPs. For any scaling function

the functions $p_{n}(x)\chi(x)$ are orthogonal with respect to the density $\omega/\chi^{2}$ at every time $t$. Thus, we can choose this alternative basis and measure to perform the projections.

To formalize this tilting with $\chi$, define $\nu^{(t)}$ to be the normalized measure with density proportional to $\omega^{(t)}/(\chi^{(t)})^{2}$.

We will calculate the normalized measure and the orthonormal basis for it. Let

be the normalization constant, so that $\nu^{(t)}$ has density $\frac{\omega^{(t)}}{\zeta(t)(\chi^{(t)})^{2}}$. If $\chi(t,x)=1$ (no tilting), this constant is $\zeta(t)=1$. In general, we assume that $\zeta$ is constant for all $t$; if not, it can be folded into $\chi$ directly.

Next, note that (dropping the dependence on $x$ inside the integral for shorthand)

Thus we define the orthogonal basis for $\nu^{(t)}$

We let each element of the basis be scaled by a $\lambda_{n}$ scalar, for reasons discussed soon, since arbitrary scaling does not change orthogonality:

Note that when $\lambda_{n}=\pm 1$, the basis $\{g_{n}^{(t)}\}$ is an orthonormal basis with respect to the measure $\nu^{(t)}$, at every time $t$. Notationally, let $g_{n}(t,x):=g_{n}^{(t)}(x)$ as usual.

We will only use this tilting in the case of Laguerre (Section D.2 and Chebyshev (Section D.5).

Note that in the case $\chi=1$ (i.e., no tilting), we also have $\zeta=1$ and $g_{n}=\lambda_{n}p_{n}$ (for all $t,x$).

C.2 The Projection and Coefficients

Given a choice of measures and basis functions, we next see how the coefficients $c(t)$ can be computed.

The function $f$ can be approximated by storing its coefficients with respect to the basis $\{g_{n}\}_{n<N}$. For example, in the case of no tilting $\chi=1$, this encodes the optimal polynomial approximation of $f$ of degree less than $N$. In particular, at time $t$ we wish to represent $f_{\leq t}$ as a linear combination of polynomials $g_{n}^{(t)}$. Since the $g_{n}^{(t)}$ are orthogonal with respect to the Hilbert space defined by $\langle\cdot,\cdot\rangle_{\nu^{(t)}}$, it suffices to calculate coefficients

At any time $t$, $f_{\leq t}$ can be explicitly reconstructed as

Equation (19) is the $\operatorname{proj}_{t}$ operator; given the measure and basis parameters, it defines the optimal approximation of $f_{\leq t}$.

The $\operatorname{coef}_{t}$ operator simply extracts the vector of coefficients $c(t)=(c_{n}(t))_{n\in[N]}$.

C.3 Coefficient Dynamics: the hippohippo\operatorname{hippo} Operator

In our framework, we will compute these coefficients over time by viewing them as a dynamical system. Differentiating (18),

Here we have made use of the assumption that $\zeta$ is constant for all $t$.

The key idea is that if $\frac{\partial}{\partial t}P_{n}$ and $\frac{\partial}{\partial t}\frac{\omega}{\chi}$ have closed forms that can be related back to the polynomials $P_{k}$, then an ordinary differential equation can be written for $c(t)$. This allows these coefficients $c(t)$ and hence the optimal polynomial approximation to be computed online. Since $\frac{d}{dt}P_{n}^{(t)}$ is a polynomial (in $x$) of degree $n-1$, it can be written as linear combinations of $P_{0},\dots,P_{n-1}$, so the first term in Eq. 20 is a linear combination of $c_{0},\dots,c_{n-1}$. For many weight functions $\omega$, we can find scaling function $\chi$ such that $\frac{\partial}{\partial t}\frac{\omega}{\chi}$ can also be written in terms of $\frac{\omega}{\chi}$ itself, and thus in those cases the second term of Eq. 20 is also a linear combination of $c_{0},\dots,c_{N-1}$ and the input $f$. Thus this often yields a closed-form linear ODE for $c(t)$.

Our purpose of defining the free parameters $\lambda_{n}$ was threefold.

First, note that the orthonormal basis is not unique, up to a $\pm 1$ factor per element.

Second, choosing $\lambda_{n}$ can help simplify the derivations.

Third, although choosing $\lambda_{n}=\pm 1$ will be our default, since projecting onto an orthonormal basis is most sensible, the LMU used a different scaling. Section D.1 will recover the LMU by choosing different $\lambda_{n}$ for the LegT measure.

Suppose that equation (20) reduced to dynamics of the form

Then, letting $\Lambda=\operatorname*{diag}_{n\in[N]}\{\lambda_{n}\}$,

Therefore, if we reparameterize the coefficients ($\Lambda^{-1}c(t)\to c(t)$) then the normalized coefficients projected onto the orthonormal basis satisfy dynamics and associated reconstruction

These are the $\operatorname{hippo}$ and $\operatorname{proj}_{t}$ operators.

C.4 Discretization

As defined here, $\operatorname{hippo}$ is a map on continuous functions. However, as $\operatorname{hippo}$ defines a closed-form ODE of the coefficient dynamics, standard ODE discretization methods (Section B.3) can be applied to turn this into discrete memory updates. Thus we overload these operators, i.e. $\operatorname{hippo}$ either defines an ODE of the form

Section F.5 validates the framework by applying (20) and (19) to approximate a synthetic function.

Appendix D Derivations of HiPPO Projection Operators

We derive the memory updates associated with the translated Legendre (LegT) and translated Laguerre (LagT) measures as presented in Section 2.3, along with the scaling Legendre (LegS) measure (Section 3). To show the generality of the framework, we also derive memory updates with Fourier basis (recovering the Fourier Recurrent Unit ) and with Chebyshev basis.

The majority of the work has already been accomplished by setting up the projection framework, and the proof simply requires following the technical outline laid out in Appendix C. In particular, the definition of the coefficients (18) and reconstruction (19) does not change, and we only consider how to calculate the coefficients dynamics (20).

For each case, we follow the general steps:

define the measure $\mu^{(t)}$ or weight $\omega(t,x)$ and basis functions $p_{n}(t,x)$,

compute the derivatives of the measure and basis functions,

plug them into the coefficient dynamics (equation (20)) to derive the ODE that describes how to compute the coefficients $c(t)$,

provide the complete formula to reconstruct an approximation to the function $f_{\leq t}$, which is the optimal projection under this measure and basis.

The derivations in Sections D.1 and D.2 prove Theorem 1, and the derivations in Section D.3 prove Theorem 2. Sections D.4 and D.5 show additional results for Fourier-based bases.

Figure 5 illustrates the overall framework when we use Legendre and Laguerre polynomials as the basis, contrasting our main families of time-varying measures $\mu^{(t)}$.

This measure fixes a window length $\theta$ and slides it across time.

We use a uniform weight function supported on the interval $[t-\theta,t]$ and pick Legendre polynomials $P_{n}(x)$, translated from $$to$[t-\theta,t]$, as basis functions:

Here, we have used no tilting so $\chi=1$ and $\zeta=1$ (equations (15) and (16)). We leave $\lambda_{n}$ unspecified for now.

At the endpoints, these basis functions satisfy

The derivative of Legendre polynomials can be expressed as linear combinations of other Legendre polynomials (cf. Section B.1.1).

As a special case for the LegT measure, we need to consider an approximation due to the nature of the sliding window measure.

When analyzing $\frac{d}{dt}c(t)$ in the next section, we will need to use the value $f(t-\theta)$. However, at time $t$ this input is no longer available. Instead, we need to rely on our compressed representation of the function: by the reconstruction equation (19), if the approximation is succeeding so far, we should have

We are ready to derive the coefficient dynamics.

Plugging the derivatives of this measure and basis into equation (20) gives

Now we consider two instantiations for $\lambda_{n}$. The first one is the more natural $\lambda_{n}=1$, which turns $g_{n}$ into an orthonormal basis. We then get

The second case takes $\lambda_{n}=(2n+1)^{\frac{1}{2}}(-1)^{n}$. This yields

By equation (19), at every time $t$ we have

D.2 Derivation for Translated Laguerre (HiPPO-LagT)

The result in Theorem 1 for HiPPO-LagT is for the case $\alpha=0,\beta=1$, corresponding to the basic Laguerre polynomials and no tilting.

We flip and translate the generalized Laguerre weight function and polynomials from $[0,\infty)$ to $(-\infty,t]$. The normalization is found using equation (9).

to be the norm of the generalized Laguerre polynomial $L_{n}^{(\alpha)}$, so that $\lambda_{n}p_{n}^{(t)}=L_{n}^{(\alpha)}(t-x)$, and (following equation (17)) the basis for $\nu^{(t)}$ is

The derivative of Laguerre polynomials can be expressed as linear combinations of other Laguerre polynomials (cf. Section B.1.2).

Plugging these derivatives into equation (20) (obtained from differentiating the coefficient equation (18)), where we suppress the dependence on $x$ for convenience:

By equation (19), at every time $t$, for $x\leq t$,

Finally, following equations (21) and (22) to convert these to dynamics on the orthonormal basis of the normalized (probability) measure $\nu^{(t)}$ leads to the following $\operatorname{hippo}$ operator

and correspondingly a $\operatorname{proj}_{t}$ operator:

D.3 Derivation for Scaled Legendre (HiPPO-LegS)

As discussed in Section 3, the scaled Legendre is our only method that uses a measure with varying width.

Here, $P_{n}$ are the basic Legendre polynomials (Section B.1.1). We use no tilting, i.e. $\chi(t,x)=1$, $\zeta(t)=1$, and $\lambda_{n}=1$ so that the functions $g_{n}(t,x)$ are an orthonormal basis.

We first differentiate the measure and basis:

Now define $z=\frac{2x}{t}-1$ for shorthand and apply the properties of derivatives of Legendre polynomials (equation (8)).

where we have used $g_{n}(t,t)=(2n+1)^{\frac{1}{2}}P_{n}(1)=(2n+1)^{\frac{1}{2}}$. Vectorizing this yields equation (3):

where $D:=\operatorname*{diag}\left[(2n+1)^{\frac{1}{2}}\right]_{n=0}^{N-1}$, $\mathbf{1}$ is the all ones vector, and the state matrix $M$ is

Equation (29) is a linear dynamical system, except dilated by a time-varying factor $t^{-1}$, which arises from the scaled measure.

By equation (19), at every time $t$ we have

D.4 Derivation for Fourier Bases

In the remainder of Appendix D, we consider some additional bases which are analyzable under the HiPPO framework. These use measures and bases related to various forms of the Fourier transform.

Similar to the LMU, the sliding Fourier measure also has a fixed window length $\theta$ parameter and slides it across time.

For each $t$, we will use a sliding measure uniform on $[t-\theta,t]$ and rescale the basis as $e^{2\pi in\frac{t-x}{\theta}}$ (so they are still orthonormal, i.e., have norm 1):

Note that $p_{n}(t,t)=p_{n}(t,t-\theta)=1$. Additionally, we no longer have access to $f(t-\theta)$ at time $t$, but this is implicitly represented in our compressed representation of the function: $f=\sum_{k=0}^{N-1}c_{k}(t)p_{k}(t)$. Thus we approximate $f(t-\theta)$ by $\sum_{k=0}^{N-1}c_{k}(t)p_{k}(t,t-\theta)=\sum_{k=0}^{N-1}c_{k}(t)$. Finally, this yields

Hence $\frac{d}{dt}c(t)=Ac(t)+Bf(t)$ where

D.4.2 Fourier Recurrent Unit

Using the HiPPO framework, we can also derive the Fourier Recurrent Unit (FRU) .

For each $t$, we will use a sliding measure uniform on $[t-\theta,t]$ and the basis $e^{2\pi in\frac{x}{\theta}}$:

In general the basis is not orthogonal with respect to the measure $\omega(t,x)$, but orthogonality holds at the end where $t=\theta$.

We no longer have access to $f(t-\theta)$ at time $t$, but we can approximate by ignoring this term (which can be justified by assuming that the function $f$ is only defined on $[0,\theta]$ and thus $f(x)$ can be set to zero for $x<0$). Finally, this yields

Applying forward Euler discretization (with step size = 1), we obtain:

Taking the real parts yields the Fourier Recurrent Unit updates .

Note that the recurrence is independent in each $n$, so we don’t have the pick $n=0,1,\dots,N-1$. We can thus pick random frequencies $n$ as done in Zhang et al. .

D.5 Derivation for Translated Chebyshev

The final family of orthogonal polynomials we analyze under the HiPPO framework are the Chebyshev polynomials. The Chebyshev polynomials can be seen as the purely real analog of the Fourier basis; for example, a Chebyshev series is related to a Fourier cosine series through a change of basis .

Note that at the endpoints, these evaluate to

We also choose $\lambda_{n}=1$ for the canonical orthonormal basis, so

We consider differentiating the polynomials separately for $n=0$, $n$ even, and $n$ odd, using equation (11). Defined $z=\frac{2(x-t)}{\theta}+1$ for convenience. First, for $n$ even,

where we take $f(t-\theta)=0$ as we no longer have access to it (this holds when $t<\theta$ as well).

In the usual way, we can write this as linear dynamics

Appendix E HiPPO-LegS Theoretical Properties

The second-to-last equality uses the change of variables $x\mapsto\frac{x}{\alpha}$. ∎

E.2 Speed

In this section we work out the fast update rules according to the forward Euler, backward Euler, bilinear, or generalized bilinear transform discretizations (cf. Section B.3). Recall that we must be able to perform matrix-vector multiplication by $I+\delta A$ and $(I-\delta A)^{-1}$ where $\delta$ is some multiple of the step size $\Delta t$ (equation (13)).

It is easily seen that the LegS update rule involves a matrix $A$ of the following form (Theorem 2): $A=D_{1}(L+D_{0})D_{2}$, where $L$ is the all $1$ lower triangular matrix and $D_{0},D_{1},D_{2}$ are diagonal. Clearly, $I+\delta A$ is efficient (only requiring $O(N)$ operations), as it only involves matrix-vector multiplication by diagonals $D_{0},D_{1},D_{2}$, or multiplication by $L$ which is the $\mathsf{cumsum}$ operation.

Now we consider multiplication by the inverse $(I+\delta A)^{-1}$ (the minus sign can be absorbed into $\delta$). Write

Since diagonal multiplication is efficient, the crucial operation is inversion multiplication by a matrix of the form $L+D$.

Consider solving the equation $(L+D)x=y$. This implies $x_{0}+\dots+x_{k-1}=y_{k}-(1+d_{k})x_{k}$. The solution is

Finally, consider how to calculate a recurrence of the following form efficiently.

Evidently $x$ can be computed in a vectorized way as

E.3 Gradient Norms

We analyze the discrete time case under the Euler discretization (Section B.3), where the HiPPO-LegS recurrent update is equation (4), restated here for convenience:

These gradient asymptotics hold under other discretizations.

The problem reduces to showing that $\rho=\Theta(1/l)$.

We will use the following facts about the function $\log\left(1-\frac{1}{x}\right)$. First, it is an increasing function, so

Finally, note that $x\log\left(1-\frac{1}{x}\right)$ is an increasing function, and bounded from above since it is negative, so it is $\Theta(1)$ (this can also be seen from its Taylor expansion). Thus we have

E.4 Function Approximation Error

Since $p_{n}^{(t)}$ forms an orthonormal basis of the Hilbert space defined by the inner product $\langle\cdot,\cdot\rangle_{\mu^{(t)}}$ , by Parseval’s identity,

To bound the error $\left\|{f_{\leq t}-g^{(t)}}\right\|_{\mu^{(t)}}$, it suffices to bound the sum of the squares of the high-order coefficients $c_{n}(t)$ for $n=N,N+1,\dots$. We will bound each coefficient by integration by parts.

We first simplify the expression for $c_{n}(t)$. For any $n\geq 1$, we have

As $P_{n}(x)=\frac{1}{2n+1}\frac{d}{dx}(P_{n+1}(x)-P_{n-1}(x))$ (cf. Section B.1.1), integration by parts yields:

Notice that the boundary term is zero, since $P_{n+1}(1)=P_{n-1}(1)=1$ and $P_{n+1}(-1)=P_{n-1}(-1)=\pm 1$ (either both 1 or both $-1$ depending on whether $n$ is odd or even). Hence:

Now suppose that $f$ is $L$-Lipschitz, which implies that $\left\lvert f^{\prime}\right\rvert\leq L$. Then

We then obtain that $\left\|{f_{\leq t}-g^{(t)}}\right\|_{\mu^{(t)}}=O(tL/\sqrt{N})$ as claimed.

Now supposed that $f$ has $k$ derivatives and the $k$-th derivative is bounded. The argument is similar to the one above where we integrate by parts $k$ times. We sketch this argument here.

Take $k$ to be a constant, and let $n\geq k$. Applying integration by parts $k$ times, noting that all the boundary terms are zero, gives:

We then obtain that $\left\|{f_{\leq t}-g^{(t)}}\right\|_{\mu^{(t)}}=O(t^{k}N^{-k+1/2})$ as claimed. ∎

Remark. The approximation error of Legendre polynomials reduces to how fast the Legendre coefficients decay, subjected to the smoothness assumption of the input function. This result is analogous to the classical result in Fourier analysis, where the $n$-th Fourier coefficients decay as $O(n^{-k})$ if the input function has order-$k$ bounded derivatives . That result is also proved by integration by parts.

Appendix F Experiment Details and Additional Results

Given inputs $x_{t}$ or features thereof $f(x_{t})$ in any model, the HiPPO framework can be used to memorize the history of features $f_{t}$ through time. As the discretized HiPPO dynamics form a linear recurrent update similar in style to RNNs (e.g., Theorem 2), we focus on these models in our experiments.

Thus, given any RNN update function $h_{t}=\tau(h_{t-1},x_{t})$, we simply replace the previous hidden state with a projected version of its entire history.

Equations (31) lists the explicit update equations and Figure 6 illustrates the model. In our experiments, we choose a basic gated RNN update

We consider the following instantiations of our framework HiPPO.

HiPPO-LegT, LagT, and LegS, use the translated Legendre, and tilted Laguerre, and scaled Legendre measure families with update dynamics (1), (2), and (3). As mentioned, LegT has an additional hyperparameter $\theta$, which should be set to the timescale of the data if known a priori. We attempt to set it equal to its ideal value (the length of the sequences) in every task, and also consider $\theta$ values that are too large and small to illustrate the effect of this hyperparameter.

Our derivations in Sections D.1, D.2, D.3, D.4 and D.5 show that there is a large variety of update equations that can arise from the HiPPO framework—for example, the tilted generalized Laguerre polynomials lead to an entire family governed by two free parameters (Section D.2)—many of which lead to linear dynamics of the form $\frac{d}{dt}c(t)=-Ac(t)+Bf(t)$ for various $A,B$. Given that many different update dynamics lead to such dynamical systems that give sensible results, we additionally consider the HiPPO-Rand baseline that uses random $A$ and $B$ matrices (normalized appropriately) in its dynamics.

We additionally compare against the following standard RNN baselines. The RNN is a vanilla RNN. The MGU is a minimal gated architecture, equivalent to a GRU without the reset gate. The HiPPO architecture we use is simply the MGU with an additional $\operatorname{hippo}$ intermediate layer. The LSTM is the most well-known and popular RNN architecture, which is a more sophisticated gated RNN. The expRNN is the state-of-the-art representative of the orthogonal RNN family of models designed for long-term dependencies . The LMU is the exact same model as in Voelker et al. ; it is equivalent to HiPPO-LegT with a different RNN architecture.

All methods have the same hidden size in our experiments. In particular, for simplicity and to reduce hyperparameters, HiPPO variants tie the memory size $N$ to the hidden state dimension $d$. The hyperparameter $N$ and $d$ is also referred to as the number of hidden units.

The model (31) we use is a simple RNN that bears similarity to the classical LSTM and the original LMU cell. In comparison to the LSTM, HiPPO can be seen as a variant where the memory $m_{t}$ plays the role of the LSTM’s hidden state and $h_{t}$ plays the role of the LSTM’s gated cell state, with equal dimensionalities. HiPPO updates $m_{t}$ using the fixed $A$ transition matrix instead of a learned matrix, and also lacks “input” and “output” gates, so for a given hidden size, it requires about half the parameters.

The LMU is a version of the HiPPO-LegT cell with an additional hidden-to-hidden transition matrix and memory-to-memory transition vector instead of the gate $g$, leaving it with approximately the same number of trainable parameters.

Unless stated otherwise, all methods use the Adam optimizer with learning rate frozen to $0.001$, which has been a robust default for RNN based models .

All experiments use PyTorch 1.5 and are run on a Nvidia P100 GPU.

F.2 Permuted MNIST

The input to the sequential MNIST (sMNIST) task is an MNIST source image, flattened in row-major order into a single sequence of length 784. The goal of the model is to process the entire image sequentially before outputting a classification label, requiring learning long-term dependencies. A variant of this, the permuted MNIST (pMNIST) task, applies a fixed permutation to every image, breaking locality and further straining a model’s capacity for long-term dependencies.

Models are trained using the cross-entropy loss. We use the standard train-test split (60,000 examples for training and 10,000 for testing), and further split the training set with 10% to be used as validation set.

Table 1 is duplicated here in Tables 4 and 5, with more complete baselines and hyperparameter ablations.

Table 4 consists of our implementations of various baselines related to our method, described in Section F.1. Each method was ran for 3 seeds, and the maximum average validation accuracy is reported.

All methods used the same hidden size of $512$; we found that this gave better performance than $256$, and further increasing it did not improve more. All methods were trained for 50 epochs with a batch size of 100.

Table 5 directly shows the reported test accuracy of various methods on this data (Middle and Bottom). Table 5 (Top) reports the test accuracy of various instantations of our methods. We additionally include our reproduction of the LMU, which achieved better results than reported in Voelker et al. (possibly due to a larger hidden size). We note that all of our HiPPO methods are competitive; each of them (HiPPO-LegT, HiPPO-LagT, HiPPO-LegS) achieves state-of-the-art among previous recurrent sequence models. Note that differences between our HiPPO-LegT and LMU numbers in Table 5 (Top) stem primarily from the architecture difference (Section F.1).

Table 4 also shows ablations for the HiPPO-LegT and HiPPO-LagT timescale hyperparameters. HiPPO-LagT sweeps the discretization step size $\Delta t$ (Sections 2.4 and B.3). For LegT, we set $\Delta t=1.0$ without loss of generality, as only the ratio of $\theta$ to $\Delta t$ matters. These timescale hyperparameters are important for these methods. Previous works have shown that the equivalent of $\Delta t$ in standard RNNs, i.e. the gates of LSTMs and GRUs (Section 2.4), can also drastically affect their performance . For example, the only difference between the URLSTM and LSTM in Table 5 is a reparametrization of the gates.

F.3 Copying

In the Copying task , the input is a sequence of $L+20$ digits where the first 10 tokens $\left(a_{0},a_{1},\dots,a_{9}\right)$ are randomly chosen from $\left\{1,\dots,8\right\}$, the middle N tokens are set to , and the last ten tokens are $9$. The goal of the recurrent model is to output $\left(a_{0},\dots,a_{9}\right)$ in order on the last 10 time steps, whenever the cue token $9$ is presented. Models are trained using the cross-entropy loss; the random guessing baseline has loss $\log(8)\approx 2.08$. We use length $L=200$. The training and testing examples are generated in the same way.

Our motivation of studying the Copying task is that standard models such as the LSTM struggle to solve it. We note that the Copying task is much harder than other memory benchmarks such as the Adding task , and we do not consider those.

The HiPPO-LegS method solves this task the fastest. The LegT method also solves this task quickly, only if the parameter $\theta$ is initialized to the correct value of $200$. Mis-specifying this timescale hyperparameter to $\theta=20$ or $\theta=2000$ drastically slows down the convergence of HiPPO-LegT. The LMU (at optimal parameter $\theta=200$) solves this task at comparable speed; like in Section F.2, differences between HiPPO-LegT $(\theta=200)$ and LMU here arise from the minor architecture difference in Section F.1.

The HiPPO-Rand baseline (denoted “random LTI” system here) does much worse than the updates with the dynamics derived from our framework, highlighting the importance of the precise dynamics (in contrast to just the architecture).

Standard methods such as the RNN and LSTM are also nearly stuck at baseline.

F.4 Trajectory Classification

The Character Trajectories dataset from the UCI machine learning repository consists of pen tip trajectories recorded from writing individual characters. The trajectories were captured at 200Hz and data was normalized and smoothed. Input is 3-dimensional ($x$ and $y$ positions, and pen tip force), and there are 20 possible outputs (number of classes). Models are trained using the cross-entropy loss. The dataset contains 2858 time series. The length of the sequences is variable, ranging up to $182$. We use a train-val-test split of 70%-15%-15%.

RNN baselines include the LSTM , GRU , and LMU . Our implementations of these used 256 hidden units each.

The GRU-D is a method for handling missing values in time series that computes a decay between observations. The ODE-RNN and Neural CDE (NCDE) baselines are state-of-the-art neural ODE methods, also designed to handle irregularly-sampled time series. Our GRU-D, ODE-RNN, and Neural CDE baselines used code from Kidger et al. , inheriting the hyperparameters for those methods.

The goal of this experiment is to investigate the performance of models when the timescale is mis-specified between train and evaluation time, leading to distribution shift. We considered the following two standard types of time series:

Irregularly-sampled time series (i.e., missing values) with timestamps

Timescale shift is emulated in the corresponding ways, which can be interpreted as different sampling rates or trajectory speeds.

Either the train or evaluation sequences are downsampled by a factor of 2

The train or evaluation timestamps are halved.Instead of the train timestamps being halved, equivalently the evaluation timestamps can be doubled.

The first scenario in each corresponds to the original sequence being sampled at 100Hz instead of 200Hz; alternatively, it is equivalent to the writer drawing twice as fast. Thus, these scenarios correspond to a train $\to$ evaluation timescale shift of 100Hz $\to$ 200Hz and 200Hz $\to$ 100Hz respectively.

Note that models are unable to obviously tell that there is timescale shift. For example, in the first scenario, shorter or longer sequences can be attributed to the variability of sequence lengths in the original dataset. In the second scenario, the timestamps have different distributions, but this can correspond to different rates of missing data, which the baselines for irregularly-sampled data are able to address.

F.5 Online Function Approximation and Speed Benchmark

The task is to reconstruct an input function (as a discrete sequence) based on some hidden state produced after the model has traversed the input function. This is the same problem setup as in Section 2.1; the online approximation and reconstruction details are in Appendix C. The input function is randomly sampled from a continuous-time band-limited white noise process, with length $10^{6}$. The sampling step size is $\Delta t=10^{-4}$, and the signal band limit is 1Hz.

We compare HiPPO-LegS, LMU, and LSTM. The HiPPO-LegS and LMU model only consists of the memory update and not the additional RNN architecture. The function is reconstructed from the coefficients using the formula in Appendix D, so no training is required. For LSTM, we use a linear decoder to reconstruct the function from the LSTM hidden states and cell states, trained on a collection of 100 sequences. All models use $N=256$ hidden units. The LSTM uses the $L2$ loss. The HiPPO methods including LMU follow the fixed dynamics of Theorem 1 and Theorem 2.

We measure the inference time of HiPPO-LegS, LMU, and LSTM, in single-threaded mode on a server Intel Xeon CPU E5-2690 v4 at 2.60GHz.

F.6 Sentiment Classification on the IMDB Movie Review Dataset

The IMDB movie review dataset is a standard binary sentiment classification task containing 25000 train and test sequences, with sequence lengths ranging from hundreds to thousands of steps. The task is to classify the sentiment of each movie review into either positive or negative. We use 10% of the standard training set as validation set.

RNN baselines include the LSTM , vanilla RNN, LMU , and expRNN . Our implementations of these used 256 hidden units each.

As shown in Table 6, our HiPPO-RNNs have similar and consistent performance, on par or better than LSTM. Other long-range memory RNN approaches that constrains the expressivity of the network (e.g. expRNN) performs worse on this more generic task.

F.7 Mackey Glass prediction

The Mackey-Glass data is a time series prediction task for modeling chaotic dynamical systems. We build on the implementation of Voelker et al. . The data is a sequence of one-dimensional observations, and models are tasked with predicting 15 time steps into the future. The models are 4-layer stacked recurrent neural networks, trained with the mean squared error (MSE) loss. Voelker et al. additionally consider a hybrid model with alternating LSTM and LMU layers, which improved on either by itself. We did not try this approach with our method HiPPO-LegS such as combining it with the LSTM or other HiPPO methods, but such ideas could further improve our performance. As a baseline method, the identity function does not simulate the dynamics, and simply guesses that the future time step is equal to the current input.

F.8 Additional Analysis and Ablations of HiPPO

To further analyze the tradeoffs of the memory updates derived from our framework, in Fig. 9 we plot a simple input function $f(x)=1/4\sin x+1/2\sin(x/3)+\sin(x/7)$ to be approximated. The function is subsampled on the range $x\in$, creating a sequence of length $1000$. This function is simpler than the functions sampled from white noise signals described in Section F.5. Given this function, we use the same methodology as in Section F.5 for processing the function online and then reconstructing it at the end.

In Figure 9(a, b), we plot the true function $f$, and its absolute approximation error based on LegT, LagT, and LegS. LegS has the lowest approximation error, while LegT and LagT are similar and slightly worse than LegS. Next, we analyze some qualitative behaviors.

In Figure 9(c), shows that the approximation error of LegT is sensitive to the hyperparameter $\theta$, the length of the window. Specifying $\theta$ to be even slightly too small (by $0.5\%$ relative to the total sequence length) causes huge errors in approximation. This is expected by the HiPPO framework, as the final measure $\mu^{(t)}$ is not supported everywhere, so the projection problem does not care that the reconstructed function is highly inaccurate near $x=0$.

Our LagT method actually comprises a family of related transforms, governed by two parameters $\alpha,\beta$ specifying the original measure and the tilting (Section D.2). Fig. 10 shows the error as these parameters change. Fig. 10(a) shows that small $\alpha$ generally performs better. Fig. 10(b, c) show that the reconstruction is unstable for larger $\beta$, but small values of $\beta$ work well. More detailed theoretical analysis explaining these tradeoffs would be an interesting question to analyze.