FlashFFTConv: Efficient Convolutions for Long Sequences with Tensor Cores

🎯 Pitch

Long convolutions theoretically beat attention on long sequences, but in practice they’re crippled by FFTs that barely touch tensor cores—FlashFFTConv rewrites the FFT as matrix multiplies and fuses entire convolutions into a single kernel, delivering up to 7.9× speedup. This lets Hyena language models and M2-BERT match the quality of models twice their size, while enabling the first DNA model to process the longest human genes at 2.3 million base pairs.

1. Executive Summary

This paper introduces FlashFFTConv, a system that optimizes the Fast Fourier Transform (FFT) convolution for long sequences on modern GPU tensor cores using a Monarch decomposition — rewriting the FFT as a series of matrix-matrix multiply operations that better utilize specialized hardware. The authors study convolutions across a range of sequence-modeling benchmarks (MATH, Path-512 from the Long Range Arena, language modeling on the PILE, and DNA modeling with HyenaDNA) using PaLM 2-S* and Hyena-family models, and develop two key mechanisms: the order-p Monarch decomposition (trading off FLOP cost against I/O cost by adjusting matrix sizes at different sequence lengths) and kernel fusion techniques (eliminating expensive trips to HBM by keeping intermediate results in SRAM for sequences up to 32K). FlashFFTConv achieves up to 7.93× speedup over PyTorch FFT convolutions and enables up to 4.4× end-to-end speedup in convolutional sequence models, while reducing memory footprint by up to 5.60× through recomputation and fusion — establishing that FFT convolutions can match or exceed the wall-clock efficiency of highly optimized Transformers using FlashAttention-v2 at sequence lengths of 2K and longer, but only when the Monarch decomposition is tuned to balance tensor core utilization against memory bandwidth constraints.

2. Context and Motivation

The Core Problem: FFT Convolutions Are Theoretically Efficient but Practically Slow

The paper addresses a stark contradiction that has held back convolutional sequence models from competing with Transformers in real-world deployment. On paper, long convolutions — where the kernel is as long as the input sequence — enjoy an asymptotic complexity of $O(N \log N)$ in sequence length $N$ when computed via the Fast Fourier Transform (FFT) convolution algorithm (Equation 1). This compares favorably to the self-attention mechanism in Transformers, which costs $O(N^2)$ in the naive case. Theoretically, convolutions should be faster, especially as sequence length grows.

In practice, however, convolutional sequence models "still lag behind Transformers in wall-clock time" (Section 1, paragraph 2). The paper identifies a specific, concrete reason: the FFT convolution algorithm has poor hardware utilization on modern accelerators, despite its asymptotic efficiency. The gap between theory and practice is not marginal — it is large enough that even convolutional models demonstrating state-of-the-art reasoning abilities on long-sequence tasks (language modeling, DNA analysis, high-resolution vision, audio generation) cannot match the throughput of highly optimized Transformer implementations.

This gap is particularly frustrating because the systems community has already demonstrated that Transformers can be pushed to the hardware limits. The paper explicitly cites FlashAttention-v2, which achieves "more than 72% FLOP utilization end-to-end" (Section 1, paragraph 3). The implication is clear: if attention — a quadratic-cost operation — can be optimized to near-peak hardware efficiency, then there is no fundamental reason why a fundamentally $O(N \log N)$ algorithm like the FFT convolution cannot achieve similar or better utilization. The problem must be in how the FFT convolution is mapped to hardware, not in any inherent limitation of the convolution itself.

Why This Problem Matters: Beyond a Single Benchmark

The paper is not merely optimizing a low-level primitive for its own sake. The stakes are high across multiple dimensions of machine learning research and deployment:

1. Enabling longer-context applications. The paper highlights three domains where longer sequences directly translate to better capabilities:

DNA modeling (Section 4.3, Section 5): The longest human genes, such as the dystrophin gene, span approximately 2.3 million base pairs. Modeling DNA at single-nucleotide resolution — which recent work has suggested improves downstream quality compared to tokenization or downsampling — requires processing sequences of this length. Prior to FlashFFTConv, no model could embed these genes at full resolution; the paper explicitly claims to achieve this for the first time (Section 4.3, Table 8).
High-resolution computer vision (Section 4.1): The Path-512 task from the Long Range Arena benchmark requires classifying whether two dots are connected in a $512 \times 512$ image flattened to a sequence of length 262,144. The paper states that "no model had previously achieved better than 50%" accuracy on this task (Abstract), and that existing implementations "fail to achieve better-than-random (50%) accuracy on Path-512 due to out of memory errors and a lack of support for such long sequences" (Section 4.1).
Raw audio generation (Section 5): Modeling speech directly from raw waveforms at native sampling rates (e.g., 64 kHz) inherently requires long sequences. The paper benchmarks SaShiMi, an audio generation model operating on 1-second audio clips at 64 kHz — a sequence length of 64,000.

These are not contrived benchmarks. They represent real scientific and engineering problems where longer sequences are a hard requirement, not merely a convenience. If the FFT convolution bottleneck prevents convolutional models from operating at these sequence lengths, it effectively blocks a whole class of architectures from being applied to high-impact problems.

2. The efficiency-quality feedback loop. The paper articulates a mechanism that goes beyond raw speedup (Section 4.1): given a fixed compute budget, a model with higher training throughput can simply see more data during pretraining. When FlashFFTConv speeds up convolution operations, the downstream models trained with it can process more tokens in the same wall-clock time. The paper quantifies this: Hyena-GPT-s achieves 2.3 points better perplexity on the PILE, and M2-BERT-base achieves 3.3 points higher average GLUE score — improvements that the authors explicitly equate to "matching models with twice the parameter count" (Abstract, Section 4.1, Table 1).

This feedback loop is significant because it means that systems optimizations are not merely "making the same model faster" — they are changing the effective model capacity under a realistic compute constraint. A team with a fixed GPU budget can either train a larger model with a slower convolution implementation, or train a smaller model with FlashFFTConv for more steps. The paper provides evidence that the latter can win.

3. Length generalization for pretrained models. The paper introduces partial convolutions — zeroing out later portions of the convolution kernel — as a mechanism to extend pretrained models to longer sequences without retraining. This matters because retraining a model from scratch for a new sequence length is cost-prohibitive. Partial convolutions provide a "sliding window" approach (Section 3.3) where a model trained at, say, 1 million sequence length can be applied to 4 million sequence length inputs with minimal quality degradation — the perplexity of HyenaDNA-1M drops only from 2.91 to 2.90 when extended to 4M (Table 8). This capability is purely a consequence of the efficiency gains from FlashFFTConv enabling the memory footprint reduction to be practically useful.

Where Prior Approaches Fall Short

The paper identifies two specific bottlenecks that existing FFT convolution implementations fail to address. Understanding these requires some background on GPU architecture, which the paper provides (Section 2.2, Figure 1 left):

Bottleneck 1: FFTs don't use tensor cores effectively. Modern GPUs (since the NVIDIA V100) contain specialized matrix-matrix multiply units called tensor cores. These units are dramatically faster than general-purpose compute units for their target operation — the paper cites the H100 tensor core achieving 1.0 PetaFLOP/s for matrix multiplication compared to 67 TeraFLOP/s for general arithmetic (Section 2.2). This is roughly a $15\times$ difference in peak throughput.

The traditional FFT convolution pipeline works as follows (Equation 1, Section 2.1):

Compute the FFT of the input: $\mathcal{F}u$
Compute the FFT of the kernel: $\mathcal{F}k$ (often precomputed, since kernels are shared across batches)
Compute elementwise multiplication in frequency space: $\mathcal{F}u \odot \mathcal{F}k$
Compute the inverse FFT: $\mathcal{F}^{-1}(\mathcal{F}u \odot \mathcal{F}k)$

Steps 1 and 4 involve butterfly operations — recursive combinations of complex multiplications and additions that recursively decompose the sequence into smaller pieces and recombine them. These operations are fundamentally not matrix multiplications at the scale that tensor cores are designed for. As a result, existing FFT implementations (including PyTorch's, which wraps highly tuned libraries like cuFFT) rely primarily on the GPU's general-purpose compute units, leaving tensor cores severely underutilized. The paper states this explicitly: "FFT convolutions do not effectively use the specialized matrix-matrix multiply units available on modern accelerators" (Section 1, paragraph 5).

This is not a failure of implementation quality — cuFFT is heavily optimized — but rather a mismatch between the algorithm's computational structure and the hardware's specialized throughput. The butterfly network of an FFT requires data-dependent addressing patterns and operations on small sub-vectors that do not naturally batch into the $16 \times 16$ matrix multiply operations that tensor cores accelerate.

Bottleneck 2: I/O costs dominate at long sequence lengths. The GPU memory hierarchy (Figure 1 left, Section 2.2) has three levels with dramatically different characteristics:

Level	Typical Size (H100)	Bandwidth
HBM (global memory)	40-80 GB	~1.5 TB/s
SRAM (shared memory)	~64 KB per SM	~19 TB/s
Registers	Tiny per thread	Fastest

For short sequences, data fits in SRAM, and operations can be fused — multiple computations are performed on data while it resides in fast memory, without writing intermediate results back to slow HBM. This kernel fusion is "common (and can be automated) for pointwise operations" (Section 2.2), such as the elementwise multiplication in step 3 of the FFT convolution, or additive gating operations commonly used in convolutional language models.

The critical problem emerges as sequence length increases: "sequences become too large to fit in SRAM, and kernel fusion fails, resulting in expensive I/O costs" (Section 1, paragraph 5). When intermediate FFT results exceed the 64 KB SRAM capacity, they must be written to HBM and read back later. HBM bandwidth (~1.5 TB/s) is roughly an order of magnitude slower than SRAM bandwidth (~19 TB/s). For long sequences, these I/O operations — not the floating-point arithmetic — become the dominant cost, turning a compute-bound operation into a memory-bound one.

The paper identifies two exacerbating factors for this I/O bottleneck:

Padding for causality: For causal convolutions (where output position $i$ can only depend on input positions $\leq i$ , as required in autoregressive language modeling), the input and kernel must be zero-padded to twice the sequence length before the FFT. This padding operation involves reading the input from HBM, extending it with zeros, writing the result to HBM, and reading it again for the FFT — pure I/O overhead that does no useful computation.
Real-to-complex conversion: The FFT convolution mathematically operates on complex numbers — the frequency-domain representations $\mathcal{F}u$ and $\mathcal{F}k$ are complex-valued even when $u$ and $k$ are real. Storing and transferring complex numbers ( $2 \times$ the size of real numbers) doubles the I/O cost through the memory hierarchy, and the conversion between real inputs and complex intermediates adds further overhead.

The fundamental tension: These two bottlenecks create a cruel trade-off. Using larger matrix multiply operations (lower-order Monarch decomposition) increases tensor core utilization (addressing Bottleneck 1) but requires larger matrices that consume more SRAM, triggering Bottleneck 2 at longer sequence lengths. Using smaller matrix multiply operations (higher-order decomposition) keeps SRAM usage low and enables kernel fusion for longer sequences, but reduces tensor core utilization because the matrices become too small to effectively fill the tensor core units. The paper's central technical contribution — the cost model for selecting the order $p$ of the Monarch decomposition — directly addresses this tension by characterizing precisely where these trade-off points occur.

Prior Work on FFT Optimization and Why It's Insufficient

The paper positions itself within a rich history of FFT algorithm design while arguing that existing approaches fail to address the specific challenges of modern GPU hardware for machine learning workloads:

Classical FFT algorithms (Cooley-Tukey and variants, Section 5). The Cooley-Tukey algorithm (1965) established the fundamental divide-and-conquer structure that all modern FFT implementations build on. Bailey's FFT algorithm (1989) specifically addressed the memory hierarchy problem for external or hierarchical memory — which is conceptually similar to the GPU memory hierarchy problem FlashFFTConv tackles. However, these classical algorithms were designed for CPU architectures where the primary bottleneck was the gap between main memory and cache, and where specialized matrix multiply units did not exist. They do not account for tensor cores at all, and their data layout strategies are optimized for cache line sizes rather than for warp-level tensor core operations.

NVIDIA's cuFFT and cuFFTdx. cuFFT is NVIDIA's highly optimized FFT library, and cuFFTdx is a more recent library (referenced in Section 4.2) that "recovers the strong baseline" of kernel fusion for FFTs. The paper's ablation study in Table 3 includes a "Fusion-Only/cuFFTdx" row, which represents using kernel fusion without the Monarch decomposition for tensor cores. This baseline achieves meaningful speedup over PyTorch but has two critical limitations:

It cannot use tensor cores for the matrix multiply operations, since it relies on the standard butterfly decomposition that maps to general-purpose arithmetic units.
It "does not support sequences longer than 32K due to a lack of SRAM space" (Section 4.2, Table 3 commentary), because without the Monarch decomposition's ability to decompose the FFT into smaller SRAM-resident blocks, the entire working set must fit in SRAM to remain fused.

The comparison makes the paper's contribution clear: kernel fusion alone is necessary but insufficient. The Monarch decomposition is what enables both tensor core utilization and SRAM-friendly working sets simultaneously.

Prior tensor core FFT work (Section 5). The paper cites tcFFT (Li et al., 2021), which explored using tensor cores for half-precision FFTs, and notes that "our work continues a line of work exploring how to use tensor cores for the FFT convolution [43, 44, 69], and extends the algorithmic capabilities to much longer sequences." The key distinction is scope: prior tensor core FFT work focused on the FFT operation itself, for relatively short transforms that fit in SRAM. FlashFFTConv integrates this with the full convolution pipeline (FFT → pointwise multiply → iFFT) and handles the memory hierarchy challenges that emerge at the million-length sequences required for DNA modeling and high-resolution vision.

The Monarch matrix decomposition (Dao et al., 2022, Section 5). The paper builds directly on the Monarch decomposition framework, which shows that the FFT matrix can be factorized into a product of structured sparse matrices that are composed primarily of block-diagonal matrix multiplies and permutations. The Monarch paper (ICML 2022) established the theoretical expressiveness and trainability of these structured matrices. FlashFFTConv's contribution is applying this decomposition specifically for hardware-efficient inference of convolutions, including the critical adaptation of broadcasting matrix operations across the sequence dimension rather than the batch dimension to enable fusion for longer sequences (Section 3.1, Figure 3).

How FlashFFTConv Positions Itself

The paper explicitly frames itself as doing for convolutions what FlashAttention did for attention mechanisms (Section 1, paragraph 3):

"Just as systems advances such as FlashAttention yielded improvements in modeling quality [1, 70] and the development of new attention algorithms [2, 66, 73, 92], we hope that understanding how to optimize the FFT convolution can also inspire algorithmic innovation, thus improving the quality of convolutional sequence models."

This positioning carries specific implications:

1. It frames the work as infrastructure, not an architecture. FlashFFTConv is not a new convolution variant, a new state-space model, or a new way to parameterize kernels. It is an implementation technique that accelerates an existing mathematical operation (the FFT convolution) without changing its output. This means any model that uses FFT convolutions — Hyena, M2-BERT, S4 variants, CKConv, etc. — can benefit from FlashFFTConv without architectural modification. The paper demonstrates this breadth by benchmarking across five different model families spanning four modalities and four orders of magnitude in sequence length (Table 5).

2. It claims that efficiency improvements enable qualitatively new capabilities, not just faster training. The paper's most striking results are not the speedup numbers themselves but the downstream consequences: solving Path-512 for the first time (from 50% to 96.1% accuracy), embedding the longest human genes, and achieving quality improvements equivalent to doubling model parameters. This echoes how FlashAttention's memory efficiency enabled longer-context Transformers and spurred the development of new attention algorithms — the efficiency gain changes what is possible, not just what is affordable.

3. It introduces a new design axis: sparsity in convolution filters. Beyond the core optimization, the paper proposes partial convolutions and frequency-sparse convolutions (Section 3.3) as architectural modifications that map naturally onto the Monarch decomposition's compute model. These are presented as "analogues to sparse/approximate attention in Transformers" (Section 1, paragraph 8), drawing a parallel to the extensive literature on sparse attention mechanisms. The key insight is that these sparsity patterns are not merely parameter-count reductions — they directly reduce runtime in FlashFFTConv because they can be implemented by "skipping blocks in the matrix decomposition" (Abstract), translating structural sparsity into wall-clock speedup. This contrasts with naive sparsity implementations where zeroing out weights does not necessarily yield runtime improvement due to the overhead of sparse matrix operations.

The Unstated Assumption Worth Noting

Throughout the paper, there is an implicit bet: that the FFT convolution will remain a dominant primitive for sequence modeling, making it worth the substantial engineering investment to optimize it to the hardware limit. The authors do not entertain the possibility that some other $O(N \log N)$ or $O(N)$ sequence mixing primitive (e.g., state-space models that avoid the FFT entirely, or novel linear attention mechanisms) might render the FFT convolution obsolete. Instead, they cite the widespread adoption of convolutions across "language modeling, time-series analysis, computer vision, DNA modeling, and more" (Section 1, paragraph 1) as evidence that the primitive is here to stay. This assumption justifies the depth of the engineering — FlashFFTConv involves custom CUDA kernels tuned for specific sequence lengths, not a generic library — and represents a bet that the field's architecture trajectory will continue to include FFT-based convolutions as a core building block.

3. Technical Approach

3.1 Reader Orientation

FlashFFTConv is a GPU kernel implementation that computes the FFT convolution — the mathematical operation $y = u * k$ where both the input $u$ and the kernel $k$ are sequences of length $N$ — by decomposing the Fast Fourier Transform into a series of matrix-matrix multiply operations that run efficiently on tensor cores, while fusing intermediate operations to avoid expensive memory transfers. The system solves the fundamental hardware mismatch that makes FFT convolutions slow despite their $O(N \log N)$ asymptotic complexity: standard FFT implementations use general-purpose arithmetic units (leaving tensor cores idle) and require writing intermediate results to slow global memory for sequences too long to fit in fast on-chip SRAM, whereas FlashFFTConv restructures the computation so that both problems are addressed simultaneously through a single decomposition whose order $p$ can be tuned to the sequence length.

3.2 Big-Picture Architecture (Diagram in Words)

The FlashFFTConv system has five major components that transform an input sequence into a convolved output through frequency-domain multiplication:

Input and Kernel Preparation — The real-valued input $u \in \mathbb{R}^{B \times H \times N}$ (batch size $B$ , hidden dimension $H$ , sequence length $N$ ) and complex-valued kernel $k_f \in \mathbb{C}^{H \times N}$ (the pre-computed FFT of the convolution filter, reused across batch items) are loaded into GPU memory. For causal convolutions, implicit zero-padding is handled without an explicit memory copy.
Real-to-Complex Conversion via Decimation-in-Time — Rather than running a complex FFT of length $N$ on zero-padded real data, a classical signal-processing algorithm transforms the real input into a complex sequence of length $N/2$ , halving the FFT cost. The even-indexed samples go into the real part and odd-indexed samples into the imaginary part of a complex vector $z(n)$ .
Order-p Monarch FFT Decomposition — The core computational engine: the $N$ -point FFT matrix $\mathbf{F}_N$ is factorized into a product of $p$ structured matrices, each of which is a block-diagonal matrix multiply (efficiently executed on tensor cores) followed by a permutation (implemented as an in-SRAM matrix transpose). The order $p$ controls the trade-off: higher $p$ means smaller individual matrix multiplies (lower FLOPs but more I/O between SRAM and registers) while lower $p$ means larger matrix multiplies (higher tensor core utilization but greater SRAM pressure).
Frequency-Domain Pointwise Multiplication — The transformed input $\mathcal{F}u$ is multiplied elementwise with the pre-computed kernel $k_f$ . In the Monarch decomposition, this multiplication is performed directly in the accumulator registers between the forward and inverse FFT, avoiding a round-trip through SRAM.
Inverse Monarch FFT Decomposition — The same order-p decomposition applied in reverse transforms the frequency-domain product back to the time domain, producing the convolved output $y \in \mathbb{R}^{B \times H \times N}$ . Intermediate results are recomputed in the backward pass rather than stored, trading additional compute for reduced memory footprint.

Information flows sequentially: input $\rightarrow$ decimation-in-time $\rightarrow$ forward Monarch FFT (outermost matrix multiply $\rightarrow$ twiddle correction $\rightarrow$ inner decomposition layers with fusion) $\rightarrow$ pointwise multiply with $k_f$ $\rightarrow$ inverse Monarch FFT (matching decomposition order) $\rightarrow$ decimation-in-time inverse $\rightarrow$ output. For gated convolutions (common in Hyena and M2 architectures), elementwise multiplicative gating operations are fused into this pipeline before the FFT and after the iFFT, eliminating separate kernel launches.

3.3 Roadmap for the Deep Dive

First, the FFT convolution algorithm and why it resists hardware-efficient implementation — this establishes what FlashFFTConv must overcome and explains the key properties of the GPU memory hierarchy that any solution must respect.
Second, the Monarch decomposition of the FFT matrix — the mathematical factorization that turns butterfly operations into block-diagonal matrix multiplies, which is the enabling insight for tensor core utilization. We examine how the decomposition order $p$ controls matrix sizes and why the decomposition broadcasts across the sequence dimension (not batch/hidden) to enable fusion.
Third, the GPU cost model for order-p decompositions (Equation 2) — the quantitative framework that predicts whether a given order $p$ will be compute-bound or memory-bound at a given sequence length, and the heuristic for selecting $p$ based on hardware constants. This is the intellectual core connecting the mathematical decomposition to wall-clock performance.
Fourth, kernel fusion and recomputation strategies — how the decomposition enables keeping intermediate data in SRAM (avoiding HBM trips) for sequences up to 32K on A100/H100, and how recomputation in the backward pass trades FLOPs for memory savings. This includes the handling of twiddle factors, padding, and the real-valued FFT algorithm.
Fifth, domain-specific optimizations — the decimation-in-time algorithm that halves FFT length by exploiting real-valued inputs, implicit zero-padding that eliminates outermost matrix multiply operations for causal convolutions, and gating fusion that avoids separate kernel launches for elementwise operations.
Sixth, architectural extensions: partial and frequency-sparse convolutions — how zeroing out portions of the convolution kernel (in time or frequency domain) maps directly to skipping blocks in the Monarch matrix decomposition, turning structural sparsity into actual runtime speedup.

3.4 Detailed, Sentence-Based Technical Breakdown

This is primarily a systems paper whose core idea is that the FFT convolution can achieve near-peak hardware utilization on modern GPUs by restructuring the computation into tensor-core-friendly matrix multiplies whose sizes are tuned to the sequence length, and by exploiting the resulting block-diagonal structure to keep working sets resident in fast on-chip memory.

The FFT Convolution and Its Hardware Mismatch

The convolution of two sequences $u$ (the input) and $k$ (the kernel) is defined as:

$(u * k)[i] = \sum_{j} u[j] \cdot k[i - j]$

where the sum runs over all positions $j$ where both sequences are defined. Computing this directly requires $O(N \cdot N_k)$ operations for sequence length $N$ and kernel length $N_k$ . When the kernel is as long as the input ( $N_k = N$ , which is standard in modern convolutional sequence models like Hyena and S4), direct computation costs $O(N^2)$ — prohibitive for the million-length sequences required in genomics and high-resolution vision.

The FFT convolution algorithm exploits the Convolution Theorem, which states that convolution in the time domain equals pointwise multiplication in the frequency domain:

$u * k = \mathcal{F}^{-1}(\mathcal{F}u \odot \mathcal{F}k)$

where $\mathcal{F}$ is the Discrete Fourier Transform operator, $\mathcal{F}^{-1}$ is its inverse, and $\odot$ denotes elementwise (Hadamard) multiplication. The FFT computes $\mathcal{F}$ in $O(N \log N)$ time, reducing the total convolution cost from $O(N^2)$ to $O(N \log N)$ .

What this equation computes: Given an input sequence $u$ and a kernel $k$ , the pipeline (1) transforms both to the frequency domain via FFT, (2) multiplies them pointwise (each frequency bin of $\mathcal{F}u$ is multiplied by the corresponding bin of $\mathcal{F}k$ ), and (3) transforms the result back via inverse FFT. The output is the circular convolution of $u$ and $k$ , where indices wrap around modulo $N$ . For causal convolutions (where output position $i$ cannot depend on input positions $> i$ ), both sequences are zero-padded to at least $2N$ before the FFT to prevent wrap-around interference.

Why this form: The FFT convolution decomposes an $O(N^2)$ operation into three $O(N \log N)$ operations plus one $O(N)$ pointwise multiply. The $\log N$ factor comes from the recursive divide-and-conquer structure of the Fourier transform, where a length- $N$ transform is computed from two length- $N/2$ transforms. This is asymptotically efficient but, critically, the standard implementation of this divide-and-conquer (the Cooley-Tukey butterfly) does not map to the tensor core matrix-multiply primitive that provides the highest GPU throughput.

The hardware mismatch in detail. A GPU tensor core multiplies two $16 \times 16$ matrices (or $32 \times 8 \times 16$ shapes) and accumulates the result — it is designed for dense linear algebra with predictable, regular memory access patterns. The Cooley-Tukey FFT butterfly, by contrast, involves:

Strided memory access: At stage $s$ of the decomposition, the butterfly combines elements separated by $2^s$ positions — a non-contiguous access pattern that prevents coalesced memory reads.
Small arithmetic intensity: Each butterfly node computes one complex multiplication and one complex addition, reading two complex numbers, performing two floating-point operations per real component, and writing two complex numbers. The ratio of FLOPs to bytes is too low to keep tensor cores fed.
Bit-reversal permutation: The final step reorders outputs into natural order, requiring a global data shuffle that does no arithmetic at all — pure I/O.

These properties mean that cuFFT and PyTorch's FFT implementations run primarily on the GPU's general-purpose CUDA cores (67 TFLOP/s on H100) rather than tensor cores (1,000 TFLOP/s on H100) — a potential $15\times$ throughput gap.

The memory hierarchy bottleneck. The GPU memory hierarchy (Section 2.2, Figure 1) has three levels: HBM (global memory, ~40 GB, ~1.5 TB/s bandwidth), SRAM (shared memory, ~64 KB per streaming multiprocessor, ~19 TB/s bandwidth), and registers (per-thread, fastest). For short sequences (up to ~2K on A100), the entire working set of the FFT fits in SRAM, enabling kernel fusion: the FFT, pointwise multiply, and inverse FFT can be executed as a single GPU kernel that loads input from HBM once, performs all computations in SRAM and registers, and writes the final output back to HBM once. Each intermediate trip to HBM is eliminated.

For long sequences, the working set exceeds SRAM capacity. A length- $N$ complex FFT requires $O(N)$ intermediate storage. When $N$ is large enough that $O(N)$ exceeds 64 KB, the implementation must spill to HBM — writing intermediate butterfly results, reading them back for the next stage, writing the final FFT output, reading it for the pointwise multiply, writing that result, and reading it for the inverse FFT. Each HBM round-trip costs $\sim$ 1.5 TB/s of bandwidth, and for $N = 1\text{M}$ , the total data movement can exceed the arithmetic cost by an order of magnitude, turning a compute-bound operation into a memory-bound one.

The two bottlenecks are coupled. The standard FFT implementation suffers from both poor compute unit utilization (wasting tensor cores) AND poor memory hierarchy utilization (excessive HBM traffic for long sequences). A solution that addresses only one — e.g., kernel fusion without using tensor cores (the cuFFTdx baseline in Table 3) — still leaves the other as a bottleneck and cannot support sequences beyond 32K because the fused kernel's working set eventually exceeds SRAM.

The Monarch FFT Decomposition

The core mathematical insight of FlashFFTConv is that the $N \times N$ Discrete Fourier Transform matrix $\mathbf{F}_N$ can be factorized into a product of structured matrices, each composed primarily of block-diagonal matrix-matrix multiplies. This is the Monarch decomposition, introduced by Dao et al. (2022) for training efficient neural networks and here repurposed for hardware-efficient inference.

Order-2 Monarch decomposition. Let $N = N_1 \cdot N_2$ be a factorization of the sequence length into two factors. The order-2 Monarch decomposition factorizes the Fourier matrix as:

$\mathbf{F}_N = \mathbf{P} (\mathbf{I}_{N_2} \otimes \mathbf{F}_{N_1}) \mathbf{D} \mathbf{P}^{-1} (\mathbf{I}_{N_1} \otimes \mathbf{F}_{N_2}) \mathbf{P}$

where:

$\mathbf{F}_N$ is the $N \times N$ discrete Fourier matrix (element $(j, k)$ is $e^{-2\pi i j k / N}$ ),
$\mathbf{F}_{N_1}$ and $\mathbf{F}_{N_2}$ are smaller Fourier matrices of sizes $N_1 \times N_1$ and $N_2 \times N_2$ ,
$\mathbf{I}_{N_2} \otimes \mathbf{F}_{N_1}$ is the Kronecker product that replicates $\mathbf{F}_{N_1}$ as $N_2$ independent blocks along the diagonal of a larger matrix — a block-diagonal matrix where each $N_1 \times N_1$ block is an independent FFT,
$\mathbf{P}$ is a permutation matrix that reshapes the input vector of length $N$ into an $N_1 \times N_2$ matrix, transposes it to $N_2 \times N_1$ , and reshapes back to a vector of length $N$ ,
$\mathbf{D}$ is a diagonal matrix of "twiddle factors" — complex corrections $e^{-2\pi i j k / N}$ that account for the phase shifts introduced by decomposing the full FFT into smaller independent FFTs.

What this factorization computes: The FFT of a length- $N$ vector is computed in five steps, illustrated in Figure 2:

Reshape and transpose: The input vector of length $N$ is reshaped into an $N_1 \times N_2$ matrix (reading row-major), then transposed to $N_2 \times N_1$ . This is the $\mathbf{P}$ permutation.
FFT on columns: An independent FFT of length $N_2$ is applied to each of the $N_1$ columns. This is $\mathbf{I}_{N_1} \otimes \mathbf{F}_{N_2}$ : $N_1$ parallel $N_2$ -point FFTs. In the block-diagonal matrix view, the $N \times N$ matrix is block-diagonal with $N_1$ blocks, each of size $N_2 \times N_2$ .
Twiddle factor correction: Each element is multiplied by a complex correction factor $\mathbf{D}$ . This is a pointwise (diagonal) operation.
Inverse permutation: The reverse reshape-transpose is applied ( $\mathbf{P}^{-1}$ ), returning to $N_2 \times N_1$ layout.
FFT on columns: An independent FFT of length $N_1$ is applied to each of the $N_2$ columns. This is $\mathbf{I}_{N_2} \otimes \mathbf{F}_{N_1}$ : $N_2$ parallel $N_1$ -point FFTs.
Final permutation: $\mathbf{P}$ restores the original ordering.

Why this form matters for hardware: Each $\mathbf{I} \otimes \mathbf{F}$ term is a block-diagonal matrix — a large matrix where the only non-zero entries are in $N_1 \times N_1$ (or $N_2 \times N_2$ ) blocks along the diagonal. Multiplying by a block-diagonal matrix is equivalent to performing independent matrix-matrix multiplications on each block. When these blocks are large enough (at least $16 \times 16$ , the tensor core tile size), each block multiplication can be executed on a tensor core at the full 1.0 PFLOP/s rate. The permutation operations $\mathbf{P}$ and $\mathbf{P}^{-1}$ , when the computation is broadcast across the sequence dimension (as explained below), become simple matrix transposes that can be done by renaming indices in SRAM without any data movement.

The critical transformation is that the FFT — originally a recursive butterfly network with irregular memory access — has been rewritten as a sequence of (block-diagonal matrix multiplies + transposes + pointwise multiplies), where the matrix multiplies dominate the FLOP count and are perfectly suited for tensor cores.

Higher-order decompositions. An order- $p$ Monarch decomposition recursively applies the order-2 factorization to the inner FFT matrices. If $N = \prod_{i=1}^p N_i$ (where each $N_i$ is typically chosen to be equal for balanced decomposition, so $N = N_1^p$ ), the FFT is expressed as a product of $p$ matrix operations, each of which is a block-diagonal matrix with blocks of size $N_i \times N_i$ . Specifically:

$p = 2$ : The FFT is decomposed into blocks of size $\sqrt{N} \times \sqrt{N}$ . For $N = 1\text{M}$ , blocks are $1024 \times 1024$ — large matrix multiplies that use tensor cores efficiently, but the intermediate working set (the full $\sqrt{N} \times \sqrt{N}$ matrix) is $\sim$ 1M elements, exceeding SRAM.
$p = 3$ : Blocks are $\sqrt[3]{N} \times \sqrt[3]{N}$ . For $N = 1\text{M}$ , blocks are $100 \times 100$ — smaller, fitting in SRAM but below optimal tensor core utilization (since tensor cores prefer multiples of 16).
$p = 4$ : Blocks are $\sqrt[4]{N} \times \sqrt[4]{N}$ . For $N = 1\text{M}$ , blocks are $32 \times 32$ — near-optimal for tensor cores and with a much smaller SRAM footprint, but requiring more permutation (transpose) operations between layers, increasing I/O between SRAM and registers.

The order $p$ controls a fundamental trade-off: higher $p$ reduces the size of each matrix multiply (lower FLOPs, since the total FLOPs of an order-p Monarch FFT is $O(N^{(p+1)/p})$ ) but increases the number of permutation/transpose steps (more I/O between decomposition layers). This trade-off is formally characterized by the cost model in the next section.

Why broadcast over sequence, not batch/hidden. Traditional FFT implementations — including classical parallel FFT algorithms and earlier Monarch-based implementations — parallelize by assigning different batch elements or hidden dimensions to different compute units. Each streaming multiprocessor computes the full FFT for one or a few batch items independently. This is natural because the FFT operates along the sequence dimension, and sharing the sequence across compute units would require communication.

FlashFFTConv inverts this: it broadcasts the matrix multiply operations across the sequence dimension and parallelizes across batch and hidden dimensions (Figure 3, top panel). The rationale is subtle but essential for kernel fusion:

When an $N_1 \times N_1$ matrix multiply is parallelized across batch items, each SM processes the entire $N_1 \times N_1$ block for its assigned batch item. To fuse multiple operations (e.g., the forward FFT, pointwise multiply with $k_f$ , and inverse FFT) without writing to HBM, the SM must hold the entire $N_1 \times N_1$ block for a given batch item in SRAM for the duration of all fused operations. This requires the block to be small enough to fit in SRAM alongside the kernel weights and twiddle factors.
When instead the $N_1 \times N_1$ matrix multiply is parallelized across the sequence (i.e., different SMs handle different portions of the input sequence for the same batch item), each SM only needs to hold a fraction of the sequence in SRAM — specifically, the rows or columns of the $N_1 \times N_1$ block that it is responsible for. This reduces the per-SM SRAM requirement, allowing kernel fusion to remain viable at longer sequence lengths. The paper states that broadcasting along the sequence "reduces the SRAM requirements for kernel fusion, since we only need to load a single sequence into SRAM at a time — allowing us to fuse the entire kernel for sequences up to 32K on A100 and H100" (Section 3.1, "Adapting Monarch for Fusion").
An additional benefit: when the computation is structured this way, the permutations $\mathbf{P}$ in the Monarch decomposition (which would otherwise require explicit data shuffling) become simple matrix transposes. The data for one SM, organized as a tile of the $N_1 \times N_2$ matrix, is transposed by simply reading it in column-major order instead of row-major order — a zero-cost operation in SRAM using established on-chip transpose routines.

Algorithm 1 (core loop). The complete FlashFFTConv algorithm for order-2 is given in Algorithm 1. The nested loop structure is:

Outer loop: Across SMs, tiled by $B_{\text{tile}} \times H_{\text{tile}}$ (batch tiles and hidden dimension tiles). Each SM loads the FFT matrices $\mathbf{F}$ and $\mathbf{F}^{-1}$ (the $N_1 \times N_1$ block FFT matrices, which are shared across all batch items) and the twiddle factors $\mathbf{t}$ and $\mathbf{t}_{\text{inv}}$ once from HBM.
Middle loop: Across hidden dimensions within the tile. For each $h$ , the kernel $\mathbf{K}_f \leftarrow k_f[h]$ is loaded and reshaped to $N_1 \times N_1$ . This kernel is reused across all batch items in the tile, amortizing the HBM load cost.
Inner loop: Across batch items within the tile. For each $b$ , the input $\mathbf{X} \leftarrow u[b,h]$ is loaded, reshaped to $N_1 \times N_1$ .

The core computation for each batch item is a sequence of six operations performed entirely in registers (with SRAM used only for the transpose between the forward and inverse FFT):

$\mathbf{X} \leftarrow \mathbf{F}^\top \mathbf{X}$ — forward FFT first step (matrix multiply, tensor core)
$\mathbf{X} \leftarrow \mathbf{X} * \mathbf{t}$ — twiddle factor correction (elementwise, fused in register)
$\mathbf{X} \leftarrow \mathbf{X} \mathbf{F}$ — forward FFT second step (matrix multiply, tensor core)
$\mathbf{X} \leftarrow \mathbf{X} * \mathbf{K}_f^\top$ — pointwise multiply with kernel in frequency domain (elementwise, fused in register)
$\mathbf{Y} \leftarrow ((\mathbf{X} \mathbf{F}^{-1})^\top * \mathbf{t}_{\text{inv}}) \mathbf{F}^{-1}$ — inverse FFT (two matrix multiplies with a twiddle correction and transpose)
Write $\mathbf{Y}^\top$ to HBM

The key implementation detail is that operations 1-4 execute without writing intermediate results to SRAM — the output of each matrix multiply is held in the accumulator register fragment and directly reused as the input operand for the next operation (Section A.2, "Register Reuse"). Only the transpose between the forward FFT and the twiddle correction for the inverse FFT requires a trip through SRAM, since the data layout must change from row-major to column-major (or vice versa) for the next matrix multiply.

Generalization to higher orders (Algorithms 3 and 4). For $p=3$ (Algorithm 3), an outer loop iterates $N_1$ times over the rows of the $N_1^3 = N$ reshaped input. Each iteration performs an inner order-2 decomposition (two matrix multiplies with twiddle correction, pointwise multiply with the corresponding slice of $\mathbf{K}_f$ , and two inverse matrix multiplies). The intermediate results between iterations can remain in SRAM if $N$ is small enough, making this a fully fused kernel.

For $p=4$ (Algorithm 4), the outermost level treats the 4-way decomposition as: one matrix multiply for the outermost FFT step, a call to the fully fused 3-way decomposition as a subroutine, and one matrix multiply for the outermost inverse FFT step. Intermediate results between these three stages are written to HBM because the total working set exceeds SRAM capacity — this is the regime where kernel fusion is partial: the inner decomposition is fused, but the outermost layers require HBM I/O.

The real-valued FFT optimization (decimation-in-time). The convolution used in ML operates on real-valued inputs and produces real-valued outputs (the kernel weights are real numbers). The standard FFT, however, operates on complex numbers — if you feed a real sequence of length $N$ into a complex FFT of length $N$ , you are computing with twice the necessary precision and bandwidth.

FlashFFTConv uses a classic algorithm called one-stage decimation in time (Appendix A.1, following Sorensen et al., 1987) to compute a real-to-real FFT of length $N$ using a complex FFT of length $N/2$ . The procedure is:

Pack the real input into a complex vector of half length: Create $z(n) = x(2n) + i \cdot x(2n+1)$ for $n = 0, \ldots, N/2-1$ , where $x$ is the real input sequence. Even-indexed samples go into the real part; odd-indexed samples go into the imaginary part.
Compute the complex FFT of length $N/2$ : $Z(k) = \text{FFT}_{N/2}(z)$ .
Recover the even and odd parts of the full FFT: From the Hermitian symmetry of real-valued FFTs, the even-indexed frequency bins $X_e[k]$ and odd-indexed frequency bins $X_o[k]$ of the full $N$ -point FFT can be recovered from $Z[k]$ using: $X_e[k] = \frac{Z[k] + Z^*[N/2 - k]}{2}, \quad X_o[k] = -i \cdot \frac{Z[k] - Z^*[N/2 - k]}{2i}$ where $Z^*$ denotes complex conjugation.
Combine to form the full FFT: $X[k] = X_e[k \bmod N/2] + X_o[k \bmod N/2] \cdot W_N^k$ , where $W_N^k = e^{-2\pi i k / N}$ is the twiddle factor.

The inverse procedure (real iFFT from complex iFFT of length $N/2$ ) reverses these steps.

What this achieves: The FFT cost is cut in half — from an $N$ -point complex FFT to an $N/2$ -point complex FFT, plus $O(N)$ bookkeeping operations. Since the $N/2$ -point FFT costs $O(\frac{N}{2} \log \frac{N}{2}) \approx \frac{1}{2} O(N \log N)$ , the total savings approach 50% for large $N$ . Critically, this optimization is applied as bookkeeping at the boundaries of the Monarch decomposition — the input is packed before entering the first matrix multiply, and unpacked after the final matrix multiply — so the core matrix multiply operations see complex data of length $N/2$ rather than $N$ .

Implicit zero-padding for causal convolutions. For causal convolutions, the input $u$ of length $N$ is padded with $N$ zeros to length $2N$ before the FFT, and the kernel $k$ is zero-padded similarly. In a naive implementation, this padding requires allocating a $2N$ -length buffer, copying the $N$ input elements, zero-filling the rest, and then loading the full $2N$ buffer for the FFT. The pad operation itself is pure I/O overhead.

FlashFFTConv eliminates this by recognizing that the zero-padded regions correspond to specific parts of the $N_1 \times N_1$ block structure in the Monarch decomposition. Specifically, for the outermost matrix multiply operations in the forward and inverse FFT, half of the matrix columns (or rows) correspond to the zero-padded region and can be entirely skipped — the corresponding matrix multiply operations produce zero outputs. This "eliminate[s] half of the outermost matrix multiply operations in the FFT and iFFT" (Section 3.1, "Domain-Specific Optimizations"), reducing FLOPs and memory traffic without any explicit padding kernel.

The GPU Cost Model for Order-p Decompositions

The paper introduces a cost model (Equation 2) that predicts the total time of an order-p Monarch convolution as a function of sequence length, decomposition order, and hardware parameters. This model is not used at runtime — it is an offline design tool for selecting $p$ given the target sequence length and GPU.

The cost model accounts for both compute and I/O, similar to a roofline analysis. The total cost $C$ is:

$C = BH \sum_{i=1}^{p} \left( \frac{16 N N_i}{\gamma(N_i)} + \frac{4N}{\omega(i)} \right)$

where:

$B$ is the batch size,
$H$ is the hidden dimension (number of channels),
$N$ is the sequence length, factored as $N = \prod_{i=1}^p N_i$ (typically $N_i = N^{1/p}$ for balanced decomposition),
$N_i$ is the block size at decomposition step $i$ ,
$\gamma(N_i)$ is the achievable FLOP rate for the operation at step $i$ , which depends on whether $N_i$ is large enough to use tensor cores: $\gamma(N_i) = \tau_M$ (tensor core FLOP rate, ~234 TFLOP/s on A100) if $N_i \geq \mu$ (where $\mu = 16$ , the tensor core tile size), and $\gamma(N_i) = \tau_G$ (general arithmetic FLOP rate, ~17.6 TFLOP/s on A100) if $N_i < \mu$ ,
$\frac{16 N N_i}{\gamma(N_i)}$ is the compute time for step $i$ : the FLOP count $O(N N_i)$ divided by the achievable throughput. The factor 16 accounts for the $16 \times 16$ tensor core tile size and complex arithmetic overhead,
$\omega(i)$ is the bandwidth of the memory where intermediate results at step $i$ are stored: $\omega(i) = \sigma_S$ (SRAM bandwidth, ~9.5 TB/s on A100) if the working set fits in SRAM, and $\omega(i) = \sigma_H$ (HBM bandwidth, ~1.35 TB/s on A100) if it requires HBM,
$\frac{4N}{\omega(i)}$ is the I/O time for step $i$ : reading and writing $4N$ bytes (complex values, 4 bytes per real/imaginary half-precision component) at the achievable bandwidth.

What this equation computes: For each step $i$ of the $p$ -step decomposition, we sum (a) the time to compute $N$ independent $N_i \times N_i$ matrix multiplies, each requiring $O(N_i^3)$ FLOPs but with $N/N_i$ such multiplies at this step giving $O(N N_i)$ total FLOPs, divided by the relevant throughput $\gamma(N_i)$ , and (b) the time to read the input and write the output for this step, $O(N)$ data movement, divided by the relevant bandwidth $\omega(i)$ . The total is multiplied by $BH$ since the convolution is computed independently for each batch item and hidden channel.

Why this form: The $16$ and $4$ constants come from half-precision complex arithmetic accounting: each complex multiply-add involves 4 real multiplies and 4 real adds, and storing/loading a complex number in fp16 requires 4 bytes (2 for the real part, 2 for the imaginary part, each stored as a 16-bit float). The piecewise function $\gamma(N_i)$ captures the threshold effect: if a matrix multiply block is smaller than the tensor core tile size ( $16 \times 16$ for A100/H100), the tensor core cannot be used efficiently and the operation falls back to general-purpose CUDA cores at $\sim 13\times$ lower throughput.

Figure 4 interpretation. The cost model is plotted for $p \in \{2, 3, 4\}$ on A100, for sequence lengths from 256 to 4M (Figure 4). The y-axis is cost per token ( $C/N$ ), normalized — lower is better. Key observations:

At short sequences (256–1K): $p=4$ is actually more expensive than $p=3$ or $p=2$ , because the $N_i$ values are so small (e.g., $N_i = \sqrt[4]{256} = 4$ ) that the matrix multiplies are smaller than the tensor core tile size ( $\mu = 16$ ) and run at general arithmetic speed ( $\tau_G$ ). The bump in the $p=4$ curve corresponds to the "Matrices Too Small for Tensor Cores" region in Figure 4.
At medium sequences (4K–32K): $p=2$ is optimal — the blocks are large enough for tensor cores ( $N_i = \sqrt{N}$ ranges from ~63 to ~181) and the working set still fits in SRAM (enabling fully fused execution). $p=3$ and $p=4$ have smaller blocks that use tensor cores less efficiently and incur more I/O between decomposition layers.
At long sequences (64K–256K): $p=3$ becomes optimal. The $p=2$ working set ( $\sim \sqrt{N}$ elements) now exceeds SRAM capacity (the "SRAM Limit" region in Figure 4), forcing HBM spills, while $p=3$ blocks ( $\sim \sqrt[3]{N}$ ) still fit in SRAM. The $p=3$ curve shows a bump between 32K and 64K — this is the transition where $p=3$ also begins to spill to HBM.
At very long sequences (1M–4M): $p=4$ becomes necessary. The $p=3$ working set now also exceeds SRAM, and $p=4$ provides the smallest SRAM-resident blocks. Despite requiring more HBM I/O between outer decomposition layers, the inner fused kernels avoid the catastrophic HBM bandwidth bottleneck that would plague $p=2$ . The paper explicitly notes that for $p=4$ , intermediate results between the outermost matrix multiply and the inner 3-way decomposition are written to HBM (Appendix A.3, Algorithm 4).

Heuristic for selecting $p$ . The paper does not formulate a closed-form optimization over $p$ but rather uses this cost model to pre-select $p$ offline for each target sequence length. Table 3 shows the mapping: $p=2$ for sequences up to 32K, $p=3$ for 1M–2M, and $p=4$ for 4M (though the body text in Section 4.2 groups these into ranges). The selection is based on the hardware constants $\tau_M$ , $\tau_G$ , $\sigma_H$ , $\sigma_S$ , and $\mu$ , which the paper measures empirically (Appendix C, Table 19) rather than using theoretical peak values — this is important because achievable bandwidth and FLOP rates are typically 70–85% of theoretical peaks due to overhead.

Kernel Fusion and Recomputation Strategies

The Monarch decomposition enables two memory-saving strategies that are central to FlashFFTConv's efficiency: kernel fusion (keeping intermediate data in SRAM) and recomputation (trading FLOPs for memory in the backward pass).

Kernel fusion enabled by the Monarch decomposition. The key property that enables fusion is that inner layers of the decomposition do not require the entire sequence (Section 3.1, "Kernel Fusion and Recomputation"). In an order-p decomposition:

The outermost matrix multiply operates on $N/N_1 \times N_1$ blocks — it sees a fraction $1/N_1$ of the sequence at a time.
The next layer operates on $N/(N_1 N_2) \times N_2$ blocks — an even smaller fraction.
The innermost layers operate on $N_1 \times N_1$ blocks — the smallest working set.

For a fused kernel, all operations from the innermost layer up to the layer where the total SRAM requirement (blocks + weights + twiddle factors + kernel) exceeds SRAM capacity can be executed without HBM I/O. The outermost layers that do not fit write their intermediate results to HBM. This is a partial fusion strategy: the inner decomposition is fully fused, while the outer layers take explicit HBM I/O. Compare this to the standard FFT implementation, where every stage of the butterfly writes to HBM for long sequences.

Concretely, for $p=3$ at sequence length 1M, $N_1 = 100$ . The inner two matrix multiply operations (the order-2 decomposition within each iteration of Algorithm 3's outer loop) operate on $100 \times 100$ blocks — each block is $10,000$ complex elements $\times 4$ bytes = 40 KB, which fits comfortably in the 64 KB SRAM per SM. The outer loop over $N_1 = 100$ iterations can keep the working set for one or a few iterations in SRAM, processing the sequence in tiles. For $p=4$ at 4M, $N_1 = 32$ , and the inner 3-way decomposition (Algorithm 3) is fully fused with $32 \times 32$ blocks, while the outermost matrix multiply (Algorithm 4) writes its $32 \times 131072$ intermediate to HBM.

Recomputation in the backward pass. Deep learning training requires both the forward pass (computing the output $y$ ) and the backward pass (computing gradients with respect to the input $u$ and kernel $k$ ). The backward pass needs intermediate values from the forward pass — for example, the gradient through the pointwise multiply $\mathcal{F}u \odot k_f$ requires knowing $\mathcal{F}u$ . The standard approach is to store these intermediates in memory during the forward pass and read them during the backward pass — this is called "checkpointing" or "save for backward."

For long sequences, storing the full set of intermediates — $\mathcal{F}u$ (complex, length $N$ ), intermediate results between FFT stages — can consume enormous memory. FlashFFTConv instead recomputes these values during the backward pass: "Instead of storing intermediate results on HBM for the backward pass (e.g., the intermediate result of $\mathcal{F}u$ ), we simply recompute them in the backward pass" (Section 3.1, "Kernel Fusion and Recomputation").

What this trades: The forward pass is computed twice — once during the forward pass (to produce the output) and once during the backward pass (to produce the intermediates needed for gradient computation). This costs additional FLOPs (roughly 50% more forward-pass FLOPs) but saves memory proportional to the size of the intermediates. For the FFT convolution, these intermediates include the full complex-valued $\mathcal{F}u$ (size $N \times 2 \times 2$ bytes in fp16 = $4N$ bytes), intermediate Monarch decomposition results (additional $O(pN)$ storage), and twiddle-corrected values. At $N = 1\text{M}$ , the saved memory is on the order of tens of megabytes per batch item — enough to increase the maximum batch size by a factor of 2–4× (as shown in Tables 16–17 in Appendix B, where memory savings of 2.6× are reported for long sequences).

The paper's design choice is to always use recomputation rather than selectively checkpointing — a simpler strategy that avoids the complexity of deciding which intermediates to store and which to recompute. This is possible because the FFT convolution's arithmetic intensity (FLOPs per byte of intermediate storage) is relatively low, meaning that recomputation is cheaper than the HBM I/O that storing the intermediates would require.

Tiling across batch and hidden dimensions. The outer loops in Algorithm 1 tile across $B_{\text{tile}}$ and $H_{\text{tile}}$ (batch and hidden dimension tiles) to amortize the cost of loading the FFT matrices $\mathbf{F}$ and twiddle factors $\mathbf{t}$ from HBM. Since these matrices are shared across all batch items and hidden channels for a given sequence length, loading them once per tile of $B_{\text{tile}} \times H_{\text{tile}}$ items reduces HBM traffic proportionally. This is a standard GPU optimization, but it interacts with the sequence-length broadcasting: the tile sizes must be chosen so that the aggregate SRAM requirement (matrices $\mathbf{F}$ , twiddle factors, and the working sets for $B_{\text{tile}} \times H_{\text{tile}}$ active batch/channel items) does not exceed SRAM capacity.

Domain-Specific Optimizations

Beyond the core Monarch decomposition and fusion strategy, FlashFFTConv incorporates several optimizations specific to convolutional sequence models:

Gating operation fusion. Convolutional language models such as Hyena and M2-BERT use multiplicative gating: the convolution output is multiplied elementwise by a gate signal before being passed to the next layer. The operation is $y = v \odot ((u \odot w) * k)$ , where $v$ and $w$ are linear projections of the input $u$ (computed separately, usually with matrix multiplies that already run on tensor cores). In PyTorch, this requires three separate kernel launches: the convolution, the elementwise multiply with $w$ , and the elementwise multiply with $v$ . Each kernel launch incurs overhead, and each elementwise multiply reads the data from HBM and writes it back.

FlashFFTConv fuses the gating operations into the convolution kernel: the $u \odot w$ multiply is performed on the input before it enters the FFT pipeline, and the $v \odot$ multiply is performed on the output after the iFFT, all within the same kernel that computes the convolution. The data stays in registers or SRAM between these operations, avoiding two HBM round-trips. Table 4 quantifies the benefit: at sequence length 1K, FlashFFTConv achieves 7.93× speedup over PyTorch for gated convolutions, compared to 6.54× for standard convolutions (Table 3) — the additional speedup comes from eliminating the gating I/O overhead.

Register-level data management for tensor core operations. Appendix A.2 provides implementation detail on how matrix multiply results are held in registers between operations. CUDA tensor cores are accessed through the Warp Level Matrix Multiply Accumulate (WMMA) API, where operands are loaded into "fragments" (collections of registers with an unspecified mapping to threads). The critical optimization is that "the mapping of items to threads in the wmma::accumulator fragment exactly matches that for the wmma::matrix_a fragment read row-major, allowing us to directly copy the results of a matrix-matrix multiplication and use as the operand for another matrix-matrix multiply" (Appendix A.2). This means the output of one tensor core operation can be directly reused as the input to the next without writing to SRAM, as long as the data layout (row-major vs. column-major) matches. When the layout needs to change (e.g., for a transpose), an explicit SRAM round-trip is required.

Aggressive constant tuning. The paper states that they "aggressively tune our kernel hyperparameters such as block and tile dimensions, and loop unrolling factors for the best performance on the specific underlying hardware" (Appendix A.2). This includes double-buffering I/O (overlapping data movement with computation to hide latency) and using vector intrinsics for fp16/bf16 arithmetic to execute non-tensor-core operations at twice the normal throughput. These are standard CUDA optimization techniques, but their application is specific to each sequence length and decomposition order — FlashFFTConv generates specialized kernels per configuration rather than using a generic parameterized kernel.

Architectural Extensions: Partial and Frequency-Sparse Convolutions

The Monarch decomposition's structure — where the convolution is computed as a series of block-diagonal matrix multiplies — presents a natural interface for implementing sparsity in the convolution kernel. When portions of the kernel are zero, entire blocks in the matrix decomposition can be skipped, turning parameter sparsity into actual runtime speedup.

Partial convolutions (sparsity in time domain). In a partial convolution, later portions of the convolution kernel $k[t]$ for $t > T$ (where $T < N$ ) are set to zero. This means the convolution only looks at the most recent $T$ input positions for each output — analogous to local/sliding-window attention in Transformers.

Implementation in FlashFFTConv: Because the kernel $k_f$ in frequency domain is the FFT of the time-domain kernel $k$ , zeroing out the tail of $k$ corresponds to a specific structured modification of $k_f$ (not simply zeroing out a contiguous region). However, the benefit is simpler: for a kernel of effective length $T$ , the input $u$ only needs to be processed in chunks of length $T$ , and the sequence can be streamed through in overlapping windows. The memory footprint is reduced from $O(N)$ to $O(T)$ because only a window of the input needs to be in GPU memory at once. Table 7 quantifies this: for Hyena-s-8K with a 2K effective kernel length, the memory footprint drops from 32.5 GB to 11.8 GB.

Frequency-sparse convolutions (sparsity in frequency domain). In a frequency-sparse convolution, portions of $k_f$ (the frequency-domain kernel) are set to zero. This is the convolutional analogue of band-limiting or low-pass filtering: high-frequency components of the convolution are discarded. The sparsity pattern can be chosen to map efficiently onto the Monarch decomposition's block structure.

The paper describes a specific structured sparsity pattern (Appendix A.4): reshape the 2M-length kernel $k_f$ as $32 \times 32 \times 32 \times 64$ (corresponding to a $p=4$ decomposition with $N_1 = 32$ , except the last factor is 64). Sparsity is applied along each of the four dimensions sequentially: zero out dimensions $a, b, c, d$ where $a$ is the number of rows of the $32 \times 32 \times 32 \times 64$ tensor to zero, $b$ is columns within each row, etc. The sparsity fraction is:

$S = 1 - \frac{(32-a)(32-b)(32-c)(64-d)}{32 \cdot 32 \cdot 32 \cdot 64}$

Why this specific pattern: Each dimension's sparsity translates directly to skipped computation in the Monarch decomposition:

First dimension ( $a$ ): Allows skipping computation in the innermost matrix multiply (matrix $B$ in Appendix A.4). Each row eliminated removes one $N_1 \times N_1$ block multiply.
Second dimension ( $b$ ): Allows skipping computation in the second-innermost matrix multiply (matrix $A$ ).
Third dimension ( $c$ ): Reduces the number of iterations $\beta$ in the inner loop of the 3-way decomposition — each eliminated row skips one complete inner iteration.
Fourth dimension ( $d$ ): Reduces $\alpha$ , the number of iterations in the outer loop — each eliminated row skips an outer iteration.

This is a case where the sparsity pattern is designed to match the decomposition structure, not chosen arbitrarily. The result is that sparsity translates to actual speedup: Table 9 shows that at 79% sparsity, the convolution achieves 1.4× speedup, and even at 50% sparsity (1.2× speedup), perplexity does not degrade — it actually improves slightly (2.91 to 2.90), which the authors hypothesize may be due to "removing high-frequency noise" (Section 4.3).

The key design principle: Both partial and frequency-sparse convolutions are implemented by "skipping blocks in the matrix decomposition" (Abstract), meaning the structural sparsity of the kernel directly reduces the number of tensor core matrix multiply operations that must be executed. This is fundamentally different from unstructured weight sparsity, where zeroing individual weights rarely translates to speedup because the overhead of sparse matrix formats (indexing, irregular memory access) dominates the savings. By co-designing the sparsity pattern and the compute decomposition, FlashFFTConv achieves the elusive goal of "sparsity that actually makes things faster."

4. Key Insights and Innovations

Innovation 1: Reframing the FFT Convolution as a Hardware-Mapping Problem Rather than an Algorithmic Efficiency Problem

The paper's most fundamental conceptual move is not the Monarch decomposition itself — that was introduced in prior work (Dao et al., 2022) — but rather the diagnostic reframing that identifies why FFT convolutions are slow and what kind of solution is needed. Before this work, the dominant assumption in the sequence modeling community was that the FFT convolution's $O(N \log N)$ asymptotic complexity made it inherently efficient, and that any wall-clock gap relative to Transformers was due to implementation details that could be resolved with better engineering of the standard butterfly algorithm. This assumption was reinforced by the existence of NVIDIA's cuFFT library — a heavily optimized FFT implementation that one might reasonably assume extracts near-peak hardware performance.

The paper systematically demolishes this assumption by identifying two coupled, hardware-specific mismatches that no amount of butterfly optimization can fix. The first is that the butterfly's computational pattern — small complex multiplies with strided memory access — does not map to the tensor core matrix-multiply primitive that provides the GPU's highest throughput (a ~15× gap between tensor core and general-purpose FLOP rates on the H100). The second is that, independently of compute unit utilization, the butterfly's memory access pattern forces HBM spills for long sequences because the working set at each decomposition stage exceeds SRAM capacity, turning a compute-bound operation into a memory-bound one.

The key intellectual insight is that these two bottlenecks are coupled through the decomposition structure itself. Using larger matrix multiplies (lower decomposition order) improves tensor core utilization but increases SRAM pressure; using smaller matrix multiplies (higher order) fits in SRAM but underutilizes tensor cores. This coupling means that no single FFT implementation can be optimal across all sequence lengths — a finding that contrasts sharply with the standard approach of using a one-size-fits-all FFT library (cuFFT) for all input sizes.

This reframing is significant beyond the specific solution because it establishes a design methodology for hardware-efficient deep learning primitives: rather than optimizing an existing algorithm for the hardware, start from the hardware's ideal computational pattern (matrix multiply on tensor cores, SRAM-resident working sets) and restructure the algorithm to match that pattern. This is the same methodology that produced FlashAttention (restructuring attention to be SRAM-friendly), and FlashFFTConv demonstrates that it generalizes beyond attention to frequency-domain methods. The paper implicitly argues that the FFT convolution was never fundamentally slow — it was merely mapped to hardware through an algorithmic lens (Cooley-Tukey butterfly) rather than a hardware lens (tensor core matrix multiplies with SRAM locality).

Evidence for the diagnostic framing comes from the ablation study in Table 3, where the "Fusion-Only/cuFFTdx" baseline — which applies kernel fusion to the standard butterfly but does not use tensor cores — achieves meaningful speedup over PyTorch but (a) cannot use tensor cores and (b) fails entirely beyond 32K sequence length because fusion alone cannot keep the entire butterfly working set in SRAM. This confirms that both aspects of the reframing (tensor cores AND SRAM locality) are necessary, and neither alone suffices.

Innovation 2: The Order-p Cost Model as a Design Tool Connecting Mathematical Decomposition to Wall-Clock Time

While the Monarch decomposition provides a family of mathematically equivalent FFT factorizations (parameterized by the order $p$ ), the paper's critical practical contribution is a quantitative cost model (Equation 2) that predicts which $p$ will be fastest at a given sequence length on a given GPU. This transforms the Monarch decomposition from a mathematical curiosity into an engineering tool.

Prior to this work, the relationship between decomposition order and hardware efficiency was understood only qualitatively — higher order means smaller matrices (less compute, more I/O), lower order means larger matrices (more compute, less I/O). The paper's cost model makes this operational by introducing two specific thresholds that determine the optimal $p$ : the tensor core tile size $\mu = 16$ (below which matrix multiplies fall back to an order of magnitude slower general-purpose arithmetic) and the SRAM capacity (above which intermediate results spill to an order of magnitude slower HBM bandwidth). These thresholds are not theoretical — they are measured empirically for the target GPU (Appendix C, Table 19 for A100), using achievable rather than peak throughput numbers, which makes the model predictive rather than merely illustrative.

The conceptual significance is that the cost model explains the non-monotonic behavior visible in Figure 4: $p=4$ is actually worse than $p=3$ at short sequence lengths because the blocks are too small for tensor cores (the "Matrices Too Small for Tensor Cores" bump), while $p=2$ becomes worse than $p=3$ at longer lengths because its working set exceeds SRAM (the "SRAM Limit" threshold). This non-monotonicity means that a practitioner cannot simply choose "the highest $p$ that fits in SRAM" or "the lowest $p$ that uses tensor cores" — the optimal $p$ changes as sequence length crosses these hardware thresholds, and the cost model is what identifies the crossing points.

This is a fundamental contribution rather than an incremental refinement because it establishes a design principle that extends beyond FFT convolutions: for any algorithm that can be expressed as a family of decompositions with varying block sizes, there exist hardware-determined thresholds where the optimal decomposition changes, and a cost model that accounts for both compute throughput (as a piecewise function of block size relative to the tensor core tile size) and I/O bandwidth (as a piecewise function of working set size relative to SRAM capacity) can identify those thresholds. The paper does not explore this generalization, but the framework is directly applicable to other structured matrix decompositions (Butterfly matrices, Kaleidoscope matrices, low-rank factorizations) that present similar block-size trade-offs.

The evidence is Figure 4 itself, which maps the cost model predictions to sequence length and shows the predicted crossing points consistent with the empirical $p$ -selection used in Tables 3–4 ( $p=2$ for 256–32K, $p=3$ for 1M–2M, $p=4$ for 4M). The fact that these selections — derived from the cost model — produce monotonically improving speedup across four orders of magnitude validates the model's predictive power.

Innovation 3: The Discovery That Broadcasting Along Sequence Enables Fusion for Long Sequences

One of the paper's most subtle but consequential technical insights is the decision to broadcast matrix multiply operations along the sequence dimension rather than the batch/hidden dimensions (Section 3.1, Figure 3). This is not an obvious design choice — indeed, the natural parallelization strategy for FFTs (inherited from classical parallel FFT algorithms) is to distribute batch and channel items across compute units, since each item's FFT is independent and requires no communication between items. Broadcasting along the sequence, by contrast, means different SMs process different portions of the same sequence, requiring coordination that seems counterproductive.

The insight is that this inversion transforms the permutation operations in the Monarch decomposition from explicit data shuffles (which require HBM I/O) into in-SRAM matrix transposes (which are zero-cost). When the matrix multiply operates on an $N_1 \times N_1$ block that spans the sequence rather than the batch, the output layout is already organized such that a "permutation" is just reading the data in column-major rather than row-major order — a pure indexing operation with no data movement. This is what enables the inner layers of the decomposition to be fully fused (no HBM writes between matrix multiply and pointwise multiply operations), which in turn is what makes the approach viable for sequences up to 32K on current hardware.

Additionally, this broadcasting strategy reduces per-SM SRAM requirements, enabling fusion at longer sequence lengths than would be possible with batch-parallel execution. When parallelized over the batch, each SM must hold the entire $N_1 \times N_1$ working set for its assigned batch item; when parallelized over the sequence, each SM holds only a fraction of that working set. This is the mechanism that pushes the SRAM limit from ~2K to ~32K, a critical range for practical sequence modeling tasks.

This contribution is conceptually distinctive because it inverts the standard parallelization strategy for a counterintuitive reason: not for better parallel scalability (the batch-parallel approach already scales to arbitrary batch sizes), but for better memory locality within a single sequence. It is an instance of a broader principle: when the goal is kernel fusion (keeping intermediate data on-chip), the parallelization strategy should minimize the per-compute-unit working set for a single problem instance, even if that increases the number of compute units that must coordinate. This principle is not stated in the paper but is implicit in the design choice and represents a reusable insight for other sequence-level operations (e.g., parallel prefix scans, state-space model recurrences) where the natural parallelization axis may conflict with memory locality goals.

The evidence is indirect but compelling: the comparison between FlashFFTConv and the Fusion-Only/cuFFTdx baseline (Table 3) shows that cuFFTdx — which uses standard batch-parallel FFT fusion — cannot support sequences beyond 32K, while FlashFFTConv with its sequence-broadcast approach scales to 4M. The paper does not provide an ablation that isolates the broadcasting choice (e.g., a version of FlashFFTConv with batch-parallel rather than sequence-parallel execution), so the exact contribution of this choice versus the Monarch decomposition itself is not cleanly separated, but the mechanism is clearly explained in Section 3.1.

Innovation 4: Structural Sparsity as a First-Class Consequence of the Decomposition, Not an Afterthought

The paper introduces partial convolutions and frequency-sparse convolutions not as standalone architectural innovations but as natural consequences of the Monarch decomposition's compute structure — and this is what makes the contribution conceptually interesting. In most deep learning systems work, sparsity is an optimization applied after the algorithm is designed: train a dense model, prune weights, and then struggle to translate the reduced parameter count into actual speedup because unstructured sparsity maps poorly to GPU hardware. The typical outcome is that sparse models have fewer FLOPs on paper but are not faster in practice.

FlashFFTConv inverts this relationship: the sparsity pattern is designed to match the decomposition structure, and the decomposition structure is what makes the sparsity operational. When a portion of the frequency-domain kernel $k_f$ is zeroed along a dimension that corresponds to an iteration variable in the Monarch decomposition's nested loops (Appendix A.4), that iteration can be completely skipped — not just a few FLOPs saved, but an entire $N_1 \times N_1$ matrix multiply eliminated, along with its associated I/O. The sparsity pattern is co-designed with the compute decomposition so that structural zeros translate directly to structural compute savings.

This is a fundamental conceptual contribution rather than an incremental optimization because it demonstrates a design principle: when an algorithm is restructured to expose its computation as a series of independent blocks (the Monarch decomposition's block-diagonal matrix multiplies), sparsity in the coefficients becomes sparsity in the block schedule, which GPU hardware can exploit trivially (by skipping kernel launches or loop iterations). This principle is not specific to convolutions — it applies to any structured matrix operation where sparsity can be aligned with the block structure of the decomposition, such as sparse attention patterns aligned with block-sparse matrix multiply routines.

The evidence in Table 9 is striking: at 75% sparsity, the convolution is 1.3× faster and the model's perplexity actually improves (from 2.91 to 2.90). This is the opposite of the usual sparsity-quality trade-off where performance degrades as sparsity increases. The authors hypothesize that this improvement comes from "removing high-frequency noise," but regardless of the explanation, the result demonstrates that structural sparsity can be a quality-improving regularization rather than merely a compute-saving approximation — a finding with implications for how we think about convolution kernel design beyond the systems domain.

The connection to the partial convolution extension (Table 8) reinforces the same principle in the time domain: by designing the decomposition so that the effective kernel length controls the working set size, a model trained at one sequence length can be applied to much longer sequences with no quality degradation. This capability — extending HyenaDNA from 1M to 4M sequence length while maintaining perplexity — is a direct consequence of FlashFFTConv's decomposition-based memory management, not a separate algorithmic contribution. It demonstrates that systems-level design choices (how to structure the computation for memory locality) can enable qualitatively new modeling capabilities (embedding the longest human genes) that would be architecturally impossible under the standard implementation.

Innovation 5: The Efficiency-Quality Feedback Loop as a Systematic Empirical Finding

While the idea that "faster training enables better models under a fixed compute budget" is not novel in principle, the paper provides a systematic quantification of this effect for convolutional sequence models that has significant implications for how the field evaluates architectures. The finding is not merely that FlashFFTConv is faster — it is that the speedup translates to quality improvements equivalent to doubling the model's parameter count (Table 1: Hyena-s-155M with FlashFFTConv achieves perplexity 11.1, matching Hyena-m-355M with PyTorch; M2-BERT-base-110M with FlashFFTConv achieves GLUE 80.9, matching M2-BERT-large-260M with PyTorch).

This is a diagnostic contribution because it reframes how we should compare model architectures. The standard evaluation paradigm in sequence modeling is to train models with matched parameter counts and compare their quality. But this implicitly assumes that training throughput is irrelevant — that we care about quality per parameter rather than quality per dollar or quality per GPU-hour. The paper's result shows that this assumption produces misleading comparisons: Hyena-s (155M parameters) is often compared unfavorably to larger Transformers, but when the comparison is made at fixed compute budget rather than fixed parameter count, the smaller convolutional model with an efficient implementation can match the larger model's quality. This suggests that architectural comparisons should be normalized by total compute budget, not parameter count — a methodological point that extends beyond convolutions to any comparison between architectures with different training throughputs.

The paper does not argue this methodological point explicitly, but it is the clear implication of Table 1 and the reference results in Appendix B (Table 18). The mechanism is straightforward: FlashFFTConv achieves higher training throughput, so under a fixed wall-clock time budget, the model sees more training tokens (15B vs. 5B for Hyena-s) or more training steps (70,000 vs. 16,000 for M2-BERT-base). The quality improvement comes from more data, not a better architecture — but from a practitioner's perspective, the distinction is irrelevant because the compute budget is the binding constraint.

This finding is incremental in isolation (it is an expected consequence of faster training) but fundamental in context because the paper demonstrates it across two different model families (Hyena and M2-BERT), two different pretraining objectives (causal LM and masked LM), and two different downstream evaluation protocols (perplexity on the PILE and GLUE score). The consistency of the result — efficiency improvements translating to quality gains comparable to doubling parameters — suggests a robust empirical relationship, not a one-off fluke. Combined with the FLOPs-matched comparison against Transformers using FlashAttention-v2 (Table 6), which shows that Hyena with FlashFFTConv is faster than GPT with FlashAttention-v2 at sequence lengths 2K and above, the paper makes a compelling case that convolutional sequence models, when properly optimized, are not merely competitive with Transformers in quality but can be superior in throughput at the sequence lengths that matter for long-context applications.

The limitation, which the paper does not explore, is that this efficiency-quality feedback loop is bounded: it only works when the smaller model has sufficient capacity to absorb the additional training data. There is presumably a regime where the base model is too small to benefit from more tokens regardless of throughput, and the paper does not characterize where that boundary lies. But within the range of model sizes and datasets studied (155M–355M parameters, PILE-scale data), the effect is clear and practically significant.

5. Experimental Analysis

Evaluation Methodology

Dataset. The paper evaluates on multiple datasets spanning four modalities. For convolution benchmarking and end-to-end model throughput, synthetic or standardized inputs at specific sequence lengths are used (Tables 3–6). For quality experiments: language modeling uses the PILE for Hyena-GPT-style training and C4 for M2-BERT masked language modeling (Section 4.1, Appendix C.2); the Long Range Arena (LRA) benchmark provides Path-X (sequence length 16K) and Path-512 (sequence length 256K) for high-resolution image classification (Section 4.1); DNA modeling uses the HyenaDNA validation set from Nguyen et al. (2023) (Appendix C.9); and GLUE is used for M2-BERT fine-tuning evaluation (Table 1). The DNA gene embedding visualization uses the Ensembl genome dataset (Appendix B.3).
Base model(s). The paper uses several convolutional sequence models rather than a single unified architecture. For language modeling quality: Hyena-s-155M and Hyena-m-355M (Poli et al., 2023) for GPT-style causal LM, and M2-BERT-base-110M and M2-BERT-large-260M (Fu et al., 2023) for BERT-style masked LM (Table 1, Appendix C.2). For vision: a 6-layer convolutional model with hidden dimension 256 following the architecture from Fu et al. (2023) (Appendix C.3). For audio: SaShiMi (Goel et al., 2022) with 8 layers and hidden dimension 64 (Appendix C.5). For DNA: HyenaDNA-1M and HyenaDNA-450K pretrained checkpoints from Nguyen et al. (2023) (Tables 8–9). For the FLOPs-matched Transformer comparison: a 2.7B-parameter Hyena model and a parameter-matched GPT model using FlashAttention-v2 (Table 6, Appendix C.6). The diversity of architectures demonstrates that FlashFFTConv is a general-purpose convolution primitive, not tied to a specific model design.
Metrics. Four categories of metrics are used throughout the paper:
- Throughput: sequences per second (Table 5), tokens per second (Table 6), or images per second (Table 5), measured by timing the forward pass and scaling to the batch size. All timing measurements are averaged over 30 runs (Appendix C.4).
- Quality: perplexity (PPL) for language modeling and DNA modeling (lower is better), average GLUE score for M2-BERT (higher is better), and classification accuracy for Path-X/Path-512 (higher is better).
- FLOP utilization: computed as 2 × num_tokens × num_parameters for parametric FLOPs, plus the raw FLOP count from the cost model (Equation 2) for non-parametric convolution FLOPs, divided by the measured runtime and the GPU's theoretical peak throughput (Table 6, Appendix C.6).
- Memory footprint: measured as the additional memory allocated by the convolution operation relative to the baseline, reported in GB and as a reduction factor (Tables 3–4, Tables 16–17 in Appendix B).
Baselines. The paper compares against multiple baselines:
- PyTorch FFT convolution: the standard torch.fft implementation, which wraps cuFFT. This is the primary baseline for all convolution and gated convolution benchmarks (Tables 3–4).
- Fusion-Only / cuFFTdx: an ablation that uses kernel fusion (via NVIDIA's cuFFTdx library) but does not employ the Monarch decomposition for tensor cores, effectively representing the best achievable performance without restructuring the FFT into matrix multiplies (Table 3).
- FlashAttention-v2 (Dao, 2023): the highly optimized Transformer attention implementation, used as a baseline for end-to-end model throughput comparison at 2K, 8K, and 16K sequence lengths (Table 6).
- Larger model variants: Hyena-m-355M and M2-BERT-large-260M trained with PyTorch, serving as reference points for the "twice the parameters" quality comparison (Table 1, Appendix B Table 18).
- Majority voting baselines: for Path-X and Path-512, prior work achieving 50% (random) accuracy is the implicit baseline (Table 2).
Generation budget / compute accounting. For convolution micro-benchmarks, generation budget is not the relevant unit — instead, all convolutions are benchmarked at fixed sequence lengths with identical input/output dimensions (batch size 64, hidden dimension 768 for Tables 3–4, scaled to equivalent throughput for Table 5). For the fixed-compute-budget quality experiments (Table 1), the budget is wall-clock time: the FlashFFTConv model trains for 70,000 steps (M2-BERT) or 15B tokens (Hyena) in the same time that the PyTorch model trains for 16,000 steps or 5B tokens (Appendix C.2). For the FLOPs-matched comparison against Transformers (Table 6), FLOP utilization is computed directly from the measured runtime and model parameter count (Appendix C.6). For out-of-memory cases, the batch size or hidden dimension is split across multiple forward pass calls (Appendix C.4).
Cross-validation / statistical protocol. Most experiments are deterministic given fixed model seeds. Timing measurements use 30-run averaging (Appendix C.4). There is no explicit cross-validation for quality experiments — models are trained once with fixed hyperparameters. For the frequency-sparse convolution experiment (Table 9), a single pretrained HyenaDNA-1M checkpoint is sparsified post-hoc and evaluated on the standard validation set (Appendix C.9). The lack of multiple training runs or confidence intervals on quality metrics is a limitation — particularly for Table 1, where the claim of "matching models with twice the parameter count" is based on single training runs and could be sensitive to random seed variation.

Main Quantitative Results

Convolution Micro-Benchmarks: Speedup and Memory Savings

Headline result. FlashFFTConv accelerates exact FFT convolutions by up to 7.93× over PyTorch and reduces memory footprint by up to 8.21×. The speedup varies systematically with sequence length, decomposition order, and the presence of gating operations.

Standard convolution speedup (Table 3). Across sequence lengths from 256 to 4M, FlashFFTConv achieves speedups ranging from 1.33× at 4M (using p=4 decomposition) to 6.54× at 1K (using p=2). The speedup pattern is non-monotonic: it peaks at 1K (6.54×), remains above 4× through 8K, then gradually declines to 2.85× at 32K, and drops to 1.33×–1.82× at million-length sequences. This pattern directly reflects the hardware thresholds identified by the cost model (Figure 4): the peak regime corresponds to p=2 with blocks large enough for tensor cores (≥16×16) and working sets small enough for full SRAM fusion. The decline at longer sequences corresponds to (a) the transition to p=3 at ~1M where blocks shrink modestly, reducing tensor core utilization, and (b) the transition to p=4 at 4M where the outermost matrix multiply must spill intermediate results to HBM (Algorithm 4).

The Fusion-Only/cuFFTdx baseline (kernel fusion without Monarch decomposition) achieves substantial speedups of its own (2.85×–6.54× at shorter sequences) but plateaus at 2.85× at 32K and fails entirely beyond 32K — marked as "–" in Table 3. This confirms the paper's central claim that fusion alone is necessary but insufficient: without the Monarch decomposition's SRAM-friendly block structure, the working set for a fully fused butterfly FFT exceeds the 64 KB SRAM capacity beyond 32K, making fusion impossible. The Monarch decomposition's contribution is visible in the gap between FlashFFTConv and Fusion-Only at shorter lengths (e.g., 4.78× vs. a hypothetical fusion-only speedup implied by the 0.21 ms timing at 256) — this gap represents the tensor core utilization benefit.

Memory savings (Table 3). Memory footprint reduction ranges from 8.21× at 256 to 2.63× at 4M. The mechanism is recomputation in the backward pass (trading additional forward-pass FLOPs for not storing intermediate FFT results) plus kernel fusion (eliminating separate allocation for intermediate tensors). The sharp drop in memory savings at 65K (from 6.57× to 2.64×) corresponds to the transition where the outermost decomposition layers begin writing to HBM, requiring those intermediates to be stored rather than recomputed locally. The memory savings are measured as the additional memory from calling the convolution, excluding the input tensor footprint (Appendix B).

Gated convolution speedup (Table 4). Gated convolutions — the pattern y = v ⊙ ((u ⊙ w) * k) used in Hyena and M2-BERT — see greater speedups than standard convolutions: up to 7.93× at 1K, compared to 6.54× for the standard case. The additional speedup comes from fusing the elementwise gating operations into the same kernel, eliminating the HBM I/O that PyTorch incurs from three separate kernel launches (convolution, gate with w, gate with v). At 4M, the speedup is 1.30×, reflecting the dominance of the convolution's own I/O at extreme sequence lengths where the gating overhead is a smaller fraction of total time. Memory savings for gated convolutions are 6.65× at 256, declining to 2.81× at 4M — slightly better than the standard case at short lengths (6.65× vs. 8.21×) because the absolute memory footprint of a gated convolution is larger (more intermediates to save), but the relative savings are similar.

Causal convolution with implicit padding (Tables 13–14 in Appendix B). When the input length is half the FFT size (causal convolution with zero-padding), FlashFFTConv achieves speedups very close to the non-padded case across all sequence lengths. For example, at 1K, the standard convolution achieves 6.54× speedup (Table 3) while the padded causal convolution achieves 6.45× (Table 13). The small difference reflects the overhead of the bookkeeping for the real-valued FFT (decimation-in-time) relative to the savings from halving the FFT length — the net effect is approximately neutral. This confirms that the implicit padding optimization (skipping the outermost matrix multiply for zero-padded regions, Section 3.1) successfully eliminates the pad overhead that would otherwise make causal convolutions more expensive.

Backward pass speedup (Table 15 in Appendix B). The backward pass achieves speedups ranging from 6.43× at 512 to 1.28× at 4M. These are generally lower than the forward pass speedups (e.g., 6.43× vs. 6.61× at 1K for the forward pass from the more detailed Table 11 in Appendix B), reflecting the additional FLOPs from recomputation: the backward pass must re-execute most of the forward pass to recover the intermediates that were not stored, incurring compute overhead that partially offsets the memory bandwidth savings.

End-to-End Model Throughput

Headline result. FlashFFTConv accelerates convolutional sequence models end-to-end by 1.3×–4.4×, with the speedup varying by model architecture and the fraction of total runtime spent in the convolution.

Throughput across model architectures (Table 5). The five benchmarked models span four orders of magnitude in sequence length and four modalities:

M2-BERT-base (110M, seqlen 128): 1.9× throughput improvement (4,480 → 8,580 sequences/second). The convolution is relatively fast at this short sequence length, so the speedup is modest — the bottleneck shifts to other operations (attention-like mixing, MLP layers) after the convolution is optimized.
Hyena-s-4K (155M, seqlen 4K): 1.7× throughput improvement (84.1 → 147 sequences/second). The 4K sequence length is in the sweet spot where p=2 decomposition provides near-peak tensor core utilization with full SRAM fusion, but the model has substantial non-convolution compute (gating projections, MLP layers).
Long convs, Path-X (102M, seqlen 16K): 2.4× throughput improvement (126 → 308 images/second). The higher speedup reflects the convolution being a larger fraction of total runtime at 16K, where PyTorch's FFT implementation is increasingly bottlenecked by HBM I/O while FlashFFTConv remains SRAM-resident.
SaShiMi (5.4M, seqlen 64K): 1.3× throughput improvement (38.7 → 50.3 audio clips/second). The minimal speedup is explicitly attributed by the authors to this model interleaving convolutions with "SSM-based filter generation, pooling layers, and MLPs, which reduces the relative amount of time spent computing the convolution itself" (Section 4.2). This is a negative result in the sense that it reveals the limit of convolution optimization: when the convolution is not the dominant cost, optimizing it further yields diminishing end-to-end returns.
HyenaDNA-1M (seqlen 1M): 4.4× throughput improvement (0.69 → 3.03 sequences/second). The largest speedup occurs where PyTorch is most severely bottlenecked: at 1M sequence length, the PyTorch implementation can only fit batch size 1 on an 80GB GPU, whereas FlashFFTConv's memory savings enable batch size 4 (Section 4.2). The 4.4× speedup combines the per-iteration convolution speedup (~1.57× from Table 3 for standard convolution at 1M) with a 4× increase in batch size, yielding super-linear end-to-end throughput improvement.

Comparison Against FlashAttention-v2 Transformers

Headline result. A 2.7B-parameter Hyena model using FlashFFTConv achieves higher throughput than a parameter-matched GPT model using FlashAttention-v2 at sequence lengths 2K and above, despite lower FLOP utilization, because convolutions incur fewer total FLOPs.

Throughput and FLOP utilization (Table 6). At 2K sequence length, Hyena with FlashFFTConv processes 35.2K tokens/second versus 33.8K for GPT with FlashAttention-v2 — a 1.1× speedup. At 8K, the speedup grows to 1.3× (35.2K vs. 27.8K), and at 16K, to 1.5× (32.3K vs. 21.6K). The GPT model's throughput drops with increasing sequence length (as expected from attention's quadratic scaling in practice, even with FlashAttention's optimizations), while Hyena's throughput remains nearly flat (35.2K → 35.2K → 32.3K), consistent with convolutions' O(N log N) scaling.

The FLOP utilization numbers tell a more nuanced story. FlashAttention-v2 achieves 65.7%–78.5% end-to-end FLOP utilization, while FlashFFTConv achieves 56.5%–62.3%. FlashFFTConv is 10 percentage points less efficient at extracting the GPU's theoretical peak throughput, yet still wins on wall-clock time because the convolution operation simply requires fewer total FLOPs than attention at these sequence lengths. This is a critical finding: peak hardware utilization is not the right metric when comparing operations with different asymptotic complexities. A less efficient implementation of a cheaper operation can be faster than a highly efficient implementation of an expensive operation.

The authors note that FlashFFTConv achieves "only 10% less than FlashAttention-v2" in FLOP utilization (Section 4.2) — framing the 62.3% vs. 72.1% gap at 8K as relatively small. This is fair in the context of systems optimization (62% utilization is very high for a complex FFT-based operator), but the comparison is somewhat apples-to-oranges: FlashAttention-v2's utilization is measured on matrix multiplies that naturally achieve high tensor core occupancy, while FlashFFTConv's utilization includes the FFT decomposition which has inherently lower arithmetic intensity in its permutation and twiddle correction steps.

What this comparison demonstrates and what it doesn't. Table 6 demonstrates that convolutional models can be as fast as or faster than Transformers at sequence lengths where convolution's asymptotic advantage overcomes attention's implementation maturity advantage. However, the comparison is at a single model scale (2.7B parameters) and three specific sequence lengths. The authors do not explore whether this advantage holds at smaller scales (where the constant-factor overhead of FlashFFTConv's multi-stage decomposition might dominate) or much larger scales (where both models become compute-bound and FLOP utilization becomes the binding constraint). Additionally, the GPT baseline uses FlashAttention-v2 but the Hyena baseline uses the author's own implementation — the comparison would be strengthened by also benchmarking Hyena with the standard PyTorch FFT convolution to isolate the FlashFFTConv contribution from the architectural difference (convolution vs. attention).

Quality Improvements Under Fixed Compute Budget

Headline result. Given the same wall-clock training budget, FlashFFTConv enables models to see more training data, translating to quality improvements equivalent to doubling model parameters: Hyena-s achieves 2.3 points better perplexity (13.4 → 11.1) and M2-BERT-base achieves 3.3 points higher average GLUE score (77.6 → 80.9).

Language modeling quality (Table 1). The mechanism is straightforward and explicitly quantified: the PyTorch Hyena-s model trains on 5B tokens in the available time, while the FlashFFTConv Hyena-s trains on 15B tokens — triple the data (Appendix C.2). The resulting perplexity of 11.1 matches the PyTorch Hyena-m (355M parameters, 2.3× larger) trained on 5B tokens, which achieves 11.1 (Appendix B, Table 18). Similarly, M2-BERT-base trains for 70,000 steps with FlashFFTConv versus 16,000 with PyTorch — over 4× more updates — and achieves 80.9 GLUE score, matching the PyTorch M2-BERT-large (260M parameters, 2.4× larger) at 81.0.

These comparisons do not demonstrate that convolutions are superior to Transformers, or that smaller models with more training data universally match larger models. They demonstrate a more specific and limited claim: under this particular fixed compute budget and these particular model scales, the FlashFFTConv implementation enables enough additional training to close the gap to the next model size. The result would not hold if the compute budget were larger (eventually both implementations saturate on available data), if the base model were too small to absorb additional tokens, or if the larger model also used FlashFFTConv (which would increase its training throughput as well, restoring the relative gap).

The absence of confidence intervals is notable: both numbers in Table 1 are single training runs, and the reference larger-model numbers in Table 18 are also single runs. Given known variance in large model training (sensitivity to random seed, data ordering), it is possible that the 11.1 vs. 11.1 match for Hyena is fortuitous and would not replicate exactly. The qualitative claim — "efficiency improvements translate to quality gains" — is robust even if the exact numerical match is not.

Longer Sequence Models: Solving Path-512

Headline result. FlashFFTConv enables the first model to achieve above-random accuracy on Path-512 (256K sequence length), reaching 96.1% accuracy, where all prior implementations fail with out-of-memory errors or achieve only 50% (random guessing).

Path-X and Path-512 results (Table 2). On Path-X (16K sequence length), both PyTorch and FlashFFTConv achieve 96.9% accuracy — the sequence is short enough that PyTorch's convolution implementation is functional and the model capacity is sufficient. On Path-512 (256K, a 16× increase in sequence length over Path-X), PyTorch cannot run at all ("✗" in Table 2, indicating out of memory), while FlashFFTConv achieves 96.1%.

The 96.1% accuracy on Path-512 is close to the 96.9% on Path-X, suggesting that the task does not become fundamentally harder at higher resolution — the same convolutional architecture, given the ability to process the longer sequence, solves it. This implies that the previous failure to solve Path-512 was purely a systems limitation (insufficient memory and no support for sequences of this length), not an algorithmic limitation of convolutional models. This is a strong result for the paper's motivating narrative: systems advances can unlock qualitatively new capabilities that were architecturally out of reach.

However, the Path-512 model uses a different architecture than Path-X: 4 layers vs. 6 layers, and kernel dropout 0.1 vs. 0.3 (Appendix C.3). The filter length is capped at 65,536 rather than being allowed to grow to the full 256K sequence length, which effectively makes this a partial convolution. The paper does not report an ablation showing that the full 256K filter is necessary or that the capped filter is the reason for the 0.8-point accuracy drop from Path-X. This is a minor gap — the headline result (first model to solve Path-512) stands regardless.

Partial Convolutions: Reduced Memory and Length Extension

Headline result. Partial convolutions — zeroing out later portions of the kernel — reduce memory footprint during training by up to 5.6× (32.5G → 5.8G from full 8K to effective 256) while maintaining perplexity, and enable extending a pretrained 1M-sequence model to 4M sequence length with no quality degradation.

Memory reduction during training (Table 7). A Hyena-s-8K model trained with progressively shorter effective kernels shows that the kernel can be pruned from 8K to 2K (a 4× reduction in filter length) with no perplexity degradation (13.8 at all three lengths). Further pruning to 1K causes a small degradation (13.9), and to 512 or 256 causes modest degradation (14.0, 14.2). The memory footprint drops proportionally: 32.5 GB at full 8K to 11.8 GB at 2K, and to 5.8 GB at 256.

Length extension of pretrained models (Table 8). A HyenaDNA model pretrained at 1M sequence length, when applied to 4M inputs using a sliding-window partial convolution approach, achieves perplexity 2.90 — slightly better than its 1M validation perplexity of 2.91. A HyenaDNA-450K model extended to 4M achieves 2.91. The authors describe this as yielding "the first model that can embed the longest human genes at single nucleotide resolution (2.3M base pairs)" (Section 4.3), which is supported by the perplexity numbers and the t-SNE visualization in Appendix B (Figure 5) showing the dystrophin gene (the longest human gene) clustered with other protein-coding genes.

The finding that perplexity does not degrade — and may slightly improve — when extending to longer sequences is non-obvious. It suggests that the HyenaDNA model's learned filters do not depend strongly on the absolute position beyond ~1M, or that the biological signal in DNA is sufficiently local that a 1M effective context window captures most of the relevant dependencies. The paper does not investigate this mechanistic explanation, but the empirical result is practically significant: it means a model trained at substantial computational expense on 1M sequences can be deployed on 4M sequences without retraining.

Frequency-Sparse Convolutions: Speedup Without Quality Loss

Headline result. Frequency-sparse convolutions — zeroing out portions of the frequency-domain kernel k_f — can achieve up to 1.8× convolution speedup, with model quality maintained or slightly improved up to 75% sparsity (1.3× speedup) and modest degradation at 91% sparsity.

Quality-speedup trade-off (Table 9). Starting from a pretrained HyenaDNA-1M model (PPL 2.91), structured frequency-domain sparsity is applied following the pattern in Appendix A.4 (Table 10). At 50% sparsity, perplexity is unchanged at 2.91 with 1.2× convolution speedup. At 75% sparsity, perplexity improves to 2.90 (the best in the table) with 1.3× speedup. At 79% sparsity, perplexity returns to 2.91 with 1.4× speedup. Beyond this, quality degrades: 84% sparsity → PPL 2.93, 91% sparsity → PPL 2.98.

The speedup numbers are measured on the convolution operation itself, not end-to-end model throughput. They scale non-linearly with sparsity: from 1.0× (dense) to 1.2× (50%) to 1.3× (75%) to 1.4× (79%) to 1.5× (84%) to 1.8× (91%). The non-linearity arises because the first dimensions to be zeroed (the outermost decomposition layers, controlled by parameters a and b in Table 10) allow skipping inner matrix multiply blocks, while later sparsity (parameters c and d) skips entire loop iterations in the outer decomposition layers, which provides larger per-element savings (Appendix A.4).

The quality improvement at 75% sparsity is the most intriguing result. The authors hypothesize this is "potentially as a result of removing high-frequency noise" (Section 4.3). If this hypothesis is correct, it suggests that the HyenaDNA model's pretrained kernels contain high-frequency components that are not useful (or are mildly harmful) for the DNA modeling task, and that frequency-sparse convolution acts as a form of post-hoc regularization. This has implications beyond systems optimization: it suggests that convolution kernel parameterizations may benefit from explicit frequency-domain regularization during training, not just post-hoc sparsification.

A limitation is that the sparsity is applied to a single pretrained model. The paper does not explore training from scratch with frequency-sparse kernels, which could potentially achieve the speedup without any quality loss (or with improved quality if the sparsity acts as a beneficial inductive bias). The post-hoc nature of the experiment means the speedup is "free" in the sense that it costs no additional training, but a model trained from scratch with the sparsity constraint might find different (potentially better) weight configurations.

Ablation Studies and Robustness Checks

Fusion-Only / cuFFTdx ablation (Table 3): Removing the Monarch decomposition (using only kernel fusion on the standard butterfly FFT) achieves 2.85×–6.54× speedup at sequences up to 32K but fails entirely beyond 32K due to SRAM capacity limits. This isolates the Monarch decomposition's contribution as enabling (a) tensor core utilization at shorter lengths and (b) SRAM-friendly working sets at longer lengths. The gap between FlashFFTConv and Fusion-Only at 256 (4.78× vs. an effective ~2× implied by the 0.21 ms timing) quantifies the tensor core contribution, while the ability to run at 1M–4M (where Fusion-Only cannot) quantifies the memory hierarchy contribution.

Domain-specific optimization ablation (Tables 3–4 vs. Tables 13–14): Comparing standard convolutions (Tables 3–4) against causally-padded convolutions (Tables 13–14 in Appendix B) shows nearly identical speedups, confirming that the implicit padding optimization (skipping outermost matrix multiplies for zero-padded regions) effectively eliminates the pad overhead. At 1K, the standard convolution achieves 6.54× speedup (Table 3) while the padded version achieves 6.45× (Table 13) — a negligible 1.4% difference. Similarly, gated convolution speedup is 7.93× (Table 4) vs. 6.75× for padded (Table 14) at 1K, with the difference attributable to the pad's elimination of some outer matrix multiplies slightly changing the FLOP-to-I/O ratio.

Backward pass recomputation ablation (Table 3 vs. Table 15): The backward pass achieves lower speedups than the forward pass (e.g., 4.37× vs. 6.61× at 1K) because recomputation adds forward-pass FLOPs to the backward pass. The memory savings from recomputation (Tables 16–17) are substantial — 6.40×–8.21× at short sequences — but the paper does not provide an ablation showing what memory usage would be without recomputation (i.e., storing all intermediates). This makes it difficult to isolate the recomputation contribution from the fusion contribution to memory savings.

Sequence length sweep (Tables 11–17 in Appendix B): Full results for all sequence lengths in powers of two from 256 to 4M are provided for standard convolutions, gated convolutions, causally-padded convolutions, causally-padded gated convolutions, backward pass, and memory usage. These comprehensive tables reveal that the speedup and memory savings trends are smooth and monotonic within each decomposition regime (p=2, p=3, p=4) but exhibit discontinuities at the decomposition boundaries where the cost model predicts transitions. For example, speedup for standard convolution drops from 2.84× at 32K (p=2) to 2.08× at 65K (p=2, but at the edge of SRAM capacity) to 1.57× at 1M (p=3) — the inflection points align with the cost model thresholds in Figure 4.

FP16 vs. BF16 precision: The paper mentions that FlashFFTConv supports both fp16 and bf16 (Appendix A.2), but all benchmarks are reported in fp16 (implied by the use of tensor cores, which operate on 16-bit floats). There is no ablation comparing fp16 vs. bf16 quality or throughput. For DNA modeling and language modeling tasks, bf16 is often preferred for training stability due to its larger dynamic range — the lack of bf16 benchmarks is a minor gap.

Monarch decomposition order sweep: The paper implicitly sweeps p through the cost model and the sequence length tables (p=2 for 256–32K, p=3 for 1M–2M, p=4 for 4M), but does not provide an explicit ablation showing, e.g., p=2 vs. p=3 at 1M sequence length. Such a comparison would directly validate the cost model's prediction that p=3 outperforms p=2 at this length. The implicit evidence is that the authors chose p=3 for 1M (Table 3) and that this choice produces a 1.57× speedup; a p=2 implementation at 1M would likely be slower or out of memory, but this is not shown.

Critical Assessment

Claim 1 (Abstract, Section 4.2): "FlashFFTConv speeds up exact FFT convolutions by up to 7.93× over PyTorch." This claim is well-supported by the gated convolution benchmark at sequence length 1K (Table 4). However, the "up to" qualifier is doing substantial work: the 7.93× speedup occurs at a single sweet spot (1K, gated, H100, batch size 64, hidden dimension 768) and speedups at other sequence lengths range from 1.30× to 6.54×. The claim is accurate for the specific configuration tested but should not be interpreted as a typical or expected speedup. The paper does not explore how speedup varies with batch size or hidden dimension, which could significantly affect the absolute numbers (smaller dimensions mean smaller matrix multiplies, potentially dropping below the tensor core threshold and reducing speedup).

Claim 2 (Abstract, Section 4.1): "Given the same compute budget, FlashFFTConv allows Hyena-GPT-s to achieve 2.3 points better perplexity... matching models with twice the parameter count." This claim is supported with important caveats. The 2.3-point perplexity improvement (13.4 → 11.1) is clearly measured and the mechanism (3× more training tokens) is well-documented. The claim of "matching models with twice the parameter count" relies on comparing against a PyTorch-trained Hyena-m-355M (PPL 11.1, Table 18 in Appendix B) — a single training run with no confidence interval. The match is exact (11.1 vs. 11.1), which is suspiciously precise and likely reflects fortuitous rounding rather than a fundamental equivalence. More importantly, the comparison is not symmetric: the larger model would also benefit from FlashFFTConv if implemented, and would presumably see its own throughput improvement (though likely smaller in relative terms, since larger models spend more time in non-convolution operations). The claim is best understood as "FlashFFTConv enables enough additional training on the small model to close the gap to the next model size when the larger model uses a slower implementation" — a conditional claim that the paper's language does not always make explicit.

Claim 3 (Abstract): "FlashFFTConv achieves 96.1% accuracy on Path-512, a high-resolution vision task where no model had previously achieved better than 50%." This claim is strongly supported. The 96.1% accuracy is a dramatic improvement over the 50% random-guessing baseline, and the paper is transparent that prior implementations failed due to out-of-memory errors (not poor architecture design). However, the claim that "no model had previously achieved better than 50%" is historical context, not a controlled experimental comparison — it is possible that an equivalently optimized Transformer (using FlashAttention-v2 with a sparse attention pattern) could also solve Path-512, but this experiment was not run. The paper's contribution is demonstrating that convolutional models can solve this task with sufficient systems support, not that they are uniquely capable of doing so.

Claim 4 (Abstract, Section 4.3): "Partial convolutions yield the first DNA model that can process the longest human genes (2.3M base pairs)." This claim is supported by the 4M sequence length extension results (Table 8) and the t-SNE visualization (Figure 5 in Appendix B). The ability to process a 2.3M-base-pair gene at single-nucleotide resolution follows directly from supporting 4M sequence length. However, "process" is a weak verb — the paper demonstrates that the model can embed these genes (produce a vector representation), not that it can perform specific downstream tasks (variant effect prediction, regulatory element identification) at this resolution. The quality metric is perplexity (which improves slightly), but perplexity on held-out DNA sequence is a proxy for biological usefulness, not a direct measure. The claim is accurate as stated but the practical significance depends on whether single-nucleotide-resolution embeddings of full-length genes enable downstream applications that coarser models cannot — a question the paper does not answer.

Claim 5 (Section 4.2, Table 6): "FlashFFTConv is faster in wall-clock time than FlashAttention-v2 end-to-end at sequence lengths 2K and longer." This claim is supported for the specific models and sequence lengths tested but has limited generalizability. The comparison is between a 2.7B Hyena and a 2.7B GPT, both at 2K, 8K, and 16K sequence lengths, on A100 GPUs. The speedup (1.1×, 1.3×, 1.5×) is modest and highly dependent on the model architecture: Hyena makes different accuracy-efficiency trade-offs than GPT, and the comparison conflates the architectural difference with the systems optimization. A fairer systems-only comparison would benchmark the same convolution operation implemented via FlashFFTConv vs. the best available alternative (cuFFTdx or PyTorch) and separately compare convolution-based models against attention-based models. The paper does the former (Tables 3–4) but presents the latter (Table 6) as a systems comparison when it is actually an architecture comparison mediated by systems.

Missing experiments that would strengthen the paper:

Batch size and hidden dimension sweeps. All convolution benchmarks use batch size 64, hidden dimension 768. The cost model (Equation 2) shows that cost scales linearly with BH, but the tensor core utilization threshold (N_i >= 16) depends on N_i, which is a function only of sequence length and decomposition order, not batch or hidden size. Small hidden dimensions (e.g., 64 or 128) might produce matrix multiplies too small for tensor core efficiency even at moderate sequence lengths, reducing speedup. The paper does not characterize this regime.
End-to-end training wall-clock time for the quality experiments. Table 1 reports quality improvements "given the same compute budget" but does not specify what that budget was in GPU-hours. Without this number, a practitioner cannot assess whether the additional training tokens (15B vs. 5B) represent a realistic budget constraint. If the fixed budget was, say, 24 hours on 8×A100, then FlashFFTConv's throughput advantage is sufficient to triple the tokens seen; if the budget was 1 hour, the advantage might be proportionally similar but the absolute quality would be lower.
Ablation of the Monarch decomposition broadcasting strategy. The paper claims that broadcasting along the sequence dimension (rather than batch/hidden) is critical for enabling fusion at longer sequences (Section 3.1, Figure 3), but provides no ablation comparing sequence-parallel vs. batch-parallel FlashFFTConv. The cuFFTdx baseline is batch-parallel but also lacks tensor core support, so it confounds two differences. An ablation using Monarch decomposition with batch-parallel broadcasting would isolate the contribution of this design choice.
Quality experiments with frequency-sparse convolutions trained from scratch. Table 9 applies sparsity post-hoc to a pretrained model; training from scratch with the sparsity constraint could yield different (potentially better) quality-speedup trade-offs, since the model could adapt its non-zero frequency components to the sparsity pattern during training.
Multiple training runs for the fixed-budget quality comparison. The central quality claim (Table 1) rests on single training runs. Given known variance in large model training, reporting mean and standard deviation over 3–5 seeds would substantially increase confidence that the 2.3 perplexity point improvement is reliable.
Scaling the compute budget. The fixed-budget experiment (Table 1) uses a specific budget where FlashFFTConv's throughput advantage translates to 3× more tokens. At larger budgets (more total GPU-hours), both implementations would see more data, and the relative advantage might shrink if the model saturates. At smaller budgets, the advantage might grow. The paper does not explore how the quality gap varies with total compute.

Overall assessment. The experiments convincingly demonstrate that FlashFFTConv substantially accelerates FFT convolutions (2–8× depending on configuration) and that this acceleration translates to meaningful downstream benefits: higher training throughput, longer workable sequence lengths, and reduced memory footprint. The results are comprehensive in their coverage of sequence lengths (four orders of magnitude), model architectures (five model families), and modalities (four). The primary weakness is that the quality claims — particularly the "twice the parameters" comparison — are based on single training runs with a specific compute budget, and the generalizability of these claims to other budgets, model scales, and random seeds is not established. The systems benchmarks (speedup, memory, FLOP utilization) are more robust because they measure deterministic properties of the convolution implementation; the quality benchmarks combine systems effects with training dynamics and are subject to more sources of variance that the paper does not control for.

6. Limitations and Trade-offs

6.1 The Cost Model Assumes Tuned, Sequence-Specific Kernels — But the Tuning Overhead Is Not Accounted For

The assumption or constraint. FlashFFTConv's core contribution — the order-p Monarch decomposition with p selected via a hardware cost model — requires separate CUDA kernel implementations for each sequence length and decomposition order. The paper is explicit about this in Appendix A.2: "To ensure high performance, we implement CUDA kernels for each specific sequence length, allowing us to cater to specific performance nuances that arise from the decomposition at that sequence length." The cost model (Equation 2, Section 3.2) treats p as a free variable to be optimized per sequence length, but this optimization produces a design-time decision, not a runtime adaptation. Each (sequence length, decomposition order) pair requires a hand-tuned kernel with customized block dimensions, tile sizes, loop unrolling factors, and register allocation.

The consequence. The headline speedup numbers (Tables 3–4) measure the performance of kernels that have already been tuned for their target sequence length. They do not account for the engineering cost of producing those kernels. For a practitioner deploying convolutions at a new sequence length not in the paper's benchmark set — say, 12,288 tokens for a specific application — there are two unappealing options: (1) pad to the nearest supported length (e.g., 16,384), wasting computation, or (2) develop and tune a new kernel from scratch, requiring deep CUDA expertise and substantial engineering time. The paper's cost model tells you which p to use, but it does not automate the kernel generation — unlike libraries such as cuFFT which handle arbitrary transform sizes automatically (albeit without the tensor core and fusion optimizations FlashFFTConv provides).

This limitation is particularly acute because modern sequence models increasingly use non-power-of-two sequence lengths or dynamically varying sequence lengths (e.g., packing variable-length documents, processing audio of arbitrary duration). FlashFFTConv's tight coupling between sequence length and kernel implementation makes it brittle in these settings. The paper acknowledges that it "look[s] forward to integrating more general libraries such as Cutlass to support a wider range of GPUs, and developing support for non-GPU accelerators" (Appendix A.5), but this is framed as future work, not a current capability.

What evidence exists in the paper. The comprehensive benchmarks (Tables 11–17 in Appendix B) report results only for powers of two from 256 to 4M — no intermediate sequence lengths appear anywhere in the paper. The cost model plot (Figure 4) suggests that the model could predict performance at intermediate lengths, but no such validation is provided. The architecture-specific optimizations described in Appendix A.2 — "aggressively tune our kernel hyperparameters such as block and tile dimensions, and loop unrolling factors for the best performance on the specific underlying hardware" — confirm that each kernel is manually optimized for its target configuration.

Mitigation status. The paper does not attempt to mitigate this limitation. It suggests future integration with Cutlass (a template library for GPU matrix multiplies) as a direction for broader hardware support, but Cutlass would address the GPU-generation portability issue (V100 vs. A100 vs. H100), not the sequence-length-specificity issue. The fundamental challenge — that optimal tile sizes depend on the specific matrix dimensions in the Monarch decomposition, which in turn depend on sequence length and p — is inherent to the approach and would require an auto-tuning framework (analogous to how cuBLAS auto-tunes matrix multiply kernels) to resolve. The paper does not propose such a framework.

6.2 Hard Problems Remain Essentially Unsolved — Test-Time Compute Amplifies Existing Capability but Cannot Create It

The assumption or constraint. FlashFFTConv optimizes the computation of the FFT convolution but does not change the convolution's expressiveness or what the model can represent. If a model architecture lacks the capacity to solve a task, no amount of convolution efficiency will help — the convolution computes exactly the same mathematical function, just faster and with lower memory. This limitation is parallel to the finding in the analyzed compute-optimal test-time scaling paper: test-time compute amplifies existing capability but cannot create it from nothing.

The consequence. For tasks that require capabilities beyond what the base convolutional model can represent — e.g., reasoning that requires modelling dependencies spanning hundreds of thousands of tokens at full resolution, or tasks where the convolutional inductive bias is fundamentally mismatched — FlashFFTConv provides no benefit beyond enabling the model to fail faster. The paper demonstrates this implicitly: on Path-512 (Table 2), the same convolutional architecture that achieves 96.9% on Path-X achieves 96.1% on Path-512 when FlashFFTConv enables the longer sequence, but this is because the task structure is preserved at higher resolution — the model already had the capability and was only blocked by systems constraints.

More subtly, the partial convolution results (Table 7) reveal that Hyena-s-8K's convolution kernel can be pruned from 8K to 2K (a 4× reduction) with no perplexity degradation. This suggests that the effective context length actually used by the model is substantially shorter than the nominal sequence length — the model does not benefit from the full 8K context even when it is provided. FlashFFTConv enables processing longer sequences, but if the architecture does not learn to use the additional context, the longer sequences provide no quality benefit — only additional computational cost.

What evidence exists in the paper. The paper's quality experiments (Section 4.1) show improvements from faster training (more tokens seen) rather than from longer sequence modeling per se. The Path-512 result (Table 2) is the only experiment where longer sequence length directly improves task performance, and it is a carefully controlled vision task where "longer sequence = higher image resolution" directly translates to more information. For language modeling (Tables 1, 7), the benefits come from training throughput, not from the model leveraging longer contexts. The frequency-sparse convolution results (Table 9) further suggest that the pretrained convolution kernels contain high-frequency components that can be zeroed out without quality loss — implying that the models are not effectively utilizing their full representational capacity.

Mitigation status. The paper does not address this limitation — it is outside the scope of a systems paper. However, the implication is important for practitioners: FlashFFTConv is a force multiplier for convolutional architectures, not a substitute for architectural innovation. If the goal is to solve tasks that require qualitatively different reasoning capabilities (rather than the same reasoning applied to longer inputs), FlashFFTConv alone will not suffice. The paper's positioning — "we hope that understanding how to optimize the FFT convolution can also inspire algorithmic innovation" (Section 1) — implicitly acknowledges this boundary.

6.3 The Speedup Numbers Exclude the Memory Overhead of the Cost Model's Difficulty Estimation and the Engineering Cost of Kernel Development

The assumption or constraint. The headline 7.93× speedup (Table 4) measures the convolution operation in isolation — timing a single forward pass of a convolution kernel given pre-loaded inputs, pre-computed FFT matrices, and pre-loaded twiddle factors. It does not account for the one-time costs of (1) computing the FFT of the convolution filter k_f (though this is amortized across batch items and sequence positions, as the paper notes in Section 2.1), (2) loading the Monarch decomposition matrices F, F^{-1}, and twiddle factors t, t_inv from host to device memory (Algorithm 1 shows these are loaded from HBM, implying they are already on-device), or (3) the memory footprint of storing these decomposition-specific matrices for multiple sequence lengths if a model processes variable-length inputs.

More importantly, the paper does not quantify the engineering cost of producing FlashFFTConv itself. The approach requires: implementing the Monarch decomposition in custom CUDA kernels, tuning those kernels for each (sequence length, decomposition order) pair on each target GPU architecture, implementing the real-valued FFT bookkeeping for decimation-in-time, fusing gating operations, and maintaining separate code paths for forward and backward passes with recomputation logic. This is a substantial software engineering investment that the paper's authors undertook but that a typical practitioner deploying a convolutional sequence model cannot replicate without specialized systems expertise.

The consequence. The "up to 7.93×" speedup is an idealized measurement that a practitioner will approach but not fully realize in a production deployment, for several reasons. First, the one-time costs of loading decomposition matrices and pre-computing filter FFTs reduce the effective speedup for models that process few tokens per convolution (e.g., inference with batch size 1). Second, the memory overhead of storing decomposition matrices for multiple sequence lengths increases the GPU memory footprint beyond what the convolution memory savings alone would suggest — this overhead is not measured. Third, the engineering cost means that FlashFFTConv is not a drop-in replacement for torch.fft — adopting it requires either using the paper's open-source code (which may not support the exact sequence lengths or GPU architecture needed) or reimplementing the approach, which is a multi-person-month effort requiring deep CUDA expertise.

For practitioners deciding whether to adopt FlashFFTConv, the relevant question is not "what speedup does the paper report?" but "what speedup will my model see in my deployment setting given my engineering resources?" The paper provides extensive benchmarking data to answer the first part of this question (the speedup varies from 1.30× to 7.93× depending on configuration) but no data to answer the second part.

What evidence exists in the paper. The end-to-end throughput numbers (Table 5) partially address this by measuring full model throughput rather than isolated convolution speed. These show speedups of 1.3×–4.4×, which are substantially lower than the convolution micro-benchmark speedups (1.33×–7.93×). The gap represents the fraction of total model runtime spent in the convolution — for SaShiMi, only 1.3× end-to-end speedup despite likely substantial convolution speedup, because the model interleaves convolutions with other operations (SSM filter generation, pooling, MLPs) that FlashFFTConv does not accelerate. The paper is transparent about this: "Speedup varies by the size of the models and the relative amount of time spent computing the convolution compared to other parts of the models" (Section 4.2). However, the paper does not provide a breakdown of where runtime is spent for each model, making it impossible for a practitioner to estimate their own expected speedup from the published numbers alone.

Mitigation status. The paper acknowledges the end-to-end speedup variability and provides per-model measurements (Table 5), which is appropriate for a systems paper. The engineering cost is not addressed — the paper does not estimate development time, provide a usability comparison against cuFFT, or discuss integration complexity. The open-source release (implied by the repository structure but not explicitly discussed) would partially mitigate the engineering cost for users whose needs align with the supported configurations, but the paper does not characterize which sequence lengths, GPU architectures, or precision formats the released code supports.

6.4 Single Hardware Generation and GPU-Only Scope — No Path to CPUs, Edge Devices, or Non-NVIDIA Accelerators

The assumption or constraint. FlashFFTConv is developed for and evaluated exclusively on NVIDIA A100 and H100 GPUs. The paper explicitly states: "FlashFFTConv was developed on A100 GPUs, and tested on A100 and H100 GPUs. Older generations of GPU such as V100 are not supported, since the sizes of the tensor cores are different" (Appendix A.5). The tensor core tile size μ = 16 (for the 16 × 16 matrix multiply unit) is hardcoded into the cost model (Equation 2) and the decomposition strategy — a GPU with a different tensor core size (e.g., V100's 4 × 4 or future architectures with larger tiles) would require re-deriving the optimal decomposition and re-tuning all kernels.

The consequence. FlashFFTConv's approach is fundamentally tied to a specific hardware primitive (NVIDIA tensor cores of a specific generation) and a specific memory hierarchy (HBM → SRAM → registers with the capacities and bandwidths of the A100/H100). This has several practical consequences:

No CPU deployment. Convolutional models that might be deployed for CPU inference (e.g., on-device DNA analysis, embedded audio processing) cannot benefit from FlashFFTConv at all, since CPUs lack tensor cores and have an entirely different memory hierarchy. The paper does not discuss CPU implementations even as future work.
No support for older or commodity GPUs. The V100, which remains widely deployed in cloud instances and academic clusters, is explicitly unsupported. The consumer GPU line (RTX series) has smaller tensor cores and less SRAM per SM, which would shift the optimal decomposition thresholds and potentially make the approach less effective or non-functional.
Vendor lock-in to NVIDIA. The approach relies on CUDA-specific primitives (WMMA API for tensor cores, shared memory for SRAM fusion, warp-level primitives). Porting to AMD GPUs (which have an analogous matrix core unit but a different programming model), Intel GPUs, or dedicated AI accelerators (TPUs, Cerebras, SambaNova) would require a near-complete rewrite. The paper mentions "non-GPU accelerators" as future work (Appendix A.5) but provides no abstraction layer or hardware-agnostic formulation that would facilitate porting.
Fragility to hardware evolution. The H100's tensor core supports FP8 in addition to FP16/BF16, and future architectures will likely have different tensor core sizes, different SRAM capacities, and different bandwidth ratios. Each generation shift potentially invalidates the tuned kernel parameters and the cost model's hardware constants (which the paper measures empirically, Appendix C, Table 19, rather than deriving from architecture specifications).

What evidence exists in the paper. All experiments are on A100-40GB or H100-SXM GPUs (Appendix C.1). The empirical GPU constants in Table 19 are measured specifically for the A100-40GB. The paper acknowledges the V100 incompatibility explicitly (Appendix A.5) and frames broader hardware support as future work. There are no experiments on different GPU generations, no discussion of CPU fallback paths, and no abstraction that would insulate users from hardware-specific details.

Mitigation status. The paper does not attempt to mitigate this limitation beyond mentioning future work (Appendix A.5). The cost model (Equation 2) is parameterized by hardware constants (τ_M, τ_G, σ_H, σ_S, μ), which in principle could be re-measured for different hardware, but the kernel implementations themselves are hand-tuned CUDA with hardcoded tile sizes and tensor core usage patterns. A more portable approach — e.g., expressing the Monarch decomposition in a high-level tensor compiler (TVM, Triton) that auto-generates hardware-specific code — is not explored. For practitioners committed to the NVIDIA A100/H100 ecosystem (which includes most cloud GPU offerings as of the paper's writing), this limitation is not immediately blocking, but it precludes deployment on the broader landscape of ML hardware.

6.5 The Real-Valued FFT Optimization Assumes a Specific Input Layout and Does Not Generalize to Complex-Valued or Multi-Channel Convolutions

The assumption or constraint. FlashFFTConv's decimation-in-time optimization (Section 3.1, Appendix A.1) exploits the fact that the convolution inputs and outputs are real-valued — the even-indexed samples go into the real part of a half-length complex vector, and the odd-indexed samples go into the imaginary part. This halves the FFT length from N to N/2, providing approximately a 2× reduction in the core FFT cost.

However, this optimization only works for real-valued inputs and outputs. Many convolution variants in sequence modeling — particularly those used in state-space models (S4, S4D, Mamba) and some convolutional architectures — operate on complex-valued hidden states, or interleave real and complex operations. Additionally, convolutions with complex-valued kernels (e.g., to represent phase information) would not benefit from this optimization. Even for real-valued convolutions, the optimization imposes a specific input layout: the convolution must operate along the sequence dimension with real values, which is true for the benchmarked architectures (Hyena, M2-BERT, the LRA convolutional model) but may not hold for architectures that apply convolutions along the channel dimension or that mix real and complex representations.

The consequence. A practitioner implementing a model that uses complex-valued convolutions, or convolutions along non-sequence dimensions, cannot use the decimation-in-time optimization and would need to fall back to a full-length complex FFT. This would approximately double the FFT cost relative to the optimized real-valued case, reducing the headline speedup numbers by a factor of up to 2× for the FFT portion of the computation (though the pointwise multiply and kernel loading costs would be unchanged, so the end-to-end convolution speedup would degrade by less than 2×).

More subtly, the paper does not provide an ablation measuring the contribution of the real-valued FFT optimization to the overall speedup. The implicit padding optimization (skipping outermost matrix multiplies for zero-padded causal convolutions) is similarly architecture-specific: it works for causal convolutions where half the input is zeros, but circular convolutions (used in some bidirectional models) do not benefit. The paper benchmarks both causal and non-causal convolutions (Tables 3 vs. Table 13), showing nearly identical speedups — but this comparison conflates the padding savings with the real-valued FFT savings, making it impossible to isolate either contribution.

What evidence exists in the paper. The decimation-in-time procedure is described mathematically in Appendix A.1, and its implementation is mentioned in Section 3.1 ("Domain-Specific Optimizations"). The speedup numbers in Tables 3–4, 11–14 all include this optimization — there is no ablation comparing a full-length complex FFT implementation against the optimized half-length version. The paper's statement that the real-valued FFT "cut[s] the FFT cost in half" (Section 3.1) is a theoretical claim, not an empirically validated speedup contribution. The interaction between this optimization and the Monarch decomposition order is also unexplored: for small N_i (high p), the bookkeeping overhead of the real-to-complex conversion might become non-negligible relative to the matrix multiply cost, reducing the effective savings.

Mitigation status. The paper does not address this limitation. For the models benchmarked — Hyena, M2-BERT, the LRA convolutional model — the real-valued assumption holds, so the limitation does not affect the reported results. However, the paper's framing (Section 3.1 presents the real-valued FFT as a "domain-specific optimization" alongside implicit padding and gating fusion) does not explicitly state that this optimization restricts the space of supported convolutions. A practitioner reading the paper might assume FlashFFTConv is a general-purpose FFT convolution accelerator, when in fact its best-case performance numbers assume real-valued, causally-padded, gated convolutions with length equal to a power of two — a specific, though common, configuration.

6.6 The Memory Savings Claim Includes Recomputing the Forward Pass, but the Latency Cost of This Recomputation Is Not Separately Reported

The assumption or constraint. FlashFFTConv uses recomputation in the backward pass to reduce memory footprint: "Instead of storing intermediate results on HBM for the backward pass (e.g., the intermediate result of F u), we simply recompute them in the backward pass" (Section 3.1). This means that during training, the forward pass is effectively executed twice per training step — once to compute the output (and the loss), and once during the backward pass to regenerate the intermediate activations needed for gradient computation.

The memory savings are clearly reported: up to 8.21× reduction in convolution memory footprint (Table 3, Table 16). However, the latency overhead of this recomputation — the additional wall-clock time spent in the backward pass recomputing forward-pass values — is not cleanly separated from the speedup numbers.

The consequence. The backward pass speedups reported in Table 15 (Appendix B) — e.g., 4.37× at 1K vs. 6.61× forward pass speedup at the same length (Table 11) — reflect both the Monarch decomposition's hardware efficiency and the additional FLOPs from recomputation. A practitioner cannot determine from the published numbers how much of the backward pass slowdown relative to the forward pass is due to recomputation versus inherent differences in forward/backward FLOP counts for the convolution operation. This matters because:

Training throughput depends on both forward and backward pass time. If the backward pass is disproportionately slower due to recomputation, the effective training speedup may be lower than the forward-pass speedup suggests. The paper reports end-to-end training throughput improvements (Table 5: 1.3×–4.4×), which account for this, but the forward-pass-only speedup numbers (7.93×) are prominently featured and may mislead readers who do not check the backward pass performance.
The recomputation cost scales with sequence length differently than the forward pass savings. For p=4 at 4M sequence length, the forward pass speedup is 1.33× but the backward pass speedup is 1.28× (Table 15) — the gap between forward and backward performance narrows at extreme lengths, suggesting that recomputation overhead becomes a smaller fraction of total time when HBM I/O dominates. This interaction is not analyzed.
Inference-only deployments do not benefit from recomputation at all. The memory savings from recomputation are irrelevant for inference, where there is no backward pass. The paper does not report memory footprint for inference-only mode (storing all intermediates), making it difficult to assess FlashFFTConv's advantage for inference workloads where memory capacity, not training throughput, is the binding constraint.

What evidence exists in the paper. The backward pass speedups in Table 15 (Appendix B) are reported alongside forward pass speedups (Tables 11–14), but the paper's main text (Section 4.2) focuses on forward pass numbers. The text mentions recomputation as a memory-saving strategy (Section 3.1) but does not quantify its impact on backward pass time or on the forward/backward speedup gap. The end-to-end model throughput numbers (Table 5) implicitly include both forward and backward pass effects for training, but these are model-specific and not decomposed into convolution-forward vs. convolution-backward contributions.

Mitigation status. The paper does not provide an ablation comparing training throughput with and without recomputation, nor does it report inference-only memory usage or latency. The recomputation strategy is presented as an unqualified benefit ("reduces memory footprint and I/O cost") without acknowledging that it trades increased FLOPs for reduced memory — a classic time-memory trade-off whose net effect on training throughput depends on whether the model is compute-bound or memory-bound at the relevant scale. For the benchmarked configurations, the trade-off is evidently favorable (since end-to-end speedups are positive), but the paper does not characterize the regimes where it would be unfavorable. This is a missed opportunity to provide practitioners with guidance on when to enable or disable recomputation.

7. Implications and Future Directions

How This Work Changes the Landscape

FlashFFTConv does not propose a new convolution algorithm, a new model architecture, or a new theoretical result about sequence modeling. It is a systems optimization paper — and yet its implications extend well beyond the systems community into how the field designs, evaluates, and deploys convolutional sequence models. The shift is best understood on three levels: diagnostic, methodological, and strategic.

Diagnostic: Reframing the FFT convolution as a hardware-mapping problem rather than an algorithm problem. Before FlashFFTConv, the dominant narrative around convolutional sequence models was that their $O(N \log N)$ asymptotic complexity made them inherently more efficient than Transformers at long sequence lengths, and that their wall-clock inferiority was a temporary implementation gap that would close with standard engineering. The paper demonstrates that this narrative was exactly backwards: the FFT convolution is not inherently efficient on modern hardware — its butterfly compute pattern fundamentally mismatches tensor core matrix-multiply units and its memory access pattern forces expensive HBM round-trips for sequences longer than SRAM capacity. These are not minor implementation details that better cuFFT tuning could fix; they are structural properties of the standard FFT algorithm when mapped to the GPU memory hierarchy.

By reformulating the problem as "how do we restructure the computation to match what the hardware is good at?" rather than "how do we optimize the existing algorithm for the hardware?", the paper performs the same diagnostic move that FlashAttention performed for attention mechanisms. The shift in perspective is what makes the solution possible: once you accept that the butterfly network is the wrong computational primitive for modern GPUs, you can look for decompositions (like the Monarch factorization) that express the same mathematical transform using the right primitives (dense matrix multiplies on tensor cores, SRAM-resident working sets). The cost model (Equation 2, Section 3.2) formalizes this diagnostic by quantifying precisely where the standard approach fails — the tensor core size threshold μ = 16 and the SRAM capacity threshold — and using these to select the decomposition order p per sequence length.

This diagnostic reframing resolves a tension that has existed in the sequence modeling literature since at least 2021. Several papers claimed that convolutional or state-space models were asymptotically more efficient than Transformers (citing $O(N \log N)$ or $O(N)$ complexity), while other papers showed that Transformers with FlashAttention were faster in practice at all but the most extreme sequence lengths. The apparent contradiction came from comparing theoretical FLOP counts against wall-clock measurements of different implementations. FlashFFTConv provides the missing piece: the theoretical advantage of convolutions can be realized in wall-clock time, but only when the implementation is agnostic to the butterfly and restructured for tensor cores and SRAM locality. The paper proves this with the FlashAttention-v2 comparison (Table 6): at 2K sequence length and 2.7B parameters, the convolution model achieves 1.1× higher throughput than the Transformer model despite 10 percentage points lower FLOP utilization — the convolution's lower total FLOPs finally win when the implementation is competitive.

Methodological: Establishing that hardware-aware decomposition is a design tool, not an implementation detail. The paper's most reusable conceptual contribution is the order-p cost model as a bridge between mathematical structure and hardware performance. The Monarch decomposition itself existed before FlashFFTConv (Dao et al., 2022), but was studied as a way to parameterize efficient neural network layers — its inference-time hardware implications were unexplored. FlashFFTConv demonstrates that the decomposition order p is not a fixed hyperparameter but rather a hardware-dependent tuning knob whose optimal value changes systematically with sequence length, crossing distinct thresholds defined by tensor core tile size and SRAM capacity.

This is a methodological contribution because it provides a template for how to think about other structured matrix operations. The pattern is: (1) identify a family of mathematically equivalent decompositions parameterized by block size, (2) model the compute cost as a piecewise function of block size relative to the tensor core tile size (below which throughput drops ~13×), (3) model the I/O cost as a piecewise function of working set size relative to SRAM capacity (above which bandwidth drops ~14×), and (4) select the decomposition order that minimizes the sum. This template applies directly to other structured matrices used in efficient ML — Butterfly matrices, Kaleidoscope matrices, low-rank factorizations with block structure, structured sparse attention patterns — and the paper's empirical validation of the cost model (Figure 4 predicting the transitions visible in Tables 3, 11–15) gives confidence that the template is practical, not just theoretical.

The sequence-broadcasting design choice (Section 3.1, Figure 3) is a second methodological contribution of independent interest. The standard parallelization strategy for sequence operations is to distribute batch items across compute units, since they are independent. FlashFFTConv inverts this — broadcasting the matrix multiply across the sequence dimension — for the counterintuitive reason that it reduces per-SM memory pressure, enabling kernel fusion at longer sequence lengths. This principle — that parallelization strategy should be chosen to minimize per-compute-unit working set for a single problem instance rather than to maximize independent parallelism — is not widely appreciated in the ML systems literature, where data parallelism across the batch dimension is the default. It is likely applicable to other sequence-level operations (parallel scans, recurrences, state-space model updates) where the natural parallelization axis conflicts with memory locality goals, and FlashFFTConv provides a concrete case study of the benefit (supporting sequences to 32K with full fusion vs. ~2K with batch-parallel fusion).

Strategic: Shifting the efficiency-quality trade-off from architecture design to systems implementation. The fixed-compute-budget experiment (Table 1) demonstrates something qualitatively important: the choice of convolution implementation can be as consequential for downstream model quality as doubling the parameter count. This is not because FlashFFTConv changes what the convolution computes — it computes the exact same mathematical function — but because under a realistic compute constraint, a faster implementation enables more training data in the same wall-clock time, and more training data improves quality. The Hyena-s model with FlashFFTConv achieves perplexity 11.1, matching a Hyena-m model (355M parameters, ~2.3× larger) trained with PyTorch — a gain that comes entirely from training on 15B tokens instead of 5B.

This finding has strategic implications for how the field should evaluate architectures. The standard practice is to compare models at matched parameter counts or matched FLOP counts, implicitly assuming that training throughput is irrelevant to quality. FlashFFTConv's result suggests this practice systematically undervalues architectures that can achieve high hardware utilization, because the quality benefit of that utilization — more training tokens per dollar — is not captured in parameter-matched comparisons. If we compare Hyena-s (FlashFFTConv) against a Transformer at fixed parameter count, Hyena-s might look worse; if we compare at fixed wall-clock training budget, the gap narrows or reverses.

This is not a new insight — the importance of training efficiency has been argued since at least the Chinchilla scaling laws (Hoffmann et al., 2022), which normalized by FLOPs rather than parameters — but FlashFFTConv provides a concrete, quantified case study specifically for convolutional sequence models, and at a scale (155M–355M parameters, PILE-scale data) that is representative of typical academic and small-industry training budgets. The implication is that architecture research and systems research cannot be separated: a new convolution variant that is 10% more parameter-efficient but 2× slower on hardware is, under a fixed budget, a worse model than the original — even if it looks better in a parameter-matched table. FlashFFTConv makes this argument empirically rather than rhetorically.

What becomes more attractive as a result. (1) Research into hardware-efficient structured matrix decompositions — FlashFFTConv validates the Monarch decomposition as a practical tool, opening the door for Butterfly, Kaleidoscope, and other structured matrices to be evaluated not just for parameter efficiency but for inference-time hardware utilization. (2) Co-design of sparsity patterns and compute schedules — the frequency-sparse convolution results (Table 9) show that structural sparsity designed to match the decomposition's block structure yields actual speedup, not just parameter-count reduction. This suggests a research program where sparsity patterns are chosen as much for their hardware mapping as for their statistical properties, and where sparsity is evaluated by wall-clock time, not FLOP counts. (3) Convolutional architectures for long-context applications that were previously out of reach — Path-512 (256K, Table 2) and 4M-length DNA modeling (Table 8) are now tractable, and the paper provides evidence that the quality ceiling (96.1% on Path-512, 2.90 perplexity on 4M DNA) is high enough to motivate further architectural work on these tasks.

What becomes less attractive. (1) Naive FFT convolutions as a baseline — the paper establishes that the PyTorch/cuFFT implementation is sufficiently suboptimal (up to 7.93× slower!) that comparing against it is no longer informative about the potential of convolutional architectures. Future architecture papers should benchmark with FlashFFTConv (or equivalent optimized implementations) when reporting wall-clock time, or clearly separate the systems contribution from the architectural contribution. (2) The "FFT convolution is inherently efficient" narrative — the paper's detailed bottleneck analysis (Section 3) makes it impossible to claim asymptotic complexity implies wall-clock efficiency on modern accelerators, and future work must grapple with the tensor core / memory hierarchy mismatch explicitly. (3) Unstructured sparsity as a recipe for speedup — the frequency-sparse convolution results (Table 9) demonstrate that sparsity only translates to speedup when it aligns with the computation's block structure, and that unstructured sparsity (randomly zeroing individual frequency bins) would not achieve the 1.4×–1.8× speedups reported. This reinforces the broader lesson from the sparse neural network literature: sparsity must be co-designed with the compute schedule to be operational.

Follow-Up Research This Work Enables

Auto-tuning the Monarch decomposition for arbitrary sequence lengths and hardware targets. FlashFFTConv demonstrates the cost model for selecting p at powers-of-two sequence lengths on A100/H100, but the kernels are hand-tuned for each configuration. A natural next step is to build a code generator that takes a sequence length N, a GPU architecture description (tensor core tile size μ, SRAM capacity, achievable bandwidths), and a precision (fp16/bf16/fp8), and produces an optimized FlashFFTConv kernel automatically. This would address the brittleness limitation — the current implementation supports only specific sequence lengths on specific GPUs because each kernel is manually tuned. The cost model (Equation 2) provides the objective function for selecting p and the factorization N = ∏ N_i, but the low-level CUDA parameters (tile sizes, loop unrolling factors, register allocation) require empirical search — a task well-suited for auto-tuning frameworks like TVM or Triton. A strong follow-up would auto-generate kernels for all powers-of-two from 256 to 4M on A100, H100, and at least one consumer GPU (RTX 4090), and show that the auto-tuned kernels achieve ≥90% of the hand-tuned performance across the board. The key research question is whether the cost model's piecewise-linear structure (compute-bound below μ, memory-bound above SRAM) is sufficiently predictive that auto-tuning can focus on the regime boundaries rather than exhaustively searching the entire space.

Combining the Monarch FFT decomposition with state-space model recurrences (S4, S4D, Mamba). The paper benchmarks only explicit convolution models (Hyena, M2-BERT, the LRA convolutional model), but the most prominent line of long-context sequence models — state-space models (S4, S4D, Mamba, S5) — also rely on FFT convolutions for their parallel training mode. During training, these models expand their recurrent state-space representation into an explicit convolution (the "HiPPO" convolution kernel) and compute it via FFT; during inference, they switch to a recurrent mode that is more efficient per-token. FlashFFTConv should directly accelerate the training-time FFT convolution in these models, but two complications arise: (1) the SSM convolution kernels are often complex-valued (since they represent continuous-time dynamics), so the real-valued FFT optimization would not apply, and (2) the SSM dimension (the number of independent state-space channels) is typically smaller than the hidden dimension in explicit convolution models, which could shift the tensor core utilization thresholds. A concrete experiment: benchmark S4 and Mamba training throughput with FlashFFTConv vs. the current PyTorch/cuFFT implementation at sequence lengths 2K–32K, measuring the impact on total training time for a fixed-step budget, and report whether the speedup matches the 1.7×–2.4× range observed for Hyena and the LRA model (Table 5). This would test whether the approach generalizes beyond explicit convolutions to the broader class of FFT-based sequence mixing primitives, and would identify whether the real-valued assumption can be relaxed without losing the efficiency gains.

Frequency-sparse convolutions trained from scratch with the sparsity pattern as a learned parameter. The paper demonstrates that post-hoc frequency sparsity on a pretrained HyenaDNA model improves or maintains perplexity up to 75% sparsity (Table 9). This is promising but leaves open the question: if sparsity is beneficial as a regularizer, can we train with sparsity from scratch and achieve better quality than the dense baseline, not just equal quality with speedup? A strong follow-up would parameterize the sparsity pattern during training — e.g., by learning a mask over the frequency-domain kernel k_f with a continuous relaxation (like a sigmoid temperature parameter) that is annealed toward discrete sparsity — and measure whether the resulting model achieves lower perplexity than the dense baseline at matched effective FLOPs. The key hypothesis is that frequency-sparse training acts as a spectral regularizer, forcing the model to rely on low-frequency signal components that generalize better, analogous to how weight decay or dropout regularize in the time domain. The cost model (Equation 2) and the structured sparsity patterns (Appendix A.4) provide the framework for making the sparsity operational: any learned sparsity pattern can be mapped to skipped blocks in the Monarch decomposition, translating reduced parameter count into actual speedup during training as well as inference. A negative result — that learned sparsity patterns do not improve over post-hoc sparsification — would still be informative, suggesting that the benefit is purely denoising of pretrained weights rather than a better inductive bias.

Extending the cost model to multi-GPU and distributed settings for sequences beyond 4M. FlashFFTConv's current ceiling is 4M sequence length (p=4 on H100), beyond which the outermost matrix multiply working sets exceed HBM capacity and the approach fails. However, scientific applications — whole-genome DNA modeling (human genome: ~3 billion base pairs), long-duration audio (hours at 64 kHz), full-document understanding — require or would benefit from sequences of 10M–1B tokens. A natural extension is to distribute the sequence across multiple GPUs, using the Monarch decomposition's block structure to minimize communication: since each N_i × N_i block multiply is independent, the decomposition naturally exposes sequence parallelism where different GPUs own different portions of the sequence and communicate only at permutation boundaries. A concrete proposal: implement a model-parallel FlashFFTConv where the outermost decomposition layer is split across GPUs along the sequence dimension (exploiting the block-diagonal structure), with the cost model extended to include inter-GPU communication bandwidth (NVLink, ~300 GB/s on A100) as an additional term in ω(i). The research question is whether the decomposition order p and the degree of model parallelism can be jointly optimized — adding more GPUs reduces per-GPU memory pressure (allowing lower p and better tensor core utilization) but introduces communication cost at the permutation steps. A strong follow-up would characterize this trade-off and demonstrate convolutional models operating on 32M–128M sequence lengths across 8×A100 nodes.

Characterizing the boundary where efficiency gains stop translating to quality gains. The fixed-compute-budget experiment (Table 1) shows that FlashFFTConv's 3× higher training throughput translates to 2.3 PPL improvement for Hyena-s on 15B tokens vs. 5B. But this experiment uses a specific budget, model size, and dataset — it does not characterize the shape of the throughput-to-quality curve. There must exist regimes where additional training tokens no longer improve quality (model capacity saturation), and regimes where FlashFFTConv's speedup is insufficient to close the gap to the next model size. A systematic study would sweep: model sizes from 50M to 1B parameters, total training token budgets from 1B to 100B, and throughput multipliers from 1× to 4× (simulating PyTorch vs. FlashFFTConv), measuring downstream quality (PPL or GLUE) at each point. The goal is to produce an efficiency-quality scaling law — analogous to the Chinchilla scaling laws for pretraining — that predicts, for a given model size and compute budget, whether investing in systems optimization (FlashFFTConv) or in model scaling (more parameters, PyTorch) yields better quality per dollar. This would convert the paper's point estimate ("matches models with twice the parameter count") into a predictive framework that practitioners can use to make build-vs-buy decisions about convolution implementations.

Stress-testing the real-valued FFT optimization on complex-valued and multi-channel convolution variants. The decimation-in-time optimization (Appendix A.1) provides up to 2× savings for real-valued convolutions, but state-space models, complex-valued convolutions, and architectures interleaving real and complex operations cannot use it. A useful stress test would implement FlashFFTConv with and without the real-valued optimization, benchmarking both on: (1) standard real-valued Hyena/M2 convolutions (where the optimization should help), (2) complex-valued SSM convolutions from S4/S4D (where it cannot apply), and (3) a synthetic benchmark with varying ratios of real-to-complex channels. This would quantify the real-valued optimization's contribution to the headline speedup numbers, and would provide guidance for architecture designers — if the optimization provides a 1.5–2× training throughput advantage, that may justify designing models to use real-valued convolutions even when complex-valued parameterizations have other advantages. The interaction with the Monarch decomposition order is also worth characterizing: for high p (small N_i), the relative overhead of the real-to-complex bookkeeping grows, potentially eating into the savings. A comprehensive ablation would measure speedup vs. p with and without the real-valued optimization across sequence lengths 1K–1M.

Practical Applications and Downstream Use Cases

Long-read genomic analysis and variant calling at single-nucleotide resolution. The paper demonstrates that FlashFFTConv enables HyenaDNA to process 4M-length DNA sequences — sufficient to embed the longest human genes including dystrophin (2.3M base pairs) — with perplexity 2.90, essentially identical to the 1M model's perplexity of 2.91 (Table 8). This has direct application in clinical genomics: whole-genome sequencing produces reads of 10K–100K base pairs with technologies like PacBio and Oxford Nanopore, and analyzing these reads at single-nucleotide resolution without tokenization or downsampling could improve variant calling accuracy (detecting single-nucleotide polymorphisms, insertions, and deletions). The practical workflow would use a pretrained HyenaDNA-4M model (extended from the 1M checkpoint via partial convolutions) as an embedder, with a lightweight classifier on top for variant effect prediction. The key advantage over existing approaches (DNABERT, Enformer, Nucleotide Transformer) is resolution: these models tokenize DNA into k-mers (typically k=3–6), losing single-nucleotide precision, whereas FlashFFTConv's memory efficiency makes full-resolution modeling tractable at gene-scale lengths. The speedup numbers (4.4× end-to-end for HyenaDNA-1M, Table 5) mean that inference on a full human genome (~3 billion base pairs, ~30,000 genes) becomes practical on a single GPU node — without FlashFFTConv, the PyTorch implementation achieves only 0.69 sequences/second at batch size 1 (Table 5), making genome-scale inference infeasible.

High-resolution satellite and medical image classification. The Path-512 result (96.1% accuracy at 256K sequence length, Table 2) demonstrates that convolutional models with FlashFFTConv can process images at resolutions that are simply out of memory for standard implementations. This generalizes beyond the synthetic Path-X/Path-512 tasks to real-world high-resolution imaging: satellite imagery (where identifying small objects requires processing gigapixel images at full resolution), digital pathology (where whole-slide images are routinely 100K × 100K pixels), and astronomical surveys (where detecting transient events requires scanning terabyte-scale sky maps). The practical architecture would flatten 2D images into 1D sequences (following the Path-X formulation) and apply a convolutional classifier, with the key advantage being that FlashFFTConv's $O(N \log N)$ scaling makes processing these resolutions tractable without the $O(N^2)$ cost that would make attention-based approaches prohibitive. The paper's memory reduction (up to 8.21×, Table 3) is particularly relevant here: a 256K × 256K image flattened to 65B tokens would require hundreds of GB of memory with standard FFT convolution, but FlashFFTConv's recomputation and fusion strategies reduce the per-sequence memory footprint proportionally, potentially bringing it within reach of a single GPU. The Path-512 result (96.1% at 256K) provides evidence that the quality does not degrade at high resolution — the model architecture can use the extra information.

Real-time audio processing and generation on raw waveforms. The SaShiMi benchmark (Table 5) achieves a modest 1.3× end-to-end speedup because the model interleaves convolutions with other operations, but for architectures where the convolution dominates runtime — such as WaveNet-style autoregressive models or real-time speech enhancement models operating on raw audio at 44.1–64 kHz — the convolution micro-benchmark speedups (2–6× for sequences of 16K–64K, Table 3) would translate more directly. A practical deployment scenario is on-device real-time audio denoising or source separation on a smartphone or hearing aid: the raw audio is processed in 1-second chunks (64K samples at 64 kHz), and FlashFFTConv's memory savings (up to 7.73× at 1K–4K, Table 16) mean the convolution fits in the limited memory of an edge GPU. The end-to-end model speedup depends on the architecture, but the gated convolution variant — which reflects the y = v ⊙ ((u ⊙ w) * k) pattern common in audio models — achieves 3.22×–7.93× speedup in the 16K–1K sequence length range (Table 4), suggesting that real-time processing at low latency becomes feasible where it was not before. The key practical advantage over existing FFT-based audio processing (which uses cuFFT from C libraries) is the tensor core utilization — the same audio model running FlashFFTConv consumes less power and achieves lower latency on the same hardware, which is critical for battery-powered devices.

Cost-efficient training of convolutional language models at academic compute scales. A typical academic lab might have access to 8×A100 GPUs for a fixed wall-clock duration — say, a 1-week allocation on a cluster. Given this constraint, the fixed-compute-budget result (Table 1) has a direct practical interpretation: using FlashFFTConv for Hyena-s training yields a model with perplexity 11.1 (matching a 355M-parameter model), while using PyTorch yields perplexity 13.4 — a difference of 2.3 PPL points that could determine whether a paper's results are state-of-the-art or merely incremental. The practical recommendation is that any lab training convolutional language models at the 100M–1B parameter scale should adopt FlashFFTConv (or an equivalent optimized implementation) as the default, because the throughput improvement (1.7× for Hyena-s-4K, 1.9× for M2-BERT-base, Table 5) translates directly to better models under realistic compute constraints. This is not a speculative future benefit — the paper provides the reference PyTorch numbers and the FlashFFTConv numbers for the exact architectures and datasets used in the literature (Hyena, M2-BERT), so labs can directly compute their expected quality improvement given their specific budget. The key practical caveat is the sequence-length specificity: labs must use sequence lengths supported by the released kernels (powers of two from 256 to 4M), but this covers the standard lengths used in the literature (128, 512, 1K, 2K, 4K, 8K).