# LAUGHING HYENA DISTILLERY: Extracting Compact Recurrences From Convolutions Stefano Massaroli\*,1, Michael Poli\*,2, Daniel Y. Fu\*,2, Hermann Kumbong2, Rom N. Parnichkun3, Aman Timalsina4, David W. Romero5, Quinn McIntyre2, Beidi Chen6, Atri Rudra7, Ce Zhang8, Christopher Re2,†, Stefano Ermon2,†, Yoshua Bengio1,† NeurIPS 2023, Last Compiled: October 31, 2023. ## Abstract Recent advances in attention-free sequence models rely on convolutions as alternatives to the attention operator at the core of Transformers. In particular, *long* convolution sequence models have achieved state-of-the-art performance in many domains, but incur a significant cost during auto-regressive inference workloads – naively requiring a full pass (or *caching* of activations) over the input sequence for each generated token – similarly to attention-based models. In this paper, we seek to enable $\mathcal{O}(1)$ compute and memory cost per token in any pre-trained long convolution architecture to reduce memory footprint and increase throughput during generation. Concretely, our methods consist in extracting low-dimensional linear state-space models from each convolution layer, building upon rational interpolation and model-order reduction techniques. We further introduce architectural improvements to convolution-based layers such as **Hyena**: by weight-tying the filters across channels into *heads*, we achieve higher pre-training quality and reduce the number of filters to be distilled. The resulting model achieves $10\times$ higher throughput than Transformers and $1.5\times$ higher than **Hyena** at 1.3B parameters, without any loss in quality after distillation. ## 1 Introduction Attention-free approaches such as *long convolution sequence models* (LCSMs), e.g., H3 [1], Hyena [2], have shown promise in matching Transformer [3, 4] performance across a wide range of tasks, with sub-quadratic complexity with respect to sequence length. Despite the improved efficiency during training on long sequences, unless the convolution filters are either *short* or admit a *low*-dimensional state-state-space realization, LCSMs still need to process the entire growing sequence at every step of auto-regressive generation, similarly to Transformers. In this work, we seek to refine LCSMs in both **efficiency** and **quality**. First, we study the inference stage, and propose methods to enable a *recurrent* mode for auto-regressive generation. Recurrent modes prescribe the existence of a *state* encoding the past information of the process in a fixed-dimension memory, enabling **constant per-step time** and **constant-memory** in generation. Then, we draw upon an analysis of pre-trained models to develop architectural enhancements for the **Hyena** block, simultaneously improving model quality and efficiency of the distillation procedure. **Distilling fast recurrences** We introduce **LaughingHyena**, the first distillation approach for LCSMs that enables recurrent inference without impacting downstream quality. **LaughingHyena** seeks compact recurrences in the form of *state-space models* (SSMs) [5, 6] as the solution of a nonlinear interpolation problem involving the convolution filters of a pre-trained model. Since the total memory cost of SSMs grows linearly in the state dimension $d$ , our distillation procedure enables high throughput by enabling processing of large batches during generation. --- \*Equal contribution. †Equal senior authorship. 1Mila and Université de Montréal. 2Stanford University. 3The University of Tokyo. 4Purdue University. 5Vrije Universiteit Amsterdam. 6Carnegie Mellon University and Meta AI (FAIR). 7University of Buffalo, SUNY. 8University of Chicago and Together Computer.We identify and address three core challenges related to distillation, including the identification of: - • **Target state dimension:** we identify candidate state dimensions of our distilled SSMs by analyzing the spectrum of the Hankel operator associated with each convolution [7]. - • **Parametrization:** we address issues with naive parametrizations by introducing a factorized *modal* form, inspired by barycentric [8] and Prony-like [9] methods . - • **Approximation metric:** to ensure compatibility with any downstream task, we choose discrepancy metrics on the convolution filter, rather than model outputs. In auto-regressive workloads, **LaughingHyena**-distilled models with state dimension $d$ can generate $K$ tokens in $\mathcal{O}(dK)$ time and with constant $\mathcal{O}(d)$ memory – improving over the $\mathcal{O}(K^2)$ time and $\mathcal{O}(K)$ memory usage of *kv-cached* Transformers and naively executed long convolutions. At model sizes above one billion parameters, **LaughingHyena** achieves $10\times$ higher peak throughput over comparable Transformers (Figure 1.1), and can process larger batch sizes. Constant memory generation enables larger $K$ for a given a memory constraint e.g., generating 512 tokens with **LaughingHyena** requires $3\times$ less memory than with a Transformer. At smaller batch sizes, latency of **LaughingHyena** is also competitive with Transformers, reaching $\geq 2\times$ speedups at longer prompt lengths. **Improving pre-training quality** We leverage our analysis of the distillation process to open up new avenues of improvement for LCSM architectures. Indeed, the high compression rates achievable through **LaughingHyena** hint at sub-utilization of the convolution. We revisit the multi-headed design of H3 [1]; tying weights across channels pushes long convolution filters towards larger effective dimension, and as an additional advantage reduces the runtime of post-training distillation and inference memory footprint. Further, multi-head Hyena models improve on pre-training perplexity over regular Hyena and GPT [10] architectures on the language dataset THE PILE [11]. Figure 1.1: Throughput (in generated tokens) of Transformers, H3 and Hyena models. **LaughingHyena** is a recurrent model distilled from a pre-trained Hyena. Workload involves generating 256 tokens given a prompt of length 512. ## 2 Preliminaries and Related Work We discuss convolutions, state spaces and auto-regressive generation workloads for sequence models. **Convolutions** Let $*$ denote the convolution operator. It is defined as the dual operation to point-wise multiplication under Fourier transform. In signal processing and deep learning alike, one often encounters the causal linear convolution of a filter $h$ (which may extend indefinitely) with an input $u$ of length $L$ : $$(h * u)_t = \sum_{j=0}^t h_{t-j} u_j. \quad (2.1)$$ Generally, $u_t \in \mathbb{R}^D$ where $D$ is the width of the signal – or in deep learning parlance – the number of *channels*. Without loss of generality, we specialize our analysis to *single input single output* layers, i.e. with $D = 1$ . For the input-output relations of type (2.1), we use the terms *convolution layer* and *linear system* interchangeably. Similarly, the function $t \mapsto h_t$ is referred to as both the *filter* and the *impulse response* of a linear system. Existing convolution sequence models can be classified in terms of the parametrization used for their filters. The class of *implicit* convolutions represent the filter as a parametric function $\gamma_\theta : t \mapsto h_t$ . **State-space realization** One option is to select $\gamma_\theta$ as the *impulse response* function of a discrete linear time-invariant system, $$\begin{aligned} x_{t+1} &= \mathbf{A}x_t + \mathbf{B}u_t \\ y_t &= \mathbf{C}x_t + h_0 u_t \end{aligned}, \quad t \mapsto h_t = \begin{cases} h_0 & t = 0 \\ \mathbf{C}\mathbf{A}^{t-1}\mathbf{B} & t > 0 \end{cases} \quad (2.2)$$ with *state* $x_t \in \mathbb{R}^d$ , *input* $u_t \in \mathbb{R}$ , and *output* $y_t \in \mathbb{R}$ . The matrices $\mathbf{A} \in \mathbb{R}^{d \times d}$ , $\mathbf{B} \in \mathbb{R}^{d \times 1}$ , $\mathbf{C} \in \mathbb{R}^{1 \times d}$ , and $h_0 \in \mathbb{R}$ are the learnable parameters of the model while the initial state $x_0$ is usually set to zero suchthat $u \mapsto y$ is a pure convolution. While linear systems (2.2) are the staple of signal processes and control theory, their use as implicit parametrization of convolution filters in deep neural networks have only recently emerged [12, 6]. Other parametrizations [13, 14, 2] select $\gamma_\theta(t)$ as different flavors of implicit representation neural networks [15, 16]. The latter are generally more powerful in terms of the class of filters they can represent and flexibility during training, at the cost of losing a fixed state dimension. ## 2.1 Long Convolution Sequence Models The H-family of convolution sequence models – H3 [1] and Hyena [2] – relies on a combination of long convolutions and data-controlled gating to replace attention with sub-quadratic scaling in sequence length1. We use the deep learning convention of naming different projections as *query* $q$ , *key* $k$ and *value* $v$ . Let $M_q$ and $M_k$ be the $L$ -by- $L$ diagonal matrices whose respective main diagonal entries are the respective entries of length- $L$ sequences $q$ and $k$ . A H-block realizes a surrogate attention matrix with a data-controlled, parameterized decomposition in three terms: $$(q, k, v) \mapsto H(q, k)v, \quad H(q, k) = M_q T_h M_k \quad (2.3)$$ where $T_h \in \mathbb{R}^{L \times L}$ is the Toeplitz matrix constructed from the learnable long convolution filter $h$ , i.e., $T_h = (h_{i-j})_{i,j=0}^{L-1}$ . The $qkv$ -projections are themselves the output of a convolution between the input sequence and three distinct *short* filters. The degrees of freedom in H-block design are the three short filters2 and the *long* filter $h$ . The long filter can be parameterized using an implicit neural representation [2], state-space model [1], or explicit values [17]. The threefold decomposition of the attention operator, allows evaluation of (2.3) in just $\tilde{O}(L) := \mathcal{O}(L \log_2 L)$ time (two convolutions3 and two element-wise products), $y_t = q_t (h * kv)_t$ . The overall operator acts on an input $u$ by constructing a third-order multi-variate polynomial of $u$ whose coefficients are controlled (nonlinearly) by parameters of the block. Figure 2.1: H-block. $T^{(q)}$ , $T^{(k)}$ , $T^{(v)}$ are *short*-convolution operators. ## 2.2 Auto-Regressive Generation A typical workload for sequence models is auto-regressive generation. Given a length- $T$ prompt $u \in \mathbb{R}^T$ , the model is tasked with producing the following $K$ additional outputs – one at a time – for a resulting output sequence $y$ of length $L=T+K$ . **Convolution sequence models** After processing the initial prompt in $\tilde{O}(T)$ time and obtaining a length- $T$ output $u \mapsto y_0, \dots, y_{T-1}$ , a generic convolution layer can *cache* the output sequence and generate any additional outputs using (2.1) auto-regressively, i.e. $y_{t+1} = \sum_{j=0}^t h_{t-j} y_j$ for $t=T-1, \dots, T+K-1$ . It is important to note that auto-regressive generation with generic long convolutions is expensive. It comes with a **quadratic** cost in the number $K$ of tokens to be generated and require storing a cache of length up to $L$ . **Lemma 2.1.** *Generating $K$ tokens with a long convolution layer (2.1) from a length- $T$ prompt has time complexity $\mathcal{O}(T \log_2 T + TK + K^2)$ and requires $\mathcal{O}(L)$ memory.* **State-space models** When the linear system admits a state space realization (2.2), i.e. it is able to switch between convolution and recurrent mode, the cost of auto-regressive generation can be dramatically reduced. The memory footprint is $\mathcal{O}(d)$ : all we need to cache is the state $x_t$ , a $d$ -dimensional vector. With some further machinery that we develop in next section, we can retain $\tilde{O}(T)$ time and $\mathcal{O}(T)$ memory to process the prompt4 and initialize the state $x_{T-1}$ . Each additional generation step only requires $\mathcal{O}(d)$ time. **Lemma 2.2.** *Generating $K$ tokens with a state-space model (2.2) from a length- $T$ prompt has time complexity $\mathcal{O}(T \log_2 T + dK)$ and requires $\mathcal{O}(T + d)$ memory.* Note that long filters $h$ truncated to length $d$ (i.e. $h_t=0$ for $t > d-1$ ) can also be interpreted as $d$ -dimensional SSMs (see Appendix A.7) where the state (a cache) coincides with the last $d$ inputs. 1In this work, we consider second-order Hyena blocks [2] to automatically extend our findings to H3 [1]. 2The short filters are *explicitly* parameterized, see [2]. 3The $qkv$ short convolutions can be evaluated in batch with a single pass. The second convolution is the one with the long filter $h$ and performed via Fast Fourier Transform (FFT), hence the $\mathcal{O}(L)$ complexity. 4In §3.4 we show that multiple pre-filling strategies exist, with different trade-offs in time and memory.**Transformers** Self-attention is certainly less efficient than long convolutions in processing the prompt, coming with a hefty $\mathcal{O}(T^2)$ time complexity. However, Transformers can achieve a similar efficiency in auto-regressive generation by **caching** the sequences of past keys $\{k_t\}$ and values $\{v_t\}$ . Specifically, from $t=T-1$ onward, the new projections $(q_{t+1}, k_{t+1}, v_{t+1})$ are evaluated from the current output $y_t$ , and the new output $y_{t+1}$ can be computed in linear time with two reductions $$y_{t+1} = \frac{\sum_{j=0}^{t+1} \varphi(q_{t+1}k_j)v_j}{\sum_{i=0}^{t+1} \varphi(q_{t+1}k_i)} \quad \text{where } \varphi : \mathbb{R} \rightarrow \mathbb{R} \text{ is usually chosen as } \varphi(x) = e^x.$$ **Lemma 2.3.** *Generating $K$ tokens with self-attention from a length- $T$ prompt has time complexity $\mathcal{O}(T^2+TK+K^2)$ and requires $\mathcal{O}(L)$ memory.* ### 3 The Laughing Hyena Distillery In this section, we introduce our distillation method. We discuss choosing an approximation objective, a parametrization for the approximant and setting a target state dimension. Given any pre-trained LCSM, the objective of the distillation procedure is to convert each pre-trained convolution filter into a distinct state-space model (2.2). This should be achieved with the smallest state dimension $d$ which preserves, up to a certain tolerance, the input-output characteristics of the convolution layer. Formally, given a filter $h$ the **distillation problem** is defined as follows. Given the sequence $h_1, \dots, h_L$ , find a state-space model (2.2) of dimension $d \ll L$ , whose input-output behavior *approximates* the one of the convolution with $h$ over the largest class of input sequences. The choice of approximation metrics and assumptions on the input sequences yield different *distillation objectives*. A *distillation algorithm* constitutes a systematic procedure for optimally choosing the systems matrices with respect to a particular objective. In instances where the original filter $h$ is itself the impulse response of a finite-dimensional state-space model, e.g., when attempting distillation of H3 or S4 [6] filters, the term distillation becomes analogous to *model-order reduction*. Hence, in such cases, the distillation algorithm should yield a state-space representation of a lower order state-dimension. There exist several algebraic solutions to the model reduction problem [18, 19, 20], typically seeking low-rank structures of the state space by inspecting some invariant of the system, e.g. the *Gramians* in *balanced truncation* [19, Ch. 7]. The lower-order system is then obtained as a projection of the system dynamics onto the found subspace where the system retains desired characteristics, e.g., input-output behavior, stability, etc. **Truncated filters** In theory, implicitly parameterized convolution filters can represent arbitrarily long signals. In practice, these filters are trained on a fixed *maximum length* $L$ . At inference time the model can then be evaluated for sequences longer than $L$ . During distillation it is nonetheless reasonable to treat the pre-trained filters as potentially very long (even beyond $L$ ) but *finite* impulse response functions [21, 22, 23, 24]. We show how this choice is supported by empirical evidence displaying how pre-trained filters typically decay to zero in finite time (see Appendix D). **Transfer function representation** An alternative description of the system (2.2) is its *transfer function* $H$ , defined as the $z$ -transform of the impulse response $H(z) = \sum_{t=0}^{\infty} h_t z^{-t}$ for all $z \in \mathbb{C}$ where the sum converges. The transfer function is a *proper rational function* of $z$ $$H(z) = h_0 + \mathbf{C}(z\mathbf{I} - \mathbf{A})^{-1}\mathbf{B} = h_0 + \frac{b_1 z^{-1} + \dots + b_d z^{-d}}{1 + a_1 z^{-1} + \dots + a_d z^{-d}}. \quad (3.1)$$ In the $z$ -domain, the transfer function defines the input-output map as $Y(z) = H(z)U(z)$ . Here, $H(z)$ is defined outside the $\mathbb{C}$ -plane circle of radius $\rho(\mathbf{A})$ , $\mathbb{D}_{\rho(\mathbf{A})} := \{z \in \mathbb{C} : |z| > \rho(\mathbf{A})\}$ where $\rho(\mathbf{A})$ is the spectral radius of $\mathbf{A}$ , i.e. the amplitude of its largest eigenvalue. We can recover all characteristics of a given system equivalently from either its transfer function or state-space representations (see Appendix A.3 for further details and derivations). Notably, the transfer function is an *invariant* of the system: if we apply a change of variables to the state, the transfer function remains unchanged (Lemma A.3). This alone should discourage attempts at modeling filters by learning *dense* state-space matrices $\mathbf{A}, \mathbf{B}, \mathbf{C}$ as such: there are infinitely many equivalent state-space realizations that mapFigure 3.1: The LaughingHyena long convolution sequence model distillation blueprint. to the same system. Starting from coefficients $(a_i)$ and $(b_i)$ of the rational transfer function (3.1), we can compute the impulse response in $\tilde{O}(L)$ time (Lemma A.6). Moreover, we can map back the transfer function to a special state-space realization – the *companion* canonical form – whose recurrence has time complexity $\mathcal{O}(d)$ (Lemma A.7), compared to the $\mathcal{O}(d^2)$ of dense state-space matrices. From Lemmas A.3 and A.7 we can also prove that any stable state-space model can be converted by *canonicalization* into its companion form, and thus can be equipped with an efficient recurrence (Thm. A.8). The distillation problem presents several challenges: 1. 1. **Defining the distillation objective.** A primary decision involves selecting a distillation objective. We are primarily interested in metrics of pure discrepancy between each filter of a pre-trained deep model and its approximator, rather than the expected input-output loss over a distribution of inputs. 2. 2. **Choosing a state-space parametrization.** It is crucial to determine a suitable parametrization of the distilled state-space realization. Once this is decided, the task is to identify the parameters that minimize the distillation desiderata, which can involve challenging optimization problems in itself. 3. 3. **Selecting the target state dimension.** Lastly, a challenge is to estimate the degree to which the model’s order can be reduced. In other words, we must select the target state dimension of the distillation process to identify the right trade-off between efficiency and accuracy. In the following, we address each of these challenges, and provide a comprehensive approach (summarized in Figure 3.1) to distill recurrences from convolution-based architectures. ### 3.1 Data-Free Distillation Objectives We focus on distillation objectives that are independent of the training data and the overall architecture of the neural network under consideration. The distillation loss should be chosen as a pure measure of discrepancy between each convolution filter $h_t$ of the model and their finite-dimensional approximations $\hat{h}_t = \mathbf{C}\mathbf{A}^{t-1}\mathbf{B}$ . This approach ensures that we do not require a full sequential inference pass over the pre-trained model at each step of distillation procedure and the distilled model can be more broadly applied to downstream tasks. This choice is supported by Young’s convolution inequality [25, 26], which indicates that the output approximation error has a bound $\|y - \hat{y}\|_r \leq \|h - \hat{h}\|_q \|u\|_p$ for properly chosen norms5. For maximum numerical stability and freedom of parametrization for the approximants, we favor modern unconstrained gradient-based approaches to then solve the resulting distillation program6. We design distillation algorithms which either match filters in *time domain* minimizing the $\ell_2$ error ( $\|h\|_2 := [\sum_{t \in \mathbb{Z}} |h_t|^2]^{1/2}$ ) or match their transfer functions optimally with respect to the $\mathcal{H}_2$ norm ( $\|H\|_2 := [(1/2\pi) \int_{-\pi}^{\pi} |H(e^{i\omega})|^2 d\omega]^{1/2}$ )7. As the distillation is carried out via gradient methods, $\ell_2$ is a natural candidate. $\mathcal{H}_2$ error minimization can instead be used to uniformly bound the worst-case discrepancy as $\|h - \hat{h}\|_{\infty} \leq \|H - \hat{H}\|_2$ (see Appendix A.2 for further details). ### 3.2 Making Hyena Laugh with Modal Interpolation Our degrees of freedoms to solve the distillation problem are the matrices A, B, and C of the state-space realization, which determine the filter for all $t > 0$ . In distilled SSMs, the passthrough (residual) term 5 $p, q, r > 0$ should satisfy $1/q + 1/p = 1/r + 1$ . In the case of infinite sequences defined on the all $\mathbb{Z}$ , the norms are taken in a $\ell_p, \ell_q, \ell_r$ sense, respectively. The bound is potentially sharp [27, 28] 6For completeness, we also test *balanced* and *modal* truncation techniques on a suite of pre-trained H3 and Hyena models in Appendix E.3. 7Such norms are always well-defined for finite sequences of interest which are in $\ell_{\infty}$ .cannot be freely assigned: it is simply $h_0$ , the value of the original filter at zero. Alternatively, given its appealing invariance properties, we can parametrize a proper rational function $\hat{H}(z)$ (3.1) and fit it to the (truncated) transfer function8 of the original filter $H_L(z) := \sum_{t=0}^L h_t z^{-t}$ (see Appendix B.2). **Modal canonical form** Optimizing the full transfer function can be numerically challenging for several reasons e.g., ensuring stability9, and ill-posedness for high-order polynomials. A natural solution, inspired by barycentric approaches to rational function approximation [29, 8], is to assume $d$ distinct roots $\lambda_n$ in the denominator’s polynomial, $\lambda_n \in \text{roots}(\text{poly}(a))$ . **Proposition 3.1** ([5]). *If $\text{poly}(a)$ has distinct roots $\{\lambda_n \in \mathbb{C}\}$ , then the transfer function of the system can be factorized as $\hat{H}(z) = \sum_{n=1}^d R_n / (z - \lambda_n)$ , $\forall z \in \mathbb{D}_{\rho(A)}$ where $\{R_n \in \mathbb{C}\}$ is the residue associated with the pole $\lambda_n$ .* Computing the inverse transform of the expanded transfer function via, e.g., the *Cauchy residue theorem* [30], shows that the resulting impulse response $\hat{h}$ corresponds to a truncated basis of exponentially decaying complex sinusoids $$\hat{h}_t = \sum_{n=1}^d R_n \lambda_n^{t-1}, \quad R_n, \lambda_n \in \mathbb{C}, t > 0. \quad (3.2)$$ In practice, this corresponds to the impulse response of state-space model with diagonal matrix $A = \text{diag}(\lambda_1, \dots, \lambda_d)$ and such that $B_i C_i = R_i$ for all $i = 1, \dots, d$ . The distillation problem can be then defined in terms of the $L$ -point nonlinear least squares interpolation error (squared $\ell_2$ ) between $h_1, \dots, h_L$ and (3.2) evaluated for $t=1, \dots, L$ : $\min_{\{\lambda_n, R_n\}} \|\hat{h} - h\|_2^2$ . Note that in case of the target filter $h$ being real-valued, the objective can be replaced by $\|\Re[\hat{h}] - h\|_2^2$ . Figure 3.2: Example of modal interpolation. The approximant is a linear combination of exponentially-decaying complex exponential basis functions with learned decay rate. Although we find solutions of the distillation (interpolation) problem via modern gradient-based optimization techniques, it is worth mentioning that Prony showed how the nonlinear least square solution can be computed solving two linear problems [9]. However, similar to Padé’s method for rational approximation [31], these techniques can be numerically unstable. We opt for a parametrization similar to [32, 33] where each eigenvalue is parameterized in polar form $\lambda_n := A_n e^{i\theta_n}$ and the residues in cartesian form10. Note that, with this parametrization we have $\Re[\hat{h}_t] = \sum_n A_n^{t-1} [\Re(R_n) \cos(\theta_n(t-1)) - \Im(R_n) \sin(\theta_n(t-1))]$ . We can also solve the distillation problem in the $\mathcal{H}_2$ sense by evaluating $\hat{h}_t$ and $h_t$ at $t = 0, \dots, L-1$ and taking their respective (discrete) Fourier transform before computing the objective. Efficient evaluation of (3.2) is crucial for distillation. In particular we show the following: **Lemma 3.1.** *Evaluation of $(\hat{h}_t)_{t=0}^{L-1}$ (3.2) can be done in $\mathcal{O}(dL)$ time from its modal form and in $\tilde{\mathcal{O}}(L)$ time from its proper rational form.* ### 3.3 Minimal Distillation Orders Distilling into lower-dimensional systems is always desirable as they require fewer parameters to be optimized and they yield recurrences that are (linearly) more efficient in terms of time and memory complexity in post-distillation auto-regressive inference workloads. The dimension of the *smallest possible state-space model with impulse response exactly $\{h_t\}_{t \in \mathbb{N}}$* is the so-called *McMillan degree* [34]: $$d^* = \arg \min_d : \exists A \in \mathbb{C}^{d \times d}, B \in \mathbb{C}^{d \times 1}, C \in \mathbb{C}^{1 \times d} \text{ with } h_t = CA^{t-1}B, \forall t > 0 \quad (3.3)$$ **Theorem 3.1** (Ho-Kalman [35, Theorem 2, Corollary]). *Let $S$ be the (infinite) Hankel matrix constructed with $h$ , i.e. $S := (h_{i+j})_{i,j=1}^\infty$ . Then, $d^* = \text{rank}(S)$ .* 8As already partially discussed in [6], the truncation introduces a correction term in the approximant transfer function. See Appendix A.4. 9i.e. normalizing denominator polynomial coefficients to constrain roots within the unit circle. 10We report additional details on the nuances of the parametrization in Appendix B.1.A lower bound for $d^*$ can be estimated from a truncated filter of length $L$ by constructing the $L \times L$ principal sub-matrix $S_L$ and using the fact that $\text{rank}(S) \geq \text{rank}(S_L)$ . Inspecting how fast the Hankel singular values $(\sigma_n)_{n=1}^L$ decay in pre-trained convolution models can be predictive of the approximation quality at a fixed dimension. As a rule of thumb, $d$ needs to be sufficiently large for $\sigma_{d+1}$ to be sufficiently small11. Specifically, we can prove that the *last* singular value $\sigma_d$ determines the upper bound of distillation quality with a SSM of dimension $d$ , in terms of the Hankel norm [19]. This is a direct consequence of Adamjan-Arov-Krein theorem [7] and can be informally stated as follows. **Theorem 3.2** (Informal). *Let $h$ be a length- $L$ filter, $\hat{h}$ a distilled filter of order $d < L$ and let $S_L, \hat{S}_L$ be the respective Hankel matrices. Then $\inf_{\hat{S}_L} \|S_L - \hat{S}_L\|_2 = \sigma_d$ .* ### 3.4 Deploying the Recurrence Once all the filters of a pre-trained model have been distilled with the proposed modal interpolation technique described above, the model unlocks a *recurrent mode* which allocates a state $x_t \in \mathbb{C}^d$ for each filter and enables fast auto-regressive inference. Deployment of distilled model involves two critical steps: the *pre-filling* and the recurrent update rule itself. **Fast pre-filling** During auto-regressive generation, when a length- $T$ prompt is fed to the model, we need to compute the state $x_T$ to start generating new tokens. Using the recurrence, the time complexity of initializing $x_T$ would be $\mathcal{O}(dT)$ with a $\mathcal{O}(d)$ memory footprint. One can alternatively distribute the computation on $d$ processors with a *parallel scan* operation [37, 38] to reach a parallel time complexity $\mathcal{O}(d \log_2 T)$ while incurring an increased memory requirement of $\mathcal{O}(dT)$ 12. A third option is to use a single FFT convolution to obtain $x_T$ in $\tilde{\mathcal{O}}(T)$ time and $\mathcal{O}(T)$ memory. **Proposition 3.2.** *$x_T = (\nu_T, \dots, \nu_{T-d})$ where $\nu_t = (g * u)_t$ and $g$ is the filter whose transfer function is $1/\text{den}(\hat{H})(z)$ and can be evaluated in $\tilde{\mathcal{O}}(T)$ .* Note that, the fast pre-filling algorithm established by this result requires evaluating the denominator polynomial of $\hat{H}$ from its roots before deployment. This is equivalent to converting the transfer function from its factorized representation to its rational form (3.1). **Recurrent step** The update rule is diagonal, thus efficiently evaluated in $\mathcal{O}(d)$ time and memory: **Proposition 3.3.** *The filter (3.2) has a state space matrices $A = \text{diag}(\lambda_1, \dots, \lambda_d) \in \mathbb{C}^{d \times d}$ , $B = (1, \dots, 1)^\top \in \mathbb{C}^{d \times 1}$ , $C = (R_1, \dots, R_d) \in \mathbb{C}^{1 \times d}$ whose step can be evaluated in $\mathcal{O}(d)$ time and memory.* As we generally want the output $y_t$ to be real-valued, we can simply update the complex state $x_{t+1} = Ax_t + Bu_t$ and then take the real part of the output, $y_t = \Re[Cx_t] + h_0 u_t$ . ## 4 Multi-head Long Convolutions We can leverage the Hankel spectrum analysis discussed in Section 3.3 to study the dynamics of the effective dimensionality of each convolution filter during LCSMs pre-training. We find that, at initialization, filters correspond to high-dimensional SSMs, and gradually converge to lower-dimensional representations during training. See Appendix E.2 for examples on Hyena and H3 models. This observation leads to the question: *is it advantageous to perform independent long convolutions on each channel, or can we reduce the total number of filters without loss in quality?* To answer this, we adapt the multi-head layer design proposed by H3 [1] to Hyena [2]: 1. 1. Given the projections $q, k, v \in \mathbb{R}^{L \times D}$ , we split them into $M$ chunks of size $N = D/M$ , $q^m, k^m, v^m \in \mathbb{R}^{L \times N}$ . 2. 2. Each chunk is processed by a modified Hyena operator: first, we perform the outer product of $k^m$ and $v^m$ along the spatial dimension, $z^m := k^m \otimes v^m \in \mathbb{R}^{L \times N \times N}$ , apply a long convolution with filter $h^m$ to all $N \times N$ elements independently, then compute $y_t^m = (h^m * z^m)_t q_t^m$ , $y^m \in \mathbb{R}^{L \times N}$ as shown in Figure 4. 3. 3. Finally, we compose $y^1, \dots, y^m$ into a single output $y \in \mathbb{R}^{L \times D}$ via concatenation. 11Formally, this is related to low-rank approximation characteristics of the Hankel operator; rigorous bounds can be constructed by application of the Eckart–Young–Mirsky theorem [36]. 12This strategy can also be used to evaluate the filter $\hat{h}$ alternatively to the standard $\mathcal{O}(dL)$ methodAn instance of a MultiHyena is equipped with $M < D$ distinct long convolution filters, which leads to (a) faster distillation, with less filters to approximate, (b) lower memory footprint, via a total reduction of the states to cache during generation and (c) faster filter generation, by tying the weights of filter parameters. We note that tying weights of key-value projections has also been shown to be an effective technique to reduce memory cost in Transformers [39, 40]. Figure 4.1: A single head of a multi-head Hyena. Crucially, the multi-head structure of MultiHyena enables us to prove favorable scaling in the *associative recall* synthetic task, which was shown in [2] to be predictive of performance at scale. In associative recall, the model is given a sequence of key-value pairs and a query, and is tasked with matching the query to a key in the sequence by returning its associated value. The difficulty of the task grows with the vocabulary size $s$ : larger vocabularies necessitate wider models. **Theorem 4.1.** *The MultiHyena layer, with $\mathcal{O}(\log s)$ heads and model size $\mathcal{O}(\sqrt{s} \log s)$ can solve the associative recall problem, where $s$ denotes the vocabulary size.* In Appendix E.1, we empirically verify improved scaling in vocabulary size with multiple heads. ## 5 Experiments - • **Pretraining:** We pretrain a suite of MultiHyena language models on The Pile [11], investigating scaling of perplexity with different amounts of total tokens (5, 10, 15 billion), as well as larger training runs for 300 billion tokens. MultiHyena outperforms Transformers and Hyena. - • **Distillation analysis:** We investigate the relation between optimal distillation orders, Hankel spectrum, and errors on the logits of distilled models. - • **Post-distillation downstreams:** We evaluate the downstream impact of distilling long convolutional language models, reporting HELM [41] and LM-Eval-Harness [42] results. - • **Benchmarking:** We benchmark latency, throughput and memory along the different axes of batch size, sequence length, number of generated tokens. We include base models, distilled models and equivalent Transformers. ### 5.1 Pre-training To validate the multi-head formulation, we train 150 and 350 million parameter MultiHyena models on The Pile [11] using 8 heads and otherwise the same architecture as equivalent Hyena models, following the setup of [2]. Via the multi-head structure introduced in 4, MultiHyena outperforms both Hyena and Transformers, including on data scaling runs with increasing numbers of tokens and full 300B tokens runs (Table 5.1). ### 5.2 Distillation Analysis Next, we verify whether Hankel singular values are predictive of downstream errors, and whether large models can be distilled without loss in quality. We apply LaughingHyena distillation to pre-trained MultiHyena, Hyena and H3 of different sizes. Concretely, for each layer and channel of a model, we parametrize the poles $\{\lambda_n\}$ of the modal canonical forms (Section 3.2) at different orders $d$ , and solve for each $\ell_2$ approximation problem. **Approximation errors and spectrum** We investigate the magnitude of approximation errors introduced by LaughingHyena distillation. Given a pretrained MultiHyena model, we compute the errors between original and distilled filters at each layer, averaged across channels. We repeat this process for different distillation orders (state dimension of the model form of Section 3.2). Figure 5.2 visualizes minimum, maximum and average errors, per-layer errors and the distribution of the singular values of the Hankel operator
Model PPL. Model 5B 10B 15B Model 5B 10B 15B
GPT 9.3 GPT (125M) 13.3 11.9 11.2 GPT (355M) 11.4 9.8 9.1
Hyena 9.3 Hyena (153M) 13.3 11.8 11.1 Hyena (355M) 11.3 9.8 9.2
MultiHyena 8.7 MultiHyena (153M) 12.1 11.0 10.6 MultiHyena (355M) 10.6 9.4 8.9
Table 5.1: [Left] Perplexity of small models on THE PILE, after pre-training for 300 billion tokens. [Center and Right] Perplexity on THE PILE for models trained until a total number of tokens e.g., 5 billion (different runs for each token total).Figure 5.1: Errors between logits of pretrained and distilled MultiHyena. In blue, we plot (ordered) logits, in light blue the cumulative distribution function, and in black the relative errors. The green dotted line indicates the 99.99% percentile. As the errors grow slowly as function of the percentiles, model outputs do not diverge from the base model. associated to each filter. We observe distillation orders ( $> 16$ ) that yield small errors to be predicted by the distribution of singular values. Thus, analysis of the Hankel operator’s spectrum is verified to be an effective approach to direct estimation of the optimal distillation order. We also note that the optimal order changes across layers, offering options for further optimization. **Output errors** Next, we compute relative $\ell_1$ error between output logits of pre-trained and distilled models to ensure LaughingHyena can be used in generation workloads. The optimal minimal distillation order estimated via Hankel operators (16) is sufficient to keep the output distribution over the vocabulary ( $> 50k$ entries) close to the pre-trained model, as shown in Figure 5.2. Inspecting the error profile over logits sorted by magnitude reveals our approach to be robust to different sampling strategies for generation, including greedy decoding, top- $k$ , top- $p$ [43]. Indeed, the relative errors are $< 10^{-2}$ up to and including the 99.99% percentile of the distribution, meaning e.g., a top- $p$ sampling strategy with large $p$ can be used on a distilled model without drift in outputs (mis-classified tokens). We note that the relative errors are maximum on small-norm logits, which are not required by most sampling strategies. In Appendix D.2, we provide a similar distillation error analysis for Hyena and H3 models. We find that Hyena can be distilled with less than 32 orders and H3 with less than 8. ### 5.3 Downstream Evaluation We check how distillation affects downstream performance on language benchmarks. We apply distillation of order 8, 16 and 32 to our THE PILE-pretrained MultiHyena language model and benchmark (Table 5.3) its performance on a suite of canonical (zero shot) tasks from LM-Eval-Harness [42] and HELM [41]. The results are consistent with our error analysis: distillation orders equal or greater to 16 introduce little-to-no quality degradation.
Model LAMBADA Winogrande PIQA HellaSwag OpenbookQA
acc acc acc acc norm. acc norm.
Pythia (160M) 32.8 53.1 61.6 31.6 29.2
MultiHyena (154M) 43.2 52.7 64.6 34.1 29.0
LaughingHyena-16 43.1 52.6 64.7 34.1 28.9
LaughingHyena-8 0.0 51.8 51.5 32.7 28.2
LaughingHyena-4 0.0 49.6 53.7 26.4 26.4
Table 5.2: Evaluation of LaughingHyena-distilled models pre and post modal distillation. We test on LM-Eval-Harness tasks, reporting Pythia [44] performance as a Transformer baseline trained on the same data. LaughingHyena- $d$ is a MultiHyena model with each filter distilled of order $d$ .Figure 5.2: Approximation error profiles (min, max, average) on the filters of MultiHyena model after distillation at different orders. We also visualize the distribution of Hankel singular values: if the spectrum decays after $n$ singular values, order $n$ distillation yields low errors. ## 5.4 Benchmarking We measure throughput, latency and memory usage of **LaughingHyena** for auto-regressive generation workloads, with initial prompt length $T$ and number of generated tokens $K$ . The throughput is computed as number of generated tokens over latency. For each setting (and additional benchmarks), we provide details in Appendix D.4. **Peak throughput** Distilled models do not need $kv$ -caches. This reduces memory requirement during generation, enabling higher peak throughput in large-batch workloads. We achieve $10\times$ higher throughput than Transformers at size 1.3 billion parameters (Figure 1.1). Throughput is higher than Transformers even at fixed batch sizes, indicating lower latency. **SSM state dimension and throughput** For typical distillation orders ( $< 100$ ), peak throughput is not greatly affected. We measure a 2% reduction in throughput from 32 to 64. **Prompt length** The throughput of **LaughingHyena**-distilled models is $4\times$ larger than Transformers at fixed batch size 64 and prompt length 1536 (Figure 5.3). As prompt length increases, the runtime gap between pre-filling via convolutions in LCSMs and pre-filling in Transformers widens (e.g., $\tilde{\mathcal{O}}(T)$ as detailed in Section 3.4, compared to $\mathcal{O}(T^2)$ ). **Memory footprint** Recurrent models do not require $kv$ -caches and use constant memory for generation of an arbitrary number of tokens (Figure 5.4). ## 6 Conclusion We study the efficiency and quality of state-of-the-art long convolutional sequence models. First, we introduce **LaughingHyena**, a novel distillation method inspired by rational function approximation and model-order reduction techniques. **LaughingHyena** can be applied after training to extract compact state-space models from each convolutional filter, without loss of quality. Distilled models achieve higher throughput than equivalently-sized Transformers, and can perform auto-regressive generation in constant memory by sidestepping the need to cache previous outputs. We theoretically and empirically investigate the trade-offs of different strategies for fast inference of recurrent models, and introduce architectural improvements to Hyena that improve pretraining quality. Figure 5.3: Scaling in prompt length $T$ . Figure 5.4: Peak GPU memory for generation.## Acknowledgments We would like to thank Together Computer for providing the compute used to train models in this paper. We gratefully acknowledge the support of NIH under No. U54EB020405 (Mobilize), NSF under Nos. CCF1763315 (Beyond Sparsity), CCF1563078 (Volume to Velocity), and 1937301 (RTML); US DEVCOM ARL under No. W911NF-21-2-0251 (Interactive Human-AI Teaming); ONR under No. N000141712266 (Unifying Weak Supervision); ONR N00014-20-1-2480: Understanding and Applying Non-Euclidean Geometry in Machine Learning; N000142012275 (NEPTUNE); NXP, Xilinx, LETI-CEA, Intel, IBM, Microsoft, NEC, Toshiba, TSMC, ARM, Hitachi, BASF, Accenture, Ericsson, Qualcomm, Analog Devices, Google Cloud, Salesforce, Total, the HAI-GCP Cloud Credits for Research program, the Stanford Data Science Initiative (SDSI), Department of Defense (DoD) through the National Defense Science and Engineering Graduate Fellowship (NDSEG) Program, and members of the Stanford DAWN project: Facebook, Google, and VMWare. This work is supported by NSF (1651565), AFOSR (FA95501910024), ARO (W911NF-21-1-0125), ONR, DOE (DE-SC0022222), CZ Biohub, and Sloan Fellowship. 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 NIH, ONR, or the U.S. Government. AR’s work is supported by NSF grant# CCF-2247014. ## Broader Impact In this work, we focus on advances related to efficient models for long sequences. **Efficiency** Our distillation methods for constant-memory, high throughput inference in *long convolution sequence models* (LCSMs) can lead to energy savings during model deployment, enabling processing of longer-form content at a fraction of the cost and reducing environmental impact. Improved efficiency may also affect other aspects of AI safety, as it may make it easier produce malicious or harmful content. **Accessibility** By improving the efficiency of training and generation, LCSMs and LaughingHyena may contribute to increased accessibility of large language models, lowering the hardware barrier to entry for individuals and organizations with limited resources. **Steerability** New method based on LCSMs enable sequence models to process long-form prompts previously inaccessible by Transformers, which may lead to increased control over models via e.g., conditioning on additional instructions [\[45\]](#).## References - [1] Daniel Y Fu et al. “Hungry Hungry Hippos: Towards Language Modeling with State Space Models”. In: 2023 (cit. on pp. 1–3, 7, 23). - [2] Michael Poli et al. “Hyena Hierarchy: Towards Larger Convolutional Language Models”. In: (2023). arXiv: [2302.10866](#) (cit. on pp. 1, 3, 7, 8, 31, 35, 43). - [3] Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. “Neural machine translation by jointly learning to align and translate”. In: (2014). arXiv: [1409.0473](#) (cit. on p. 1). - [4] Ashish Vaswani et al. “Attention is all you need”. In: *Advances in neural information processing systems* 30 (2017) (cit. on p. 1). - [5] Chi-Tsong Chen. *Linear system theory and design*. Saunders college publishing, 1984 (cit. on pp. 1, 6, 19, 21, 22). - [6] Albert Gu, Karan Goel, and Christopher Ré. “Efficiently modeling long sequences with structured state spaces”. In: (2021). arXiv: [2111.00396](#) (cit. on pp. 1, 3, 4, 6, 20, 25). - [7] Vadim Movsesovich Adamyan, Damir Zyamovich Arov, and Mark Grigor’evich Krein. “Analytic properties of Schmidt pairs for a Hankel operator and the generalized Schur–Takagi problem”. In: *Matemathikeskii Sbornik* 128.1 (1971), pp. 34–75 (cit. on pp. 2, 7, 29). - [8] Yuji Nakatsukasa, Olivier Sète, and Lloyd N Trefethen. “The AAA algorithm for rational approximation”. In: *SIAM Journal on Scientific Computing* 40.3 (2018), A1494–A1522 (cit. on pp. 2, 6). - [9] GRB Prony. “Essai experimental et analytique sur les lois de la dilatalrlite de fluids elastiques et sur cells de la vapeur de l’alcool, à différents tempoeratures”. In: *Journal de l’Ecole Polytechnique (Paris)* 1 (1795), pp. 24–76 (cit. on pp. 2, 6). - [10] Alec Radford et al. “Language models are unsupervised multitask learners”. In: *OpenAI blog* 1.8 (2019), p. 9 (cit. on p. 2). - [11] Leo Gao et al. “The pile: An 800gb dataset of diverse text for language modeling”. In: (2020). arXiv: [2101.00027](#) (cit. on pp. 2, 8, 35). - [12] Albert Gu et al. “Hippo: Recurrent memory with optimal polynomial projections”. In: *Advances in Neural Information Processing Systems* 33 (2020), pp. 1474–1487 (cit. on pp. 3, 25). - [13] David W Romero et al. “Ckconv: Continuous kernel convolution for sequential data”. In: (2021). arXiv: [2102.02611](#) (cit. on p. 3). - [14] David W Romero et al. “Flexconv: Continuous kernel convolutions with differentiable kernel sizes”. In: (2021). arXiv: [2110.08059](#) (cit. on p. 3). - [15] Vincent Sitzmann et al. “Implicit neural representations with periodic activation functions”. In: *Advances in neural information processing systems* 33 (2020), pp. 7462–7473 (cit. on p. 3). - [16] Rizal Fathony et al. “Multiplicative filter networks”. In: *International Conference on Learning Representations*. 2020 (cit. on p. 3). - [17] Daniel Y. Fu et al. “Simple Hardware-Efficient Long Convolutions for Sequence Modeling”. In: *International Conference on Machine Learning* (2023) (cit. on pp. 3, 25). - [18] Kemin Zhou and John Comstock Doyle. *Essentials of robust control*. Vol. 104. Prentice hall Upper Saddle River, NJ, 1998 (cit. on p. 4). - [19] Athanasios C Antoulas. *Approximation of large-scale dynamical systems*. SIAM, 2005 (cit. on pp. 4, 7, 29). - [20] Wilhelmus HA Schilders, Henk A Van der Vorst, and Joost Rommes. *Model order reduction: theory, research aspects and applications*. Vol. 13. Springer, 2008 (cit. on p. 4). - [21] Sun-Yuan Kung. “A new identification and model reduction algorithm via singular value decomposition”. In: *Proc. 12th Asilomar Conf. on Circuits, Systems and Computer*. 1978, pp. 705–714 (cit. on p. 4). - [22] D Friedman. “On approximating an FIR filter using discrete orthonormal exponentials”. In: *IEEE Transactions on Acoustics, Speech, and Signal Processing* 29.4 (1981), pp. 923–926 (cit. on p. 4). - [23] J Bednar. “On the approximation of FIR by IIR digital filters”. In: *IEEE Transactions on Acoustics, Speech, and Signal Processing* 31.1 (1983), pp. 28–34 (cit. on p. 4).- [24] Bartłomiej Beliczynski, Izzet Kale, and Gerald D Cain. “Approximation of FIR by IIR digital filters: An algorithm based on balanced model reduction”. In: *IEEE Transactions on Signal Processing* 40.3 (1992), pp. 532–542 (cit. on pp. 4, 43, 44). - [25] William Henry Young. “On the multiplication of successions of Fourier constants”. In: *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character* 87.596 (1912), pp. 331–339 (cit. on p. 5). - [26] William Beckner. “Inequalities in Fourier analysis on $\mathbb{R}^n$ ”. In: *Proceedings of the National Academy of Sciences* 72.2 (1975), pp. 638–641 (cit. on p. 5). - [27] John Fournier. “Sharpness in Young’s inequality for convolution”. In: *Pacific Journal of Mathematics* 72.2 (1977), pp. 383–397 (cit. on p. 5). - [28] Tong S Quek and Leonard YH Yap. “Sharpness of Young’s inequality for convolution”. In: *Mathematica Scandinavica* 53.2 (1983), pp. 221–237 (cit. on p. 5). - [29] Jean-Paul Berrut and Lloyd N Trefethen. “Barycentric lagrange interpolation”. In: *SIAM review* 46.3 (2004), pp. 501–517 (cit. on p. 6). - [30] Marcos Vicente Moreira and João Carlos Basilio. “Fair and Square Computation of Inverse Z-Transforms of Rational Functions”. In: *IEEE Transactions on Education* 55.2 (2011), pp. 285–290 (cit. on p. 6). - [31] Henri Padé. “Sur la représentation approchée d’une fonction par des fractions rationnelles”. In: *Annales scientifiques de l’Ecole normale supérieure*. Vol. 9. 1892, pp. 3–93 (cit. on pp. 6, 27). - [32] Ankit Gupta, Albert Gu, and Jonathan Berant. “Diagonal state spaces are as effective as structured state spaces”. In: *Advances in Neural Information Processing Systems* 35 (2022), pp. 22982–22994 (cit. on pp. 6, 25). - [33] Antonio Orvieto et al. “Resurrecting Recurrent Neural Networks for Long Sequences”. In: (2023). arXiv: [2303.06349](#) (cit. on pp. 6, 25, 26). - [34] Jeffrey M Hokanson. “A data-driven McMillan degree lower bound”. In: *SIAM Journal on Scientific Computing* 42.5 (2020), A3447–A3461 (cit. on p. 6). - [35] L Ho and Rudolf E Kalman. “Effective construction of linear state-variable models from input/output functions”. In: *at-Automatisierungstechnik* 14.1-12 (1966), pp. 545–548 (cit. on p. 6). - [36] Carl Eckart and Gale Young. “The approximation of one matrix by another of lower rank”. In: *Psychometrika* 1.3 (1936), pp. 211–218 (cit. on p. 7). - [37] Guy E Blelloch. “Prefix sums and their applications”. In: (1990) (cit. on p. 7). - [38] Jimmy TH Smith, Andrew Warrington, and Scott W Linderman. “Simplified state space layers for sequence modeling”. In: (2022). arXiv: [2208.04933](#) (cit. on p. 7). - [39] Noam Shazeer. “Fast transformer decoding: One write-head is all you need”. In: (2019). arXiv: [1911.02150](#) (cit. on p. 8). - [40] Joshua Ainslie et al. “GQA: Training Generalized Multi-Query Transformer Models from Multi-Head Checkpoints”. In: (2023). arXiv: [2305.13245](#) (cit. on p. 8). - [41] Percy Liang et al. “Holistic evaluation of language models”. In: (2022). arXiv: [2211.09110](#) (cit. on pp. 8, 9, 36). - [42] Leo Gao et al. *A framework for few-shot language model evaluation*. Version v0.0.1. Sept. 2021. DOI: [10.5281/zenodo.5371628](#). URL: (cit. on pp. 8, 9, 36). - [43] Ari Holtzman et al. “The curious case of neural text degeneration”. In: (2019). arXiv: [1904.09751](#) (cit. on p. 9). - [44] Stella Biderman et al. “Pythia: A suite for analyzing large language models across training and scaling”. In: (2023). arXiv: [2304.01373](#) (cit. on pp. 9, 36). - [45] Yuntao Bai et al. “Constitutional AI: Harmlessness from AI Feedback”. In: (2022). arXiv: [2212.08073](#) (cit. on p. 11). - [46] Alan V Oppenheim. *Discrete-time signal processing*. Pearson Education India, 1999 (cit. on p. 21). - [47] Lennart Ljung. *System identification*. Springer, 1998 (cit. on p. 21). - [48] RP Guidorzi. “Certain models from uncertain data: the algebraic case”. In: *Systems & control letters* 17.6 (1991), pp. 415–424 (cit. on p. 21). - [49] Irwin W Sandberg. “On the theory of linear multi-loop feedback systems”. In: *Bell System Technical Journal* 42.2 (1963), pp. 355–382 (cit. on p. 23).- [50] Michael Zhang et al. “Effectively Modeling Time Series with Simple Discrete State Spaces”. In: (2023). arXiv: [2303.09489](#) (cit. on p. [25](#)). - [51] Luca Perotti and Michał Wojtylak. “Matrix methods for Padé approximation: Numerical calculation of poles, zeros and residues”. In: *Linear Algebra and its Applications* 548 (2018), pp. 95–122 (cit. on p. [27](#)). - [52] Amer Abu-Omar and Fuad Kittaneh. “Estimates for the numerical radius and the spectral radius of the Frobenius companion matrix and bounds for the zeros of polynomials”. In: *Annals of Functional Analysis* 5.1 (2014), pp. 56–62 (cit. on p. [27](#)). - [53] Julia Eaton et al. “Polynomial root radius optimization with affine constraints”. In: *Mathematical Programming* 165 (2017), pp. 509–528 (cit. on p. [27](#)). - [54] Gerlind Plonka and Vlada Pototskaia. “Application of the AAK theory for sparse approximation of exponential sums”. In: (2016). arXiv: [1609.09603](#) (cit. on p. [29](#)). - [55] Jimmy Ba et al. “Using fast weights to attend to the recent past”. In: *Advances in neural information processing systems* 29 (2016) (cit. on p. [31](#)). - [56] William B Johnson. “Extensions of Lipschitz mappings into a Hilbert space”. In: *Contemp. Math.* 26 (1984), pp. 189–206 (cit. on p. [31](#)). - [57] Devdatt P Dubhashi and Alessandro Panconesi. *Concentration of measure for the analysis of randomized algorithms*. Cambridge University Press, 2009 (cit. on p. [34](#)). - [58] Ilya Loshchilov and Frank Hutter. “Decoupled weight decay regularization”. In: (2017). arXiv: [1711.05101](#) (cit. on p. [35](#)). - [59] Dale F. Enns. “Model reduction with balanced realizations: An error bound and a frequency weighted generalization”. In: *The 23rd IEEE Conference on Decision and Control* (1984), pp. 127–132 (cit. on p. [43](#)).--- # LAUGHING HYENA DISTILLERY ## *Supplementary Material* --- ### Contents
1Introduction1
2Preliminaries and Related Work2
2.1Long Convolution Sequence Models . . . . .3
2.2Auto-Regressive Generation . . . . .3
3The Laughing Hyena Distillery4
3.1Data-Free Distillation Objectives . . . . .5
3.2Making Hyena Laugh with Modal Interpolation . . . . .5
3.3Minimal Distillation Orders . . . . .6
3.4Deploying the Recurrence . . . . .7
4Multi-head Long Convolutions7
5Experiments8
5.1Pre-training . . . . .8
5.2Distillation Analysis . . . . .8
5.3Downstream Evaluation . . . . .9
5.4Benchmarking . . . . .10
6Conclusion10
ALinear Systems17
A.1Extended Notation and System Theory Preliminaries . . . . .17
A.2Systems Norms . . . . .17
A.3Transfer Function of State-Space Models . . . . .18
A.4Truncated Transfer Functions . . . . .19
A.5From Transfer Function to State-Space . . . . .21
A.5.1Isolating the h_0-term from Transfer Function by Long division . . . . .21
A.5.2Construction of the State-Space from the Transfer Function . . . . .22
A.6From State-Space to Transfer Function . . . . .23
A.7State-Space Representation of Truncated Filters. . . . .23
A.8Efficient Computation of State-Space Models . . . . .24
A.8.1Fast Evaluation of the Transfer Function . . . . .24
A.8.2Fast Companion Recurrence . . . . .24
A.8.3Canonization of State-Space Models . . . . .25
BLaughingHyena: Further Details26
B.1Parametrization of Modal Interpolators . . . . .26
B.2Distillation as Rational Interpolation . . . . .27
CProofs28
C.1Proof of Lemma 2.1 . . . . .28
C.2Proof of Lemma 2.2 . . . . .28
C.3Proof of Lemma 2.3 . . . . .28
C.4Proof of Proposition 3.1 . . . . .29
C.5Proof of Lemma 3.1 . . . . .29
C.6Proof of Theorem 3.2 . . . . .29
C.7Proof of Proposition 3.3 . . . . .30
C.8 Proof of Proposition 3.2 30
C.9 Proof of Theorem 4.1 31
D Experimental Details 35
D.1 Pre-training 35
D.2 Distillation Analysis 35
D.2.1 Pretrained Filters: Effective Dimension 35
D.3 Downstream Evaluation 36
D.4 Benchmarking 36
E Additional Experiments 43
E.1 Associative Recall with MultiHyena 43
E.2 Analysis of Hankel Singular Values of Pretrained Large Convolution Sequence Models 43
E.3 Model Order Reduction of H3 43
E.3.1 Modal Truncation 43
E.3.2 Balanced Truncation 43
## Authors Contribution
S.M. Conceptualized the research; coordinated collaborations; lead theory development; conducted distillation experiments.
M.P. Conceptualized the research; coordinated collaborations; lead the experimental (model pre-training, distillation, benchmarks, downstream evaluation) efforts; coordinated writing and conference submission; optimized inference stack.
D.Y.F. Assisted in development of MultiHyena; assisted in pre-training and subsequent benchmarking of distilled models; assisted in writing.
H.K. Developed benchmarking suite and interpreted results; assisted in writing.
R.N.P. Assisted in theory and algorithmic development; performed model-order reduction experiments of H3 models; assisted in writing.
A.T. Conceived and proved Theorem 4.1; assisted in writing.
D.W.R. Assisted in Hankel operator spectral analysis; Assisted in writing.
Q.M. Assisted in theory development.
B.C. Supervised development of benchmarking suite and model deployment.
A.R. Supervised theory development (solving associative recall with MultiHyena, Th. 4.1).
C.Z. Supervised research; secured compute resources.
C.R. Supervised research; reviewed manuscript; secured compute resources.
S.E. Supervised research; reviewed manuscript.
Y.B. Supervised research; reviewed manuscript.
*Stefano Massaroli, Michael Poli, and Dan Fu contributed equally to this work. Christopher Ré, Stefano Ermon, and Yoshua Bengio share equal senior authorship.* All authors read and approved the final manuscript.## A Linear Systems ### A.1 Extended Notation and System Theory Preliminaries We first introduce the notation and some mathematical concepts that will be used throughout the paper. By $\mathbb{Z}$ we denote the set of integers, by $\mathbb{R}$ the set of reals, and by $\mathbb{C}$ the set of complex numbers. The variable $t$ stands for *time*. $\ell_p(\mathbb{Z})$ denotes the Banach space of complex-valued sequences $(x_t)_{t \in \mathbb{Z}}$ with finite energy, i.e. $\|x\|_p := [\sum_{t \in \mathbb{Z}} |x_t|^p]^{1/p} < \infty$ for some $1 \leq p < \infty$ . $\ell_\infty(\mathbb{Z})$ is instead is the space of sequences for which $\|x\|_\infty := \sup_{t \in \mathbb{Z}} |x_t| < \infty$ . With $\mathbb{S}$ denoting the unit circle in the complex plane, $\mathbb{S} := \{z \in \mathbb{C} : |z| = 1\}$ we define $\mathcal{H}_p(\mathbb{S})$ as the space of functions $X$ from $\mathbb{C}$ to itself such that $\|X\|_p := [(1/2\pi) \int_{-\pi}^{\pi} |X(e^{i\omega})|^p d\omega]^{1/p} < \infty$ and $\mathcal{H}_\infty(\mathbb{S})$ the space for which $\|X\|_\infty := \sup_{z \in \mathbb{S}} |X(z)| < \infty$ . Particularly, $\mathcal{K}_2(\mathbb{S})$ is a Hilbert space with inner product $\langle X, Y \rangle := (1/2\pi) \int_{-\pi}^{\pi} X(e^{i\omega}) Y^*(e^{i\omega}) d\omega$ where “\*” denotes complex conjugation. Although we acknowledge we are using the same notation for norms in both $\ell_p(\mathbb{Z})$ and $\mathcal{H}_p(\mathbb{S})$ , the correct meaning will always be made clear by the context. The $\mathcal{Z}$ -transform of a sequence $x = (x_t)_{t \in \mathbb{Z}}$ is $X(z) = \mathcal{Z}[x](z) := \sum_{t \in \mathbb{Z}} x_t z^{-t}$ . We embrace the system theory convention of using capital letters to identify transformed sequences. The $\mathcal{Z}$ -transform is a projection of the sequence onto a basis of powers $e_t = r^{-t} e^{i\omega t}$ . This basis is not orthogonal unless $r = 1$ . That is the basis of the discrete-time Fourier transform $\mathcal{F}$ . Hence, $\mathcal{F}$ is defined as $\mathcal{F}[x](e^{i\omega}) = X(e^{i\omega}) := \sum_{t \in \mathbb{Z}} x_t e^{-i\omega t}$ . The discrete-time Fourier transform is an isometric isomorphism between $\ell_2(\mathbb{Z})$ and $L_2(\mathbb{S})$ . We say that sequences live in the *time domain* and their $\mathcal{Z}$ (or $\mathcal{F}$ ) transforms in the *frequency domain*. A *linear system* is a linear operator transforming an input sequence $u$ to an output sequence $y$ . If the sequences have continuous support, i.e. $t$ ranges over a continuous set (e.g. $\mathbb{R}$ ), we have a *continuous-time* system. Conversely, if the sequences have discrete support, i.e. $t$ ranges over a discrete set (e.g. $\mathbb{Z}$ ), we have a *discrete-time* or *digital* system. **In this manuscript we restrict ourselves to discrete-time systems.** Systems can be *single-input single-output* (SISO) if $u$ and $y$ are scalar functions or *multi-input multi-output* if either $u$ or $y$ are vector-valued. **We limit our discussion to SISO systems.** The *impulse response* of a system is the output sequence $y$ when the input sequence $u$ is the Kronecker delta function $\delta_t$ and is usually denoted by the letter $h$ . The values $h_t$ of the impulse response sequence are also known as the *Markov parameters* of the system. The most common mathematical representation of a linear system is its convolution form: $y = h * u$ , i.e. $y_t = \sum_{j \in \mathbb{Z}} h_{t-j} u_j = \sum_{j \in \mathbb{Z}} h_j u_{t-j}$ , $t \in \mathbb{Z}$ . In matrix form the input-output relation is given by the Toeplitz operator $\mathsf{T}_h$ corresponding to the (possibly infinitely long) sequence $h$ , i.e. $y = \mathsf{T}_h u$ . Taking the $\mathcal{Z}$ -transforms, we can write the input-output relation as $Y(z) = H(z)U(z)$ (this is just the Fourier convolution theorem extended outside the unit circle). $H(z)$ is called the *transfer function* of the system. When $z = e^{i\omega}$ , $H(e^{i\omega})$ is just the discrete-time Fourier transform of $h$ which is called the *frequency response* of the system. A linear system is *causal* if $h_t = 0$ for $t < 0$ . A system is called *stable* if the $\mathsf{T}_h$ is a bounded operator. If $u, y \in \ell_2(\mathbb{Z})$ , then stability implies $h \in \ell_\infty$ . **In the following, we mainly focus on causal stable systems.** ### A.2 Systems Norms When quantitatively characterizing linear systems, several norms play a crucial role. These norms provide measures of various characteristics of the systems, which are essential in both analysis and filter design. **The $\ell_2$ and $\mathcal{H}_2$ norms** As defined above, the $\ell_2$ norm represents the *energy* of a signal $h$ , $$\|h\|_2 := \left[ \sum_{t \in \mathbb{Z}} |h_t|^2 \right]^{1/2}$$ while $\mathcal{H}_2$ is the energy of the (continuous) spectrum of $h$ , $$\|H\|_2 := \left[ \frac{1}{2\pi} \int_{-\pi}^{\pi} |X(e^{i\omega})|^2 d\omega \right]^{1/2}$$ By Parseval’s theorem, the $\ell_2$ and $\mathcal{H}_2$ norms are equal, $\|h\|_2 = \|H\|_2$ . Further these norms are useful to study the approximation of convolutional filter. The following holds: **Lemma A.1** ( $\ell_\infty$ output error). *Consider the class of $\ell_2$ measurable inputs such that $\|u\|_2 \leq \zeta$ , then for all $H, \hat{H} \in \mathcal{H}_2$ ,* $$\|y - \hat{y}\|_\infty \leq \zeta \|H - \hat{H}\|_2$$*Proof.* $$\begin{aligned} \sup_{t>0} |y_t - \hat{y}_t| &= \sup_{t>0} \left| \frac{1}{2\pi} \int_{-\pi}^{\pi} [Y(e^{i\omega}) - \hat{Y}(e^{i\omega})] e^{i\omega t} d\omega \right| \\ &\leq \frac{1}{2\pi} \int_{-\pi}^{\pi} |Y(e^{i\omega}) - \hat{Y}(e^{i\omega})| d\omega \\ &= \frac{1}{2\pi} \int_{-\pi}^{\pi} |H(e^{i\omega}) - \hat{H}(e^{i\omega})| |U(e^{i\omega})| d\omega \\ &\leq \left[ \frac{1}{2\pi} \int_{-\pi}^{\pi} |H(e^{i\omega}) - \hat{H}(e^{i\omega})|^2 d\omega \right]^{1/2} \left[ \frac{1}{2\pi} \int_{-\pi}^{\pi} |U(e^{i\omega})|^2 d\omega \right]^{1/2} \quad \text{Hölder Inequality} \\ &\leq \left[ \frac{1}{2\pi} \int_{-\pi}^{\pi} |H(e^{i\omega}) - \hat{H}(e^{i\omega})|^2 d\omega \right]^{1/2} \|u\|_2 \quad \text{Parseval Theorem} \\ &\leq \zeta \|H - \hat{H}\|_{\mathcal{H}_2} \end{aligned}$$ □ If $u$ is the unit impulse function $u_t = \delta_t$ then $\zeta = 1$ . The results also holds for finite sequences of length $L$ using the discrete Fourier transform. **Lemma A.2** (Impulse response error on finite sequences). *Consider filters $h, \hat{h}$ with finite length $L$ . Then, the following holds.* $$\|h - \hat{h}\|_{\infty} \leq \|H - \hat{H}\|_2$$ where $H$ and $\hat{H}$ denote the discrete Fourier transforms of $h$ and $\hat{h}$ , respectively. *Proof.* $$\begin{aligned} \|y - \hat{y}\|_{\infty} &:= \sup_{t>0} |y_t - \hat{y}_t| = \sup_{t>0} \left| \frac{1}{2\pi} \sum_{n=0}^{L-1} [Y_n - \hat{Y}_n] e^{i2\pi nt/L} \right| \\ &\leq \frac{1}{2\pi} \sum_{n=0}^{L-1} |Y_n - \hat{Y}_n| \\ &= \frac{1}{2\pi} \sum_{n=0}^{L-1} |H_n - \hat{H}_n| |U_n| \\ &\leq \left[ \frac{1}{2\pi} \sum_{n=0}^{L-1} (H_n - \hat{H}_n)^2 \right]^{1/2} \left[ \frac{1}{2\pi} \sum_{n=0}^{L-1} U_n^2 \right]^{1/2} \quad \text{Hölder Inequality} \\ &\leq \left[ \frac{1}{2\pi} \sum_{n=0}^{L-1} (H_n - \hat{H}_n)^2 \right]^{1/2} \|u\|_2 \quad \text{Parseval Theorem} \\ &= \|H - \hat{H}\|_2 \quad \text{using } \|u\|_2 = 1 \end{aligned}$$ □ ### A.3 Transfer Function of State-Space Models The transfer function (3.1) is derived by taking the $z$ -transform of input and state, $U(z) = \mathcal{Z}[u](z)$ , $X(z) = \mathcal{Z}[x](z)$ . Plugging $U(z)$ , $X(z)$ in the state equation (2.2), it holds $$zX(z) = AX(z) + BU(z) \Leftrightarrow X(z) = (zI - A)^{-1}BU(z)$$ Substituting in the output equation yields $$Y(z) = C(zI - A)^{-1}BU(z) + h_0U(z)$$ The transfer function is then defined as $$H(z) = \frac{Y(z)}{U(z)} = C(zI - A)^{-1}B + h_0. \quad (\text{A.1})$$**Alternative derivation** The transfer function can also be derived by direct $z$ -transform of the impulse response $h_t$ of the system. This derivation is useful to highlight the region of convergence of the transfer function. $$\begin{aligned} H(z) &= h_0 + \sum_{t=1}^{\infty} z^{-t} \mathbf{C} \mathbf{A}^{t-1} \mathbf{B} & h_0 \text{ is pulled out via } h_0 z^0 = h_0 \\ &= h_0 + \mathbf{C} \left[ \sum_{t=1}^{\infty} z^{-t} \mathbf{A}^{t-1} \right] \mathbf{B} & \text{multiplication distributes over sum.} \\ &= h_0 + z^{-1} \mathbf{C} \left[ \sum_{t=1}^{\infty} z^{-(t-1)} \mathbf{A}^{t-1} \right] \mathbf{B} & \text{multiply by } z/z \\ &= h_0 + z^{-1} \mathbf{C} \left[ \sum_{t=0}^{\infty} (z^{-1} \mathbf{A})^t \right] \mathbf{B} & \text{change of index and collect like terms} \end{aligned} \tag{A.2}$$ We look at the convergence of the series $\sum_{t=0}^{\infty} \|z^{-1} \mathbf{A}\|_2^t$ . We have $$\begin{aligned} \|z^{-1} \mathbf{A}\|_2 &\leq \|z^{-1}\|_2 \|\mathbf{A}\|_2 \\ &= \|r^{-1} e^{-i\omega}\|_2 \|\mathbf{A}\|_2 & \text{using } z := r e^{i\omega} \in \mathbb{C}, r, \omega \in \mathbb{R} \\ &\leq r^{-1} \|\mathbf{A}\|_2 = r^{-1} \rho(\mathbf{A}) \end{aligned}$$ The series converges to $1/(1 - r^{-1} \rho(\mathbf{A}))$ if and only if $r^{-1} \rho(\mathbf{A}) < 1$ i.e. for $r > \rho(\mathbf{A})$ . Thus, in the exterior of the disk with radius $\rho(\mathbf{A})$ , $\mathbb{D}_{\rho(\mathbf{A})} := \{z \in \mathbb{C} : |z| > \rho(\mathbf{A})\}$ , $\sum_{t=0}^{\infty} (z^{-1} \mathbf{A})^t$ converges to $(\mathbf{I} - z^{-1} \mathbf{A})^{-1}$ and $$z \in \mathbb{D}_{\rho(\mathbf{A})} \Rightarrow H(z) = h_0 + z^{-1} \mathbf{C} (\mathbf{I} - z^{-1} \mathbf{A})^{-1} \mathbf{B} = h_0 + \mathbf{C} (z \mathbf{I} - \mathbf{A})^{-1} \mathbf{B}$$ The transfer function $H(z) = h_0 + \mathbf{C} (z \mathbf{I} - \mathbf{A})^{-1} \mathbf{B}$ of a stable lumped discrete-time system is defined outside the disc in the complex plane that encloses all the eigenvalues of $\mathbf{A}$ . **Invariance of the transfer function** $H(z)$ as defined in (A.1) is a *proper*13 rational function of $z$ . In case $h_0 = 0$ , $H(z)$ is strictly proper and the denominator is monic: $$H(z) = \frac{b_1 z^{-1} + \dots + b_d z^{-d}}{1 + a_1 z^{-1} + \dots + a_d z^{-d}} \tag{A.3}$$ Specifically, the denominator could be derived from $\mathbf{A}$ with $\det(z \mathbf{I} - \mathbf{A})$ , and the numerator is $\det(z \mathbf{I} - \mathbf{A} + \mathbf{B} \mathbf{C}) + \det(z \mathbf{I} - \mathbf{A})$ . We provide a detailed derivation below in Section A.6. While state-space representation involves the analysis and synthesis of model matrices $\mathbf{A}, \mathbf{B}, \mathbf{C}$ , the transfer function is entirely characterized by the coefficients $a = (a_n)_{n=1}^d$ , $b = (b_n)_{n=1}^d$ of numerator and denominator polynomials. Notably, the transfer function is an *invariant* of the system: if we apply a change of variables to the state, the transfer function remains unchanged. **Lemma A.3.** *Coefficients $a, b$ are **invariant** under any invertible change of variables.* *Proof.* The proof can be found in [5, pp.95] and follows from the definition of *equivalence transformation*. Consider the state-space matrices of under change of variables $\hat{x} = \mathbf{K}x$ , $$\hat{\mathbf{A}} = \mathbf{K} \mathbf{A} \mathbf{K}^{-1}, \quad \hat{\mathbf{B}} = \mathbf{K} \mathbf{B}, \quad \hat{\mathbf{C}} = \mathbf{C} \mathbf{K}^{-1}, \quad \hat{h}_0 = h_0.$$ The resulting transfer function $H(z)$ can then be computed as $$\hat{H}(z) = \hat{\mathbf{C}} (z \mathbf{I} - \hat{\mathbf{A}})^{-1} \hat{\mathbf{B}} + \hat{h}_0 = \mathbf{C} \mathbf{K}^{-1} [\mathbf{K} (z \mathbf{I} - \mathbf{A}) \mathbf{K}^{-1}]^{-1} \mathbf{K} \mathbf{B} + h_0 = H(z)$$ □ ## A.4 Truncated Transfer Functions In the case of generic *truncated* (finite) impulse response filters, such that $h_t = 0$ for all $t$ greater than a certain value $L$ (which we refer to as the *length* of the filter), the transfer function is simply a polynomial in the complex variable $z$ of order $L$ , i.e. $$H(z) = \sum_{t=0}^{\infty} h_t z^{-t} = \sum_{t=0}^L h_t z^{-t} = h_0 + h_1 z^{-1} + \dots + h_L z^{-L} \tag{A.4}$$ 13i.e. such that the denominator's order is not less than the numerator's one.In case the filter is generated by a finite dimensional (lumped parameters) system, i.e. $h_t = \mathbf{C}\mathbf{A}^{t-1}\mathbf{B}$ $t = 1, \dots, L$ , then (A.4) can still be represented exactly by a rational function of order $d$ . **Lemma A.4** (Truncated rational transfer functions). *Consider the $L$ -truncated impulse response $h_t \in \ell_2(\mathbb{N})$ of a lumped-parameter filter $(\mathbf{A}, \mathbf{B}, \mathbf{C}, h_0)$ ,* $$h_t = \begin{cases} h_0 & t = 0 \\ \mathbf{C}\mathbf{A}^{t-1}\mathbf{B} & 1 \leq t \leq L \\ 0 & t > L \end{cases}.$$ Then its truncated transfer function is $$H_L(z) = \mathcal{Z}\{h\}(z) = h_0 + \mathbf{C}(\mathbf{I} - z^{-L}\mathbf{A}^L)(z\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$$ *Proof.* By definition of $z$ -transform we have $$\begin{aligned} H_T(z) &= \sum_{t=0}^{\infty} h_t z^{-t} = h_0 + \sum_{t=1}^L z^{-t} \mathbf{C}\mathbf{A}^{t-1}\mathbf{B} \\ &= h_0 + \mathbf{C} \left[ \sum_{t=1}^L z^{-t} \mathbf{A}^{t-1} \right] \mathbf{B} = h_0 + z^{-1} \mathbf{C} \left[ \sum_{t=0}^{L-1} (z^{-1}\mathbf{A})^t \right] \mathbf{B} \end{aligned} \quad (\text{A.5})$$ The sum $\sum_{t=0}^{L-1} (z^{-1}\mathbf{A})^t$ is a partial Neumann series and can be manipulated as follows. $$\begin{aligned} \sum_{t=0}^{L-1} (z^{-1}\mathbf{A})^t (\mathbf{I} - z^{-1}\mathbf{A}) &= \sum_{t=0}^{L-1} (z^{-1}\mathbf{A})^t - \sum_{t=0}^{L-1} (z^{-1}\mathbf{A})^{t+1} \\ &= \mathbf{I} - (z^{-1}\mathbf{A})^L. \end{aligned}$$ Thus, $$\sum_{t=0}^{L-1} (z^{-1}\mathbf{A})^t = (\mathbf{I} - z^{-L}\mathbf{A}^L)(\mathbf{I} - z^{-1}\mathbf{A})^{-1},$$ which plugged in (A.5) gives $H_L(z) = h_0 + \mathbf{C}(\mathbf{I} - z^{-L}\mathbf{A}^L)(z\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$ , proving the result. $\square$ Because of truncation, evaluating the transfer function $H_L(z)$ on the $L$ roots of unity $z = e^{i\omega_k}$ , $\omega_k = 2\pi k/T$ for $k = 0, \dots, L$ gives the length- $L$ discrete Fourier transform (DFT) of the filter: $$\bar{H}_k := H_L(e^{i\omega_k}) = \sum_{t=0}^{L-1} h_t e^{-i2\pi k t/L}, \quad k = 0, \dots, L-1.$$ In practice, this means that $\bar{H} \in \mathbb{C}^L$ is the FFT of $h$ , $\bar{H} = \text{FFT}_L[h]$ . If we can find an efficient and stable algorithm to evaluate $\bar{H}$ from the system matrices $(\mathbf{A}, \mathbf{B}, \mathbf{C}, h_0)$ , then the FFT-based convolution of truncated filter with an input sequence $u \in \mathbb{R}^L$ can be evaluated in $\tilde{O}(L)$ time. **Reparametrization** Assume training a LCSM equipped with SSM filters with input/target sequences to be all of length $L$ (smaller sequences can be padded with zeros to the maximum length). Thus, for training purposes, we are only interested in evaluating $\bar{H}$ for the FFT-based convolution. The truncated transfer function $H_L$ is equal to the original one with a correction term $\mathbf{I} - z^{-L}\mathbf{A}^L$ on the numerator polynomial. As already noted in S4 [6], $z^{-L}$ is conveniently equal to one on the roots of unity, $z^{i\omega_k L} = e^{-i2\pi k} = 1$ for all $k = 0, \dots, L-1$ . Hence, the correction term due to truncation becomes constant: $H_k = \mathbf{C}(\mathbf{I} - \mathbf{A}^L)(\exp(-i2\pi k/L)\mathbf{I} - \mathbf{A})^{-1}\mathbf{B}$ ; in DFT domain the truncated filter behaves as the infinitely long one with a perturbed $\mathbf{C}$ matrix $$\bar{\mathbf{C}} = \mathbf{C} - \mathbf{C}\mathbf{A}^L$$ If –as assumed– the SSM is stable $\rho(\mathbf{A}) < 1$ , (i) the transfer function is defined on the unit circle, term $\mathbf{C}\mathbf{A}^L$ will go to zero exponentially fast as $L \rightarrow \infty$ and $\bar{\mathbf{C}} = \mathbf{C}$ (as expected). As advised in [6], it is desirable to parametrize directly $\bar{\mathbf{C}}$ ; the expensive computation of the correction term $\mathbf{C}(\mathbf{I} - \mathbf{A}^L)$ is never carried out during training. Instead, the real $\mathbf{C}$ matrix can be retrieved for recurrent inference by inverting the correction term $\mathbf{C} = \bar{\mathbf{C}}(\mathbf{I} - \mathbf{A}^L)^{-1}$ , always invertible for stable systems although possibly ill conditioned by eigenvalues too close to the stability margin (the unit circle).## A.5 From Transfer Function to State-Space Suppose the coefficients of the numerator and denominator polynomials of a proper transfer function $H$ is given: $$H(z) = \frac{b_0 + b_1 z^{-1} + \dots + b_d z^{-d}}{1 + a_1 z^{-1} + \dots + a_d z^{-d}}. \quad (\text{A.6})$$ A state-space representation of the form (2.2) can be rapidly realized in two steps: 1. 1. **Get delay-free path** From (A.6) we first notice that the *bias* term $h_0$ is $h_0 = b_0$ . We thus want to isolate $b_0$ from the rest of the numerator. This can be obtained via long division (see §A.5.1) and results in $$H(z) = \frac{\beta_1 z^{-1} + \dots + \beta_N z^{-d}}{1 + a_1 z^{-1} + \dots + a_d z^{-d}} + b_0, \quad \beta_n = b_n - b_0 a_n \quad (\text{A.7})$$ 1. 2. **Get state-space matrices** Given the transfer function $H(z)$ with the isolated pass-through coefficient $b_0$ as in (A.7), we can construct the state-space matrices by *companion* canonical realization: $$\left[ \begin{array}{c|c} \mathbf{A} & \mathbf{B} \\ \hline \mathbf{C} & h_0 \end{array} \right] = \left[ \begin{array}{cccccc|c} -a_1 & -a_2 & \dots & -a_{d-1} & -a_d & 1 & 0 \\ 1 & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 & 0 & 0 \\ \hline \beta_1 & \beta_2 & \dots & \beta_{d-1} & \beta_d & 0 & b_0 \end{array} \right] \quad (\text{A.8})$$ Details on the complete *a la* [5] derivation can be found in §A.5.2. A linear system with finite-dimensional state can be equivalently characterized: by its state-space matrices $(\mathbf{A}, \mathbf{B}, \mathbf{C}, h_0)$ , by its impulse response function $h$ , or by the coefficients $a, b$ (or $\beta$ ) of the transfer function. A fourth representation is its *linear-constant-coefficients difference equation* form $$y_t = \sum_{j=0}^d b_j u_{t-j} - \sum_{n=1}^d a_n y_{t-n},$$ typically used in signal processing literature in the theory of *infinite impulse response* filters (see [46]) and known, in the context of system identification of error-in-variables models, as *auto-regressive moving-average* filters [47, 48]. ### A.5.1 Isolating the $h_0$ -term from Transfer Function by Long division If the rational transfer function $H(z)$ accounts for the $h_0$ term, then it is simply proper (order of numerator equals the order of denominator), $h_0$ is necessarily $h_0 = b_0$ (the *delay-free* path). Given the transfer function in this form, we can isolate the $b_0$ term and the strictly rational term of (A.3) by long division. We start by expanding the fraction as $$H(z) = \frac{q(z)}{p(z)} = \frac{b_0}{p(z)} + \frac{b_1 z^{-1} + \dots + b_d z^{-d}}{p(z)}.$$ and $$\frac{b_0}{p(z)} = \frac{b_0 z^d}{z^d + a_1 z^{d-1} + \dots + a_d}$$ We then use the long division method to compute $b_0/p(z)$ : $$\begin{array}{r} b_0 \\ z^d + a_1 z^{d-1} + \dots + a_d \overline{) b_0 z^d} \\ \underline{b_0 z^d + b_0 a_1 z^{d-1} + \dots + b_0 a_d} \\ - b_0 a_1 z^{d-1} - \dots - b_0 a_d \quad (\text{remainder}) \end{array}$$ to finally get $$\begin{aligned} H(z) &= b_0 - \frac{b_0 a_1 z^{d-1} + \dots + b_0 a_d}{z^d + a_1 z^{d-1} + \dots + a_d} + \frac{b_1 z^{-1} + \dots + b_d z^{-d}}{p(z)} \\ &= b_0 + \frac{(b_1 - b_0 a_1) z^{-1} + \dots + (b_d - b_0 a_d) z^{-d}}{1 + a_1 z^{-1} + \dots + a_d z^{-d}} \end{aligned}$$Note that the coefficients $b_n$ in (A.3) correspond to $b_n - b_0 a_n$ in (A.6), $b_n \leftarrow b_n - b_0 a_n$ . It is indifferent to parameterize the coefficients of the transfer function in either forms. However, if we choose the simply proper representation (A.6), we need to apply the derived correction factor to the numerator coefficients when we separate the $h_0$ term and strictly proper part of $H(z)$ . ### A.5.2 Construction of the State-Space from the Transfer Function **Chen's derivation** The derivation is based on the steps reported for the continuous-time *multi-input multi-output* case in [5]. First, we define a pseudo-state $v$ such that $$p(z)V(z) = U(z) \quad \Leftrightarrow \quad V(z) = \frac{1}{p(z)}U(z). \quad (\text{A.9})$$ Then, we define the state $x_t := (x_t^1, \dots, x_t^d) \in \mathbb{R}^d$ as $$x_t = (v_{t-1}, v_{t-2}, \dots, v_{t-d}) \quad \Leftrightarrow \quad \mathcal{Z}\{x\}(z) = X(z) = \begin{bmatrix} z^{-1} \\ \vdots \\ z^{-d} \end{bmatrix} V(z). \quad (\text{A.10})$$ From (A.9) we have $$\begin{aligned} V(z) + a_1 z^{-1} V(z) + \dots + a_d z^{-d} V(z) &= U(z) \quad \Leftrightarrow \\ V(z) = -a_1 z^{-1} V(z) - \dots - a_d z^{-d} V(z) + U(z) &\quad \Leftrightarrow \\ v_t = -a_1 v_{t-1} - \dots - a_d v_{t-d} + u_t &\quad \Leftrightarrow \text{time-delay prop. of } \mathcal{Z}\text{-transform} \\ x_{t+1}^1 = -a_1 x_t^1 - \dots - a_d x_t^d + u_t &\quad \Leftrightarrow \text{by def. of state (A.10)}. \end{aligned}$$ Thus, we have the overall recurrence $$\begin{aligned} x_{t+1}^1 &= -a_1 x_t^1 - \dots - a_d x_t^d + u_t \\ x_{t+1}^2 &= x_t^1 \\ &\vdots \\ x_{t+1}^d &= x_t^{d-1} \end{aligned}$$ which can be written in matrix form as $$x_{t+1} = \begin{bmatrix} -a_1 & -a_2 & \dots & -a_N \\ 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 \end{bmatrix} x_t + \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix} u_t$$ The output spectrum is then given by $$\begin{aligned} Y(z) &= H(z)U(z) = \frac{q(z)}{p(z)}U(z) + b_0 U(z) \\ &= q(z)V(z) + b_0 U(z) \quad \text{by def. of } V(z). \end{aligned}$$ Therefore, $$\begin{aligned} Y(z) &= q(z)V(z) + b_0 U(z) = [\beta_1 \quad \beta_2 \quad \dots \quad \beta_N] \begin{bmatrix} z^{-1} \\ z^{-2} \\ \vdots \\ z^{-d} \end{bmatrix} V(z) + b_0 U(z) \\ &= [\beta_1 \quad \beta_2 \quad \dots \quad \beta_d] X(z) + b_0 U(z) \end{aligned}$$ and the output equation in time-domain is given by $$y_t = [\beta_1 \quad \beta_2 \quad \dots \quad \beta_d] x_t + b_0 u_t.$$ yielding state-space matrices (A.8).## A.6 From State-Space to Transfer Function We detail an implementation oriented method to compute the coefficients $(a_n)_{n=1}^d, (b_n)_{n=0}^d$ of a SSM's transfer function. Recall that $$H(z) = C[zI - A]^{-1}B + h_0 = \frac{C \text{Adj}(zI - A)B + \det(zI - A)h_0}{\det(zI - A)} \quad (\text{A.11})$$ Hence, the denominator coefficients $(a_n)_{n=1}^d$ are simply the coefficients of the characteristic polynomial of matrix $A$ . They can be easily obtained by 1. computing the eigenvalues of $A$ and 2. calculating the coefficients of the polynomial whose roots are such eigenvalues. On the other hand, the numerator apparently involves more complex symbolic manipulation. This can be simplified recalling a classic matrix-determinant identity: **Lemma A.5** ([49]). *Let $M$ , $B$ , and $C$ respectively denote matrices of orders $d \times d$ , $d \times 1$ , and $1 \times d$ . Then,* $$\det(M + BC) = \det(M) + C \text{Adj}(M)B.$$ Applying Lemma A.5 to (A.11) we obtain $$H(z) = \frac{\det(zI - A + BC) + \det(zI - A)(h_0 - 1)}{\det(zI - A)}.$$ Let $\text{poly}(r)$ denote the coefficients of the polynomials with roots $r = (r_1, \dots, r_d)$ . Then $a = \text{poly}(\text{eig}(A))$ . Since $A$ and $A - BC$ are of equal dimension, their characteristic polynomials have equal order and therefore $$b = \text{poly}(\text{eig}(A - BC)) + \text{poly}(\text{eig}(A))(h_0 - 1)$$ Listing 1: State-space $\rightarrow$ transfer function conversion code ``` def get_tf_from_ss(A,B,C,h0): a = poly(eig(A)) b = poly(eig(A - outer(B,C))) + (h0-1)*a return a, b ``` ## A.7 State-Space Representation of Truncated Filters. A truncated filter $h_0, \dots, h_L$ – as the ones found in any standard convolutional neural network – can be represented by a $L$ -dimensional companion canonical SSM. The filter's transfer function $H(z) = h_0 + h_1 z^{-1} + \dots + h_L z^{-L}$ is polynomial, i.e. a rational function with the denominator's coefficients set to zero. Following the canonical realization process detailed in Section A.5, the truncated filter has state-space form: $$\begin{aligned} x_{t+1} &= \begin{bmatrix} 0 & 0 & \dots & 0 & 0 \\ 1 & 0 & \dots & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix} x_t + \begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix} u_t \\ y_t &= [h_1 \quad h_2 \quad \dots \quad h_L] x_t + h_0 u_t. \end{aligned}$$ If $x_0 = \mathbb{0}_L$ and $u_t = 0$ for negative $t$ , then at each $t > 0$ the state is a shifted copy of the input sequence $x_t = (u_{t-1}, \dots, u_{t-L}) \in \mathbb{R}^L$ . Nonetheless, the asymptotic complexity of computing one recurrent step is $\mathcal{O}(L)$ as it requires only a shift operation and a length- $L$ dot product $$\begin{aligned} x_{t+1}^1 &= u_t \\ x_{t+1}^{2:L} &= \text{shift}(x_t) \\ y_t &= \langle h_{1:L}, x_t \rangle + h_0 u_t. \end{aligned} \quad (\text{A.12})$$ The memory footprint is also $\mathcal{O}(L)$ . In [1] it is proposed the use of shift-type SSMs to parametrize one of the filters of the H3 block.## A.8 Efficient Computation of State-Space Models ### A.8.1 Fast Evaluation of the Transfer Function Computing $H(z)$ at any point $z \in \mathbb{C}$ concerns the evaluation of the $d$ -order polynomial of numerator and denominator, $$H(z) = \frac{q(z)}{p(z)} = \frac{\sum_{n=1}^d b_n z^{-n}}{1 + \sum_{n=1}^d a_n z^{-n}}$$ In practice, we are mainly interested in a fast algorithm that allows computing $H$ on the $L$ roots of unity to obtain the DFT of the filter. The DFT of the filter can be then readily used to perform a FFT-based convolution with a length- $L$ input sequence $u$ or to recover the impulse response function via inverse DFT. We prove the following: **Lemma A.6.** *Given the coefficients $a, b$ of the transfer function, the frequency and impulse response of the filter can be evaluated in $\tilde{\mathcal{O}}(L)$ time.* *Proof.* The result is proven showing that the transfer function can be evaluated in $\tilde{\mathcal{O}}(L)$ time on the $L$ roots of unity. The fastest method to evaluate polynomials on $L$ arbitrary points $z$ of the complex plane is generally the Horner's scheme. This method is based on a sequence of nested multiplications and computes the polynomial from its vector of coefficients, delivering a time complexity of $\mathcal{O}(dL)$ . More explicitly, Horner's scheme determines $p(z)$ as $p(z) = ((\cdots((a_d z^{-1} + a_{d-1})z^{-1} + a_{d-2})\cdots)z^{-1} + a_2)z^{-1} + a_1)z^{-1} + 1$ . Each step involves a multiplication and an addition, making a total of $2d$ operations per evaluation point. Thus, for $L$ points, the total number of operations amounts to $\mathcal{O}(dL)$ . Effectively, Horner's approach implements the matrix-vector product of an $L$ -by- $(d+1)$ Vandermonde matrix $\mathbf{V} \in \mathbb{C}^{L \times (d+1)}$ constructed by $L$ evaluation points $(z_0, \dots, z_{L-1})$ with the vector of coefficients $a = (1, a_1, \dots, a_d)^\top$ : $$\begin{bmatrix} p(z_0) \\ p(z_1) \\ \vdots \\ p(z_{L-1}) \end{bmatrix} = \begin{bmatrix} 1 & z_0^{-1} & z_0^{-2} & \cdots & z_0^{-d} \\ 1 & z_1^{-1} & z_1^{-2} & \cdots & z_1^{-d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & z_{L-1}^{-1} & z_{L-1}^{-2} & \cdots & z_{L-1}^{-d} \end{bmatrix} \begin{bmatrix} 1 \\ a_1 \\ \vdots \\ a_d \end{bmatrix} = \mathbf{V}a$$ Significantly, if the polynomial is required to be evaluated at the roots of unity, the Vandermonde matrix simplifies corresponds to the $L \times (d+1)$ DFT matrix. Further, zero-padding the coefficient vector to length $L$ , enables the use a single length- $L$ FFT to compute the matrix-vector product in $\tilde{\mathcal{O}}(L)$ time. Thus, the numerator and denominator polynomials of the transfer function can be evaluated, on the roots of unity, in $\tilde{\mathcal{O}}(L)$ time by taking the FFT of the padded numerator / denominator coefficients $a, b$ and subsequently dividing element-wise the two sequences as $\text{FFT}_L[b]/\text{FFT}_L[a]$ . The overall time complexity to obtain the impulse response is also $\tilde{\mathcal{O}}(L)$ since $h$ can be recovered taking an inverse FFT of the frequency response. $\square$ ### A.8.2 Fast Companion Recurrence The recurrent step of a generic SSM (2.2) with dense system matrices usually requires $\mathcal{O}(d^2)$ operations due to the matrix-vector product $\mathbf{A}x_t$ . We show how the recurrence of SSMs in *companion canonical form*, i.e. with system's matrices (A.8), requires only $\mathcal{O}(d)$ operations. **Lemma A.7.** *The recurrent step of a state-space model in companion canonical form (A.8) can be evaluated in $\mathcal{O}(d)$ time and memory.* *Proof.* The companion state matrix $\mathbf{A}$ can be broken down into a lower shift matrix $\mathbf{L}_N$ and a low-rank term. Particularly, with $e_1$ the first element of the canonical basis of $\mathbb{R}^N$ and $\alpha = (a_1, \dots, a_N)$ , we have $$\mathbf{A} = \mathbf{L}_N - e_1 \otimes \alpha.$$ It follows that the recurrent update can be simplified to $$\begin{aligned} x_{t+1} &= (\mathbf{L}_N - e_1 \otimes \alpha) x_t + \mathbf{B}u_t \\ y_t &= \mathbf{C}x_t + b_0 u_t \end{aligned}$$ The peculiarity of this formulation is that we never need to construct the matrices to perform the recurrence. In particular we have: $$\begin{aligned} x_{t+1}^1 &= u_t - \alpha^\top x_t \\ x_{t+1}^{2:N} &= \text{shift}(x_t) \\ y_t &= \beta^\top x_t + b_0 u_t \end{aligned}$$Thus, each step only requires two inner products ( $d$ multiplications and $d$ sums each) and one shift operation, totaling $\mathcal{O}(d)$ operations. $\square$ The proof of Lemma A.7 yields the practical implementation of the recurrence: Listing 2: Python implementation of the companion canonical recurrence ``` def step(x, u, alpha, beta, b0): y = dot(beta, x) + b0 * u lr = u - dot(alpha, x) x = roll(x) x[0] = lr return x, y ``` ### A.8.3 Canonization of State-Space Models The companion canonical form discussed in Section A.5 is the ideal representation to deploy SSM-based convolutional layers: *i*) it comes with a $\mathcal{O}(d)$ fast recurrence and *ii*) allows to swiftly switch between time and frequency domains with a direct mapping between state-space matrices and coefficient of the transfer function (which in turn allow $\tilde{\mathcal{O}}(L)$ fast convolutions). Aside from [50], which directly parametrizes S4 layers in companion canonical form, all the other parameterizations [12, 6, 32, 17, 33] can be *converted* (*canonized*), under mild assumptions. **Lemma A.8** (Canonization of SSMs). *Any state-space model (2.2) with proper transfer function can be converted in companion canonical form.* *Proof.* The result can be proved following the two-step conversion process. 1. 1. **Get the coefficients of the transfer function:** Given the original state-space matrices $(A, B, C, h_0)$ , the transfer function is given by $H(z) = C(zI - A)^{-1}B + h_0$ . A proper rational function has the form $H(z) = q(z)/p(z)$ where the numerator $q(z)$ has coefficients $b = (b_n)_{n=0}^d$ and the denominator has coefficients $a = (a_n)_{n=0}^d$ ( $a_0 = 1$ since $p$ is monic). As shown in Section A.6, the coefficients of the transfer function can be extracted in closed-form as $b = \text{poly}(\text{eig}(A - BC)) + \text{poly}(\text{eig}(A))(1 - h_0)$ and $a = \text{poly}(\text{eig}(A))^{14}$ ; 2. 2. **Construct companion matrices** Given the coefficients $a$ and $b$ a new set of canonical state-space matrices which realize the transfer function can be obtained following the recipe of Section A.5. The resulting companion SSM is equivalent the the original one since they share the same transfer function. $\square$ 14 $\text{eig}(A)$ contains the eigenvalues of $A$ . $\text{poly}(r)$ yields the coefficients of the polynomial whose roots are the elements of $r \in \mathbb{C}^d$ .## B LaughingHyena: Further Details ### B.1 Parametrization of Modal Interpolators **Complex-conjugate states** Assuming even distilling dimension $d$ , we pick poles $\lambda_n$ and residues $\mathbb{R}_n$ in complex-conjugate pairs: $$\begin{aligned} A &= \text{diag}(\lambda_1, \dots, \lambda_{d/2}, \lambda_1^*, \dots, \lambda_{d/2}^*) \\ C &= \frac{1}{2}[R_1, \dots, R_{d/2}, R_1^*, \dots, R_{d/2}^*] \end{aligned} \quad (\text{B.1})$$ which allow partitioning the state-space matrices as $$A = \begin{bmatrix} \lambda & \\ & \lambda^* \end{bmatrix}, \quad C = \frac{1}{2} \begin{bmatrix} R & R^* \end{bmatrix}, \quad (\text{B.2})$$ where $$\lambda = \text{diag}(\lambda_1, \dots, \lambda_{d/2}) \quad \text{and} \quad R = [R_1, \dots, R_{d/2}]. \quad (\text{B.3})$$ If we also partition the state as $x = (\bar{x}, \tilde{x})$ , $\bar{x}, \tilde{x} \in \mathbb{C}^{d/2}$ , the resulting recurrence has the form $$\begin{aligned} \bar{x}_{t+1} &= \lambda \bar{x}_t + \mathbb{1}_{d/2} u_t \\ \tilde{x}_{t+1} &= \lambda^* \tilde{x}_t + \mathbb{1}_{d/2} u_t \end{aligned} \quad (\text{B.4})$$ We have $$\bar{x}_t = \lambda^t \bar{x}_0 + \sum_{j=0}^{t-1} \lambda^{t-j-1} \mathbb{1}_{d/2} u_t, \quad \tilde{x}_t = [\lambda^*]^t \tilde{x}_0 + \sum_{j=0}^{t-1} [\lambda^*]^{t-j-1} \mathbb{1}_{d/2} u_t \quad (\text{B.5})$$ Thus, if $\tilde{x}_0 = \bar{x}_0^*$ , then $\tilde{x}_t = \bar{x}_t^*$ for all $t > 0$ . Hence, at inference time we only need to propagate forward half of the state – say $\bar{x}$ – and then compute the output as $$\begin{aligned} y_t &= D u_t + \frac{1}{2}(R \bar{x}_t + R^* \bar{x}_t^*) \\ &= D u_t + \Re\{R \bar{x}_t\} \\ &= D u_t + \Re\{R\} \Re\{\bar{x}_t\} - \Im\{R\} \Im\{\bar{x}_t\} \end{aligned} \quad (\text{B.6})$$ This parametrization allows to update only half of the state, reducing the time and memory cost compared to an *unconstrained* linear system with complex coefficients. However, the *implicitly* achieved realness of the output (assuming $D = h_0 \in \mathbb{R}$ and $u_t \in \mathbb{R}$ ) comes at a cost of expressivity: such a system is equivalent to an unconstrained complex linear system of dimension $d/2$ of which we only keep the real part of the output. **Poles and residues** For the modal interpolation, the parametrization is analogous to the one of a diagonal state space model [33]. Poles $\lambda_n$ and residues $R_n$ need both to be complex numbers. In [33] the authors suggest parametrizing real and imaginary components of $B$ and $C$ matrices while representing the eigenvalues $\lambda_n$ in polar form, $\lambda_n = r_n e^{i\alpha_n}$ with $r_n$ and $\alpha_n$ being themselves exponential functions of the actual trainable parameters, $r_n = e^{-e^{\nu_n}}$ , $\alpha_n = e^{\zeta_n}$ leading to $$\lambda_n = \exp\{-\exp\{\nu_n\} + i \exp\{\zeta_n\}\}, \quad \nu_n, \zeta_n \in \mathbb{R} \quad (\text{B.7})$$ This ensures stability of the poles $|\lambda_n| < 1$ and positive-only phases $\alpha_n$ . For the purpose of distillation we propose a simplified parametrization as follows: 1. 1. We only parametrize the $C$ vector. Parametrizing both $B$ and $C$ is redundant and increases the computational cost of performing each step of the recurrence. The residues $R_n$ correspond in fact to $R_n = C_n B_n$ of a diagonal state space model. Setting $B = \mathbb{1}_d$ saves parameters without harming expressivity. Further if $B$ is different from $\mathbb{1}_k$ it needs to be multiplied to $u_t$ at each recurrence step. $C_n = R_n = \Re[R_n] + i \Im[R_n]$ and $\Re[R_n], \Im[R_n]$ are the trainable parameters of the residue. 2. 2. For the purpose of distillation we have no benefit in forcing the eigenvalues of the model to be stable, i.e. constrained to lie strictly inside the unit circle. Instead, such constrain may actually harm the expressivity of the approximant. We choose the simpler parametrization $\lambda_n = r_n e^{i\alpha_n}$ , $r_n, \alpha_n \in \mathbb{R}$ .## B.2 Distillation as Rational Interpolation **Distillation as *rational interpolation*** Approximating a filter with an SSM can be thus achieved by fitting a proper rational function to the (truncated) transfer function of the original filter $H_L(z) := \sum_{t=0}^L h_t z^{-t}$ . That is, $$\text{Find } a, b \text{ such that } h_0 + h_1 z^{-1} + \dots + h_L z^{-L} \approx h_0 + Q_b(z)/P_a(z). \quad (\text{B.8})$$ A modern15 way to solve this problem by $\mathcal{H}_2$ error minimization via gradient descent16. We can use the Fast Fourier Transform (FFT) to evaluate both the target and distilled transfer functions and solve: $$\min_{a, b \in \mathbb{R}^d} \sum_{k=0}^L |\text{FFT}_L[h]_k - h_0 - \text{FFT}_L[b]_k / \text{FFT}_L[a]_k|^2. \quad (\text{B.9})$$ To ensure stability of the distilled filters and well-conditioned gradient descent dynamics, the roots of the denominator polynomial must strictly lie inside the unit circle ( $\rho(\mathbf{A}) < 1$ ). This, in turn, requires constraining the coefficients $a$ into the region $\{a : \text{poly}(a) \text{ is stable}\}$ which is by itself an open research problem [52, 53]. Experimentally, we observe that standard coefficient normalization techniques overly restrict the parameters space and lead to poor distillation performances at reasonable order. 15In the late 19th century, Henri Padé had already proposed a closed-form solution of the above problem that achieves $o(z^{-L})$ error for $z \rightarrow \infty$ using $L=2d$ samples of the impulse response. His method [31] solves a $L$ -dimensional linear problem that, however, is known to often become numerically ill-conditioned even with small $d$ [51] 16In the case of finite sequences, the $\mathcal{H}_2$ norm becomes the standard Euclidean metric evaluated on the $L+1$ roots of unity, i.e. $(\sum_{k=0}^L |H(e^{i2\pi k/(L+1)})|^2)^{1/2}$ .## C Proofs ### C.1 Proof of Lemma 2.1 Generating $K$ tokens with a long convolution layer (2.1) from a length- $T$ prompt has time complexity $\mathcal{O}(T \log_2 T + TK + K^2)$ and requires $\mathcal{O}(L)$ memory. *Proof.* We compute the time complexity memory of a length- $T$ prompt processing (*pre-filling*) and subsequent auto-regressive decoding of $K$ tokens. The auto-regressive generation of long convolution computes the next token as by $$t = T, \dots, T + K - 1 \Rightarrow y_t = \sum_{j=0}^{t-1} h_{t-j} y_j \quad (\text{C.1})$$ The pre-filling step is needed to prime this recurrence by computing the first $T$ outputs till $y_{T-1}$ from the length- $T$ prompt $u$ . This is just a convolution between two length- $T$ signal and requires $\mathcal{O}(T \log_2 T)$ time and linear memory. The auto-regressive decoding of $K$ tokens requires $K$ steps (C.1) with the length of the sequences increasing by 1 at each step. Thus we have a total asymptotic complexity of $$\sum_{k=0}^{K-1} (T + k) = TK + \frac{1}{2}K(K + 1). \quad (\text{C.2})$$ and requires at worst ( $k = K - 1$ ) to store the length $T + K = L$ generated output sequence, i.e. $\mathcal{O}(L)$ memory. In the limit we thus have a total time complexity of $\mathcal{O}(T \log_2 T + TK + K^2)$ and $\mathcal{O}(L)$ memory. $\square$ ### C.2 Proof of Lemma 2.2 Generating $K$ tokens with a SSM (2.2) from a length- $T$ prompt has time complexity $\mathcal{O}(T \log_2 T + dK)$ and requires $\mathcal{O}(d)$ memory. *Proof.* In autoregressive mode, the cost of generating one token is the cost of evaluating the state recurrence (2.2). Each step then requires $\mathcal{O}(d)$ time and memory for the class of SSMs considered in this work (see Lemma A.8). Hence, generating $K$ tokens costs $\mathcal{O}(dK)$ time and constant $\mathcal{O}(d)$ memory (we only need to store the current state). The recurrence is initialized for autoregressive generation with the post-prompt state $x_{T-1}$ and output $y_{T-1}$ . The latter can be recovered in linear time and memory $\mathcal{O}(T)$ by definition $y_{T-1} = \sum_{j=0}^{T-1} h_{T-1-j} u_j$ (assuming to have the impulse response $h$ available) and state $x_{T-1}$ in $\mathcal{O}(dT)$ time and $d$ memory through the recurrence. The overall asymptotic cost is therefore $\mathcal{O}(dL)$ time and $\mathcal{O}(d)$ memory. $\square$ Note that, for prompts and SSMs of practical sizes we usually have $d > \log_2 T$ . In such a case the state $x_{T-1}$ can be computed in $T \log_2 T$ time rather than $dT$ by Proposition 3.2. ### C.3 Proof of Lemma 2.3 Generating $K$ tokens with self-attention from a length- $T$ prompt has time complexity $\mathcal{O}(T^2 + TK + K^2)$ and requires $\mathcal{O}(L)$ memory. *Proof.* The proof is identical to the one of Lemma 2.1, with the only difference of a quadratic asymptotic cost $\mathcal{O}(T^2)$ to process the prompt obtain the $kv$ cache. $\square$ Self-attention suffers with long contexts: it is significantly more expensive in prefilling than long convolutions and SSMs due to its quadratic cost. Nonetheless, in autoregressive mode, self-attention reaches the same overall asymptotic complexity $\mathcal{O}(TK + K^2)$ as long convolutions (with the memory overhead of having to cache $k$ and $v$ ).## C.4 Proof of Proposition 3.1 If $\mathbf{A}$ has semi-simple eigenvalues $\lambda_n \in \mathbb{C}$ , then the transfer function of the system can be decomposed as $\hat{H}(z) = \sum_{n=1}^d R_n / (z - \lambda_n)$ where $R_n \in \mathbb{C}$ is the residue associated with the pole $\lambda_n$ . *Proof.* If $\mathbf{A}$ is semi-simple, then it is diagonalized by a basis $\mathbf{V}$ of eigenvectors; it admits an eigenvalue decomposition $\text{diag}(\lambda) = \mathbf{V}\mathbf{A}\mathbf{V}^{-1}$ where $\lambda = (\lambda_1, \dots, \lambda_d) \in \mathbb{C}^d$ contains the eigenvalues of $\mathbf{A}$ . Projecting the state onto the basis of eigenvectors, $s := \mathbf{V}x$ , the state space model is transformed into modal form: $$\begin{aligned} s_{t+1} &= \mathbf{V}\mathbf{A}\mathbf{V}^{-1}s_t + \mathbf{V}\mathbf{B}u_t & \Leftrightarrow \quad s_{t+1} &= \text{diag}(\lambda)s_t + \tilde{\mathbf{B}}u_t \\ y_t &= \mathbf{C}\mathbf{V}^{-1}s_t & & y_t &= \tilde{\mathbf{C}}s_t \end{aligned}$$ where $\tilde{\mathbf{B}} := \mathbf{V}\mathbf{B} = (\tilde{b}_n)_{n=1}^d$ and $\tilde{\mathbf{C}} := \mathbf{C}\mathbf{V}^{-1} = (\tilde{c}_n)_{n=1}^d$ . In modal form, the state equations are decoupled, i.e. $$\begin{aligned} s_{t+1}^n &= \lambda_n s_t^n + \tilde{b}_n u_t \\ y_t &= \sum_{n=1}^d \tilde{c}_n s_t^n. \end{aligned}$$ Taking the $z$ -transform of the output equation and each state equation yields $$\begin{aligned} S_n(z) &= \mathcal{Z}[s^n](z) = \frac{\tilde{b}_n}{z - \lambda_n} U(z) \quad n = 1, \dots, d \\ Y(z) &= \sum_{n=1}^d \tilde{c}_n S_n(z) \end{aligned}$$ Thus, the overall transfer function is $$H(z) = \frac{Y(z)}{U(z)} = \sum_{n=1}^d \frac{\tilde{b}_n \tilde{c}_n}{z - \lambda_n}$$ Letting $R_n = \tilde{b}_n \tilde{c}_n$ , proves the result. $\square$ ## C.5 Proof of Lemma 3.1 The distilled filter $\hat{h}$ in modal form (3.2) can be computed in $\mathcal{O}(dL)$ time from its modal form and in $\tilde{\mathcal{O}}(L)$ from its rational form. Recalling (3.2), $\hat{h}_t = \sum_{n=1}^d R_n \lambda_n^{t-1}$ , $R_n, \lambda_n \in \mathbb{C}, t > 0$ , the $\mathcal{O}(dL)$ complexity of the impulse response is apparent: for each of the $t = 1, \dots, L$ , $\hat{h}_t$ can be computed in $\mathcal{O}(d)$ time. The $\tilde{\mathcal{O}}(L)$ cost from the rational form follows by Lemma A.6. ## C.6 Proof of Theorem 3.2 Let $h$ be a length- $L$ filter, $\hat{h}$ a distilled filter of order $d < L$ and let $\mathbf{S}_L, \hat{\mathbf{S}}_L$ be the respective Hankel matrices. Then $\inf_{\hat{\mathbf{S}}_L} \|\mathbf{S}_L - \hat{\mathbf{S}}_L\|_2 = \sigma_d$ . *Proof.* The theorem characterizes the best-case scenario in terms of approximation error of the distilled SSMs or a certain order $d$ where it is clear that $\text{rank } \hat{\mathbf{S}} \leq d$ . This theorem is a direct application of the Adamyan-Arov-Krein (AAK) theory of infinite Hankel operators [7]. Let $\hat{\mathbf{S}}_L^* = \arg \inf_{\hat{\mathbf{S}}_L} \|\mathbf{S}_L - \hat{\mathbf{S}}_L\|_2$ ; the AAK theorem says that every causal system can be optimally approximated by another causal system of lower dimension. Optimal here means $$\inf \|\mathbf{S}_L - \hat{\mathbf{S}}_L^*\| = \inf \|\mathbf{S}_L - \mathbf{K}\|$$ where the first infimum is taken over all Hankel matrices $\hat{\mathbf{S}}_L^*$ and the second over all arbitrary matrices $\mathbf{K}$ (see [19, Chapter 8] and [54] for further details and references). $\square$### C.7 Proof of Proposition 3.3 The filter (3.2) has a state space matrices $\mathbf{A} = \text{diag}(\lambda_1, \dots, \lambda_d) \in \mathbb{C}^{d \times d}$ , $\mathbf{B} = (1, \dots, 1)^\top \in \mathbb{C}^{d \times 1}$ , $\mathbf{C} = (R_1, \dots, R_d) \in \mathbb{C}^{1 \times d}$ , $\mathbf{D} = h_0$ whose step can be evaluated in $\mathcal{O}(d)$ time and memory. $\mathbf{D}$ is set to $h_0$ by default. The result is proven showing that (3.2) can be written in the form $\hat{h}_t = \mathbf{C}\mathbf{A}^{t-1}\mathbf{B}$ for $t > 0$ . If we choose $\mathbf{A} = \text{diag}(\lambda)$ then the impulse response becomes $$\hat{h}_t = \sum_{n=1}^d R_n \lambda_n^{t-1} = \mathbf{C}[\text{diag}(\lambda)]^{t-1}\mathbf{B} = \sum_{n=1}^d \mathbf{C}_n \mathbf{B}_n \lambda_n^{t-1}$$ The choice $\mathbf{B}_n = 1$ for all $n = 1, \dots, d$ , $\mathbf{B} = \mathbb{1}_d$ and $\mathbf{C}_n = R_n$ finalizes a modal canonical state-space realization of the distilled filter. The $\mathcal{O}(d)$ time complexity of the corresponding recurrent step is guaranteed by the decoupling of each state equation from another, $$\begin{aligned} x_{t+1}^n &= \lambda_n x_t^n + u_t \quad n = 1, \dots, d \\ y_t &= \sum_{n=1}^d R_n x_t^n + h_0 u_t. \end{aligned}$$ Each of the $d$ state equations can be computed (in parallel) in $\mathcal{O}(1)$ time. The output equation is a dot product requiring $d$ multiplications and $d$ additions, hence the $\mathcal{O}(d)$ time complexity of the recurrence. ### C.8 Proof of Proposition 3.2 $x_T = (v_T, \dots, v_{T-d})$ where $v = g * u$ and $g$ is the filter whose transfer function is $1/\text{den}(\hat{H})(z)$ and can be evaluated in $\tilde{\mathcal{O}}(T)$ . Without loss of generality, let us assume to have converted the distilled filter in canonical form (i.e. we have unrestricted access to the coefficients of the rational transfer function) and let $\mathbf{D} = 0$ . We use the notation of Section A.5. In $z$ -domain, the state-to-input relation is given by $$Y(z) = \mathbf{C}X(z) = [\beta_1 \quad \dots \quad \beta_d] X(z)$$ On the other hand $Y(z) = \hat{H}(z)U(z) = q(z)/p(z)U(z)$ . Therefore, $$\begin{aligned} [\beta_1 \quad \dots \quad \beta_d] X(z) &= \frac{q(z)}{p(z)}U(z) \\ \Leftrightarrow [\beta_1 \quad \dots \quad \beta_d] X(z) &= [\beta_1 \quad \dots \quad \beta_d] \begin{bmatrix} z^{-1} \\ \vdots \\ z^{-d} \end{bmatrix} \frac{1}{p(z)}U(z) \\ \Leftrightarrow X(z) &= \begin{bmatrix} z^{-1} \\ \vdots \\ z^{-d} \end{bmatrix} \frac{1}{p(z)}U(z) \end{aligned}$$ Let $V(z) = U(z)/p(z)$ . From the shift property of the $z$ -transform it holds, $$\mathcal{Z}\{x\}(z) = X(z) = \begin{bmatrix} z^{-1} \\ \vdots \\ z^{-d} \end{bmatrix} V(z) \quad \Leftrightarrow \quad x_t = (v_{t-1}, v_{t-2}, \dots, v_{t-d}) \quad \forall t > 0.$$ $v$ can be obtained in $\tilde{\mathcal{O}}(L)$ time via an FFT-convolution of the input $u$ and $g$ , the filter resulting from inverse transforming $1/p(z)$ . The proof is convoluted setting $t = L$