Indirect Convolution Algorithms

• Indirect convolution algorithms are techniques that compute discrete convolutions using mathematical transforms or reductions to auxiliary algebraic problems instead of direct sliding-sum implementations.
• They encompass methods like FFT-based, im2col+GEMM, pointer indirection, sparse matrix approaches, and Winograd, each balancing computation, memory, and latency trade-offs.
• These methods are pivotal in fields such as deep neural networks, real-time audio processing, statistical physics, and scientific computing, enhancing efficiency and scalability.

Indirect convolutions are algorithmic strategies that compute discrete convolutions through transformations, intermediary representations, or reduction to auxiliary algebraic problems, as opposed to direct sliding-sum implementations. The class of indirect algorithms includes FFT-based convolution, matrix-reduction methods (e.g., im2col+GEMM), pointer-indirection schemes, sparse-matrix-multiplication formulations, fast bilinear/WINOGRAD constructions, multipole-like expansions, and real-time block-free approaches (e.g., taches-algorithm). Indirect algorithms dominate modern high-performance computation of convolution, especially in contexts such as deep neural networks, statistical physics, large-scale scientific computing, and real-time audio signal processing, due to their favorable complexity, memory, or latency characteristics.

1. Fundamental Principles and Categories

Indirect convolution algorithms can be categorized based on their mathematical transformation and data-movement paradigm:

• Frequency-domain/FFT convolution: Relies on the convolution theorem, mapping convolution in the time or spatial domain to pointwise multiplication in the frequency domain using the Discrete Fourier Transform (DFT) computed efficiently by FFT. This technique underpins large-scale linear time-invariant filtering and probabilistic distribution summation (Ruckdeschel et al., 2010, McEwen et al., 2016).
• Matrix-based reductions (im2col+GEMM and its variants): Reshape (lower) the input tensor into a patch matrix such that convolution becomes a batched matrix-matrix multiply (GEMM), leveraging highly-tuned BLAS implementations; variants such as MEC (memory-efficient convolution) and pointer-indirection reduce memory overhead (Cho et al., 2017, Dukhan, 2019, Juan et al., 2020).
• Pointer indirection: Constructs a layout-agnostic buffer (array of pointers) referencing data-rows filling the "A" matrix for GEMM, avoiding materialization of the patch matrix as in classical im2col; enables lower memory workspaces and reduced data movement (Dukhan, 2019).
• Sparse matrix-vector convolution: Represents convolution operators as large, sparse matrices acting on vectorized input data via SpMV. Exploits the inherent sparsity caused by padding and stride, pruning operations involving input zeros (Chaudhry, 2024).
• Fast bilinear/Winograd algorithms: Reformulate convolution as a sequence of linear transforms (A-transform, B-transform), elementwise multiplication (Hadamard product), and an output transform (C-transform), minimizing arithmetic count—most notably for small-channel, small-filter deep neural network layers (Ju et al., 2019).
• Tree-based decompositions: In integral convection problems (e.g., in computational physics), tree-based descent-only or multipole-like expansions compress far-field contributions to achieve quasilinear complexity—e.g., the fast free memory (FFM) algorithm (Aussal et al., 2019).
• Time-domain accumulator/buffer schemes: E.g., the taches-algorithm accumulates delayed, scaled impulse responses ("taches") in a circular buffer, enabling real-time, low-latency convolution without transforming or blocking operations (Millot et al., 2024).

2. Core Algorithms and Mathematical Structures

Each major class admits distinct algorithmic structure and memory patterns:

Method	Transformation	Auxiliary Data	Workspace Overhead
FFT convolution	DFT & pointwise mult., IFFT	Frequency-domain buffers	O(m log m) for length m
Im2col+GEMM	Patch-lowering, GEMM	Col-unfolded input matrix	O(C·K_h·K_w·output elems)
Pointer Indirection	Patch-lookup buffers, GEMM	Indirection pointer array	O((K_h·K_w·output elems) / C)
SpMV approach	Sparse matrix on vectorized input	Sparse index matrices	O(output elems·kernel size)
Winograd	Bilinear transform approach	A, B, C transform matrices	O(block size) per block
Taches/circular buf.	Accumulate delayed tracks	Circular buffer (N_h size)	O(N_h) per channel
FFM/tree-compressed	Lagrange, ACA, or multipole	Interpolation + transfer	O(N), no global matrix

• FFT convolution: For length-m signals, classical FFT convolution costs O(m log m) and achieves high arithmetic efficiency for long filters or signals. Combination of zero-padding and smoothing maintains accuracy for continuous distributions (Ruckdeschel et al., 2010).
• Im2col+GEMM: Translates convolution into a dense GEMM by arranging all sliding window patches into columns; reduces convolution to matrix multiply but with a large memory prologue. The MEC variant reduces the lowering buffer size by a factor 1/K_h by exploiting the redundancy in overlapping patches (Cho et al., 2017).
• Pointer-indirection algorithms: Instead of lowering the entire patch matrix, an indirect buffer of pointers or indexes to underlying data is constructed, allowing microkernel routines to load data on-the-fly, bypassing the im2col bottleneck. This approach can reduce overall run time by up to 62% for non-trivial kernel shapes at moderate channel counts (Dukhan, 2019, Juan et al., 2020).
• SpMV-based convolution: Construction of a sparse matrix T = C·P representing the full strided, padded convolution allows highly sparse mat-vec computation for any input. Closed-form expressions for the number of non-zero multiplications allow precise flop counting and facilitate usage of optimized SpMV routines (Chaudhry, 2024).
• Fast bilinear/Winograd convolution: For small kernel and block sizes, e.g., F(2,3), a sequence of three linear transformations and elementwise multiplications reduces the multiply count sharply at the expense of more additions, outperforming direct convolution for small filters in DNN accelerators (Ju et al., 2019).
• Taches/circular-buffer convolution: Maintains an accumulator buffer of length equal to the filter, updates it with each input sample by adding e[n]·h[·], and reads out the output with sample delay. Enables real-time operation with per-sample latency and immediate response to impulse response changes (Millot et al., 2024).
• FFM and tree methods: For volume or surface convolution, split interactions using geometric multiresolution (octree), compress admissible box-box blocks via Lagrange or ACA, and only store/interpolate per required box; this allows O(N log N) complexity with O(N) peak storage for N points (Aussal et al., 2019).

3. Computational Complexity and Memory Characteristics

The computational and memory profile of indirect convolutions strongly depends on the problem structure and method:

• FFT/DFT-based methods: O(m log m), excellent for long signals but require full bufferization and are less suited for very low-latency, real-time systems. High arithmetic intensity and cache locality.
• GEMM-based reductions: Arithmetic cost is O(#outputs × kernel size × input channels), but memory can become dominant due to buffering and lowering; MEC or pointer-indirection reduces this penalty (Cho et al., 2017, Dukhan, 2019).
• Pointer-indirection: Memory overhead is O(1/C) times that of im2col; pointer arrays are cached efficiently, especially when C is large, and microkernels can be fused with GEMMs (Dukhan, 2019).
• SpMV-based approaches: Assembly of the effective sparse matrix is O(#outputs × k²) but can be amortized over repeated convolutions; running each SpMV is exactly as sparse as the nonzero coverage of the kernel on non-padded data—particularly advantageous for large paddings or small kernels (Chaudhry, 2024).
• Winograd: For F(2,3), reduces multiplies per output from 3 (direct) to 2, increasing additions to 6 (with transforms); O(N) overall but with smaller leading constant for multiplies (Ju et al., 2019).
• Taches/circular-buffer: O(N_h) per output sample, negligible sample latency, extremely favorable for moderate filter lengths and dynamic impulse response (Millot et al., 2024).
• Tree-based/FFM: O(N log N) for smooth kernels and O(N log²N) for oscillatory kernels; storage strictly O(N), enabling evaluation on up to 10⁹ unknowns (Aussal et al., 2019).

4. Applications and Real-World Benchmarking

Indirect convolution algorithms are deployed across a spectrum of domains, each leveraging unique features:

• Deep learning (CNN inference and training): Im2col+GEMM, pointer-indirection, MEC, and Winograd algorithms are foundational for modern deep neural network accelerators (server, mobile, embedded). Indirect approaches dominate for memory-limited devices and small-filter, large-batch operations (Cho et al., 2017, Dukhan, 2019, Juan et al., 2020, Ju et al., 2019).
• Probabilistic modeling and statistics: FFT-based indirect convolution supports high-precision summation of arbitrary independent distributions (e.g., in R's distr package), outperforming recursion-based methods at large support sizes; highly accurate to sub-ε levels (Ruckdeschel et al., 2010).
• Signal processing and audio: Taches-algorithm supports real-time audio convolution (≤15,000 taps at 44.1 kHz on modest CPUs), real-time impulse response variation, and low-latency plugins (Millot et al., 2024).
• Computational physics: FFM and similar tree-based indirect convolutions enable O(N log N) integral evaluations for large-scale scattering, potential, or electromagnetic boundary-integral problems—enabling scales up to one billion points (Aussal et al., 2019).
• Sparse convolution and inference: SpMV method is competitive or faster on CPU and similar in throughput on GPUs for architectures exploiting high sparsity or repeated convolution with fixed kernel (Chaudhry, 2024).
• Scientific computation and perturbation theory: FAST-PT algorithm applies indirect FFT-based convolution for nonlinear cosmological perturbation calculations, enabling routine embedding in sampling loops (McEwen et al., 2016).

Reported benchmarks show:

• Pointer-indirection outperforms im2col+GEMM by up to 62% in total runtime for non-1x1 convolutions; MEC cuts memory usage by a factor of 2–4× and can yield 1.2–1.6× speedups on common CPUs/ARMs (Cho et al., 2017, Dukhan, 2019, Juan et al., 2020).
• SpMV achieves up to 2× speedup over im2col+GEMM on CPUs in practical deep learning models such as DenseNet121 (Chaudhry, 2024).
• FFM handles 1e9 unknowns in 3.1 hours on 12 CPU cores with 100 GB RAM; kernel-independent FFM outperforms classical FMM/H-matrix for sufficiently large N (Aussal et al., 2019).

5. Algorithmic Trade-offs and Limitations

Indirect convolution algorithms present nuanced trade-offs:

• FFT-based methods: Optimal for large filters/signals but incur block and buffering latency; non-trivial for very short filters or real-time, streaming systems due to block structure.
• Im2col+GEMM: Maximal compute reuse but costly in auxiliary memory and data copying; indirect (pointer) buffer variants ameliorate this penalty at increased code complexity.
• Pointer-indirection: Best for NHWC layouts; on 1×1 stride-1 convolutions direct GEMM is slightly faster, while im2col overhead is avoided for other parameterizations. Not ideal for depthwise or transposed convolution modes (Dukhan, 2019).
• SpMV approaches: Efficient only if kernel sparsity or padding aligns to skip large numbers of zero multiplications; pointer management and sparse indexing overhead can dominate for small matrices or high-throughput GPU implementations (Chaudhry, 2024).
• Winograd/bilinear algorithms: Favorable for small kernels with block sizes matching hardware, but error grows with transform size and ill-conditioned matrices can cause numerical instability for deep transforms or large filters (Ju et al., 2019).
• Taches/circular-buffer: Per-output cost scales with impulse response length N_h; for N_h ≫ 10^4, FFT-based methods dominate (Millot et al., 2024).
• FFM/tree-based: Optimal for smooth kernels and uniform domains; setup time for the octree/scaling interplay becomes non-negligible for rapidly changing/fine-grained source-target distributions (Aussal et al., 2019).

6. Methodological and Implementation Considerations

Advanced indirect convolution methods often exploit hardware or domain-specific optimizations:

• Blocking and tiling: Vertical and horizontal tiling can further reduce buffer sizes and tune GEMM calls to hardware preferences (cf. MEC and BLIS indirect GEMM) (Cho et al., 2017, Juan et al., 2020).
• Object-oriented dispatch: FFT-based convolution is embedded in object-oriented frameworks for generalized probabilistic convolution, with automatic fallback to analytic results or specialized routines (Ruckdeschel et al., 2010).
• Batching: Processing multiple convolution outputs simultaneously (e.g., via Winograd blocks or SpMV batches) optimally utilizes vector and matrix units.
• Real-time updates: Buffer-accumulation algorithms (taches) support per-sample filter adaptation, which is nearly intractable with block-based FFT convolution.
• Low-memory regimes: FFM grants strict O(N) memory profile, meeting computational physics needs where N greatly exceeds available RAM.

7. Future Directions and Theoretical Implications

Ongoing work investigates:

• Minimal lowering for indirect GEMM: Rethinking the mapping from convolutional operations to linear algebra primitives to achieve strict lower bounds in memory use, guided by compact-lowering as in MEC (Cho et al., 2017).
• Generalized tiling and auto-tuning: Tailoring block sizes and decomposition for optimal cache performance and parallel efficiency across architectures.
• Hybrid time-frequency methods: Combining sample-wise and blockwise convolution to exploit both low-latency and high-throughput regimes, e.g., in artificial reverb and interactive audio.
• Sparse extensions in DNN kernels: Leveraging sparsity in kernel matrices (weights) for acceleration in inference and training.
• Tree-free kernel decompositions: Developing methods that forego explicit hierarchical trees in favor of Fourier or other analytic decompositions for general convolution kernels (Aussal et al., 2019).
• Dynamic indirect methods: Algorithms capable of adaptively switching between direct, indirect, and hybrid schemes based on real-time workload and hardware profiling.

Indirect convolution algorithms thus form a foundational, highly-adaptive computational toolkit with tunable memory, arithmetic, and latency characteristics, and continue to evolve to match the demands of large-scale, sparse, or dynamic computational regimes (Cho et al., 2017, Dukhan, 2019, Juan et al., 2020, Ruckdeschel et al., 2010, McEwen et al., 2016, Chaudhry, 2024, Aussal et al., 2019, Ju et al., 2019, Millot et al., 2024).