FFTMatvec: Split FFT for Toeplitz Matvec

• FFTMatvec is a split Fast Fourier Transform algorithm that computes matrix-vector products for block Toeplitz matrices without full circulant embedding.
• It interleaves lazy embedding with eager projection to reduce redundant data manipulation, thereby cutting operation counts and peak memory usage.
• The method offers improved performance and parallel scalability in electromagnetic and acoustic simulations, capitalizing on structured matrix properties.

FFTMatvec denotes a split Fast Fourier Transform algorithm for matrix-vector multiplication with block Toeplitz structure, introduced for settings in which multilevel circulant embedding is effective but overly expansive in memory and work. In numeric modeling of electromagnetics and acoustics, where block Toeplitz operators arise naturally and may not be highly sparse, the standard FFT route embeds each Toeplitz level into a circulant one, performs FFTs and diagonal multiplication, and thereby accelerates the product relative to naive multiplication. FFTMatvec modifies that paradigm through a “lazy embedding, eager projection” strategy that interleaves embedding and projection rather than materializing the full circulant extension, with asymptotic savings in both operation count and peak memory, and with a recursive branch structure that exposes straightforward parallelism (Siron et al., 2024).

1. Mathematical setting and target operator

The algorithm is formulated for a $d$-dimensional vector $v\in\mathbb C^{n^d}$ reshaped on an $n\times n\times\cdots\times n$ grid. The associated matrix $T$ is block Toeplitz in the sense that its entries depend only on relative grid offsets along each dimension; equivalently, $T$ is a $d$-level Toeplitz operator with generating “bands” $t[i_1,\ldots,i_d]$, each of size $n^d$. The computational objective is to evaluate $y=T\,v$ without explicitly forming $T$, since explicit storage and naive multiplication both scale as $v\in\mathbb C^{n^d}$0 (Siron et al., 2024).

This problem class is especially relevant when the Toeplitz structure originates from translationally invariant kernels sampled on regular grids. The motivating example given for the non-sparse regime is the electromagnetic Green function in a spatial basis. In such cases, FFT-based acceleration is already standard, but the standard route incurs a systematic overhead because the embedding step introduces coefficients that are not needed for the final projected output.

2. Classical circulant embedding and its overheads

The conventional multilevel method doubles the length in each of the $v\in\mathbb C^{n^d}$1 dimensions, embedding the Toeplitz structure into a circulant structure. After all $v\in\mathbb C^{n^d}$2 embeddings, the vector length grows from $v\in\mathbb C^{n^d}$3 to $v\in\mathbb C^{n^d}$4. The computation then proceeds by FFTs along each dimension, diagonal multiplication by the embedded spectrum, and inverse FFTs. In the notation of the paper, with $v\in\mathbb C^{n^d}$5, the total cost is approximated by

$v\in\mathbb C^{n^d}$6

and the peak memory requirement for vectors and data is

$v\in\mathbb C^{n^d}$7

These formulas make explicit that the acceleration is obtained by paying for a full $v\in\mathbb C^{n^d}$8 expansion of the problem representation (Siron et al., 2024).

The practical consequence is that the FFT-based method is often much faster than naive multiplication but not information-theoretically economical. The full embedded representation contains padding and coefficients that are subsequently discarded by projection back to the physical domain. FFTMatvec is designed precisely to remove that avoidable surplus while retaining the FFT-based computational structure.

3. Split-FFT construction: lazy embedding and eager projection

The central device is an even-odd frequency splitting induced by a phase shift. In one dimension, for $v\in\mathbb C^{n^d}$9, define the zero-padded vector $n\times n\times\cdots\times n$0 and the diagonal phase matrix

$n\times n\times\cdots\times n$1

Then $n\times n\times\cdots\times n$2 separates into two interleaved spectra: the even-index spectrum is $n\times n\times\cdots\times n$3, and the odd-index spectrum is $n\times n\times\cdots\times n$4. Conversely, if $n\times n\times\cdots\times n$5 and $n\times n\times\cdots\times n$6 are the two length-$n\times n\times\cdots\times n$7 spectra, the inverse reconstruction is

$n\times n\times\cdots\times n$8

This identity is the basis for splitting an embedded transform into branches that never require the full padded vector to be stored (Siron et al., 2024).

The multidimensional algorithm is described recursively by a routine $n\times n\times\cdots\times n$9. At level $T$0, the routine applies an FFT along dimension $T$1. If $T$2, it creates a phase-shifted child branch, recurses on the original and shifted branches, and then merges them. At the leaves, it performs diagonal multiplication with the corresponding local Toeplitz data block $T$3. If $T$4, it ends by applying an inverse FFT along dimension $T$5. In the pseudocode nomenclature of the paper, the branch split is denoted $T$6, merging is $T$7, and the leaf action is $T$8. The key structural fact is that no full $T$9 padding is ever materialized; each branch operates on the original size $T$0 vector (Siron et al., 2024).

This mechanism can be viewed as replacing one monolithic embedded FFT by a tree of smaller FFT tasks and branch-local diagonal multiplications. A plausible implication is that the algorithm gains not only from reduced asymptotic counts but also from improved locality and finer-grained scheduling.

4. Asymptotic savings and measured behavior

For the split method, the paper gives

$T$1

Defining $T$2, one obtains

$T$3

For memory, the split method stores at most $T$4 vector copies of size $T$5 together with $T$6 Toeplitz data blocks of size $T$7, so

$T$8

and

$T$9

For fully symmetric or skew-symmetric systems, the data storage can be reduced from $d$0 to $d$1, yielding

$d$2

These are the defining quantitative results of FFTMatvec (Siron et al., 2024).

For dimensions $d$3, the asymptotic ratios reported are:

$d$4	$d$5	$d$6	$d$7
2	1.33	1.14	1.25
3	1.71	1.33	1.80
4	2.13	1.52	2.83
5	2.58	1.68	4.71
6	3.05	1.80	8.13

In three dimensions, the general-case memory comparison is stated explicitly as $d$8 for embedding versus $d$9 for split-FFT, i.e. a $t[i_1,\ldots,i_d]$0 reduction. In the symmetric three-dimensional case, the comparison is $t[i_1,\ldots,i_d]$1 versus $t[i_1,\ldots,i_d]$2, i.e. $t[i_1,\ldots,i_d]$3. For a three-dimensional simulation of Green’s function, the measured split-FFT time is reported as approximately $t[i_1,\ldots,i_d]$4 faster, with memory approximately $t[i_1,\ldots,i_d]$5 smaller. The paper also notes that, for moderate dimensions $t[i_1,\ldots,i_d]$6 and large $t[i_1,\ldots,i_d]$7, measured wall-clock speedups often exceed the theoretical $t[i_1,\ldots,i_d]$8 values because smaller FFT sizes can interact more favorably with FFT libraries (Siron et al., 2024).

A recurrent point of clarification is that FFTMatvec does not remove FFTs; it removes the need to realize the full multilevel embedding explicitly. The computational kernel remains FFT plus diagonal multiplication, but reorganized so that projection is performed as early as possible.

5. Parallel structure, implementations, and GPU-oriented extensions

The recursive split exposes $t[i_1,\ldots,i_d]$9 leaf branches, each of which performs FFTs, inverse FFTs, and diagonal multiplications independently. The paper identifies two parallelization strategies: “communicate-and-merge,” in which branches are forked only up to a mid-level and partially merged before continuing serially, and “fully distributed leaves,” in which each leaf is sent to a separate accelerator and results are merged at the root. Because FFT sizes remain $n^d$0 rather than $n^d$1 at each level, performance is described as more consistent on highly tuned libraries such as FFTW and MKL; the lower memory footprint is also stated to reduce cache misses and paging. A Julia implementation at https://github.com/alsirc/SplitFFT_lazyEmbed uses multi-threading and task scheduling, and analogous C/C++/Python implementations are described as straightforward using FFTW or cuFFT together with BLAS/TBB for parallel forking. The stated limitation is that, for very small $n^d$2 or very small $n^d$3, recursive and threading overheads may offset the gains, so a cutoff $n^d$4 below which standard circulant embedding is simpler is recommended (Siron et al., 2024).

A subsequent FFTMatvec HPC application addresses block-triangular Toeplitz matrices and emphasizes performance portability and mixed precision on GPUs. In that formulation, a lower-triangular block-Toeplitz matrix with $n^d$5 block rows and columns is embedded into block-circulant form so that the matvec becomes a convolution,

$n^d$6

reducing the cost from $n^d$7 to $n^d$8. The original code was written in CUDA + cuFFT + cuBLAS; at build time, hipify-perl rewrites CUDA runtime and cuBLAS calls to HIP, after which the build links against HIP, rocFFT, and rocBLAS. No edits to the application source code were needed, because GPU-vendor specialization was pushed into rocBLAS. The mixed-precision workflow assigns each of $n^d$9 kernels a precision template $y=T\,v$0 with $y=T\,v$1, measures

$y=T\,v$2

and selects a configuration from the Pareto front for a prescribed tolerance. Reported single-GPU performance for a block-triangular problem of size approximately $y=T\,v$3M$y=T\,v$4M is approximately $y=T\,v$5, $y=T\,v$6, and $y=T\,v$7 TFLOP/s in double precision on MI250X, MI300X, and MI355X, respectively, and approximately $y=T\,v$8, $y=T\,v$9, and $T$0 TFLOP/s in a mixed configuration with fp32 FFTs and fp64 pointwise operations, corresponding to speedups of $T$1, $T$2, and $T$3. The implementation is reported to scale to $T$4 GPUs on Frontier, with weak-scaling efficiency above $T$5 and strong-scaling efficiency of about $T$6 from $T$7 to $T$8 GPUs, falling to about $T$9 from $v\in\mathbb C^{n^d}$00 to $v\in\mathbb C^{n^d}$01 GPUs because of MPI all-reduce overhead (Venkat et al., 13 Aug 2025).

6. Position within the structured-matrix literature

FFTMatvec belongs to a broader family of fast structured-matrix matvec methods in which Toeplitz, Hankel, or circulant structure is converted into convolution-like operations. Earlier work on Hankel matrix-vector multiplication in multiprecision arithmetic is directly relevant because the paper states that its algorithms are applicable to Toeplitz matrices and to circulant matrices as well. That work distinguishes a direct FFT approach from two alternatives: an FFT-based decomposition method that rewrites multiprecision numbers into standard-precision limbs and performs one enlarged standard-precision convolution, and a Karatsuba-inspired recursive method with complexity $v\in\mathbb C^{n^d}$02. It also reports that, in very high precision, direct multiprecision FFT can be slower than schoolbook multiplication for matrix sizes up to $v\in\mathbb C^{n^d}$03, while the recursive method can outperform FFT-based alternatives for moderate $v\in\mathbb C^{n^d}$04 because it avoids the large overhead of multiprecision FFTs (Beliakov, 2014).

This broader context clarifies the scope of the 2024 split-FFT contribution. FFTMatvec is not merely an assertion that Toeplitz matvec can be accelerated by FFTs; that fact is classical. Its specific contribution is to reorganize the multilevel embedding itself so that the block Toeplitz structure is exploited with less auxiliary expansion, and to do so in a way that is immediately compatible with branch-level parallelization. The subsequent GPU work suggests that the term “FFTMatvec” now functions not only as the name of a particular split-FFT algorithm, but also as a label for a software lineage of FFT-based structured matvec implementations across Toeplitz variants and heterogeneous HPC platforms.