Custom Non-Uniform FFT Implementation

• Custom non-uniform FFT is a method for fast Fourier transforms on irregularly sampled data using specialized gridding and interpolation techniques.
• It employs advanced algorithms such as multipole expansions and kernel fusion to achieve quasi-linear performance and effective error control.
• Applications span mesoscopic simulation, photoacoustic tomography, and Fourier neural operators, illustrating its broad utility in scientific computing.

A custom non-uniform fast Fourier transform (NUFFT) implementation encompasses a family of algorithmic and software strategies for evaluating Fourier transforms of data sampled at irregular (non-equispaced) points, with performance and accuracy comparable to the standard FFT on regular grids. Unlike the traditional discrete Fourier transform (DFT), whose O(N²) time complexity is prohibitive for large N and whose acceleration via the FFT requires data on uniform grids, the NUFFT and its variants are designed specifically to scale efficiently for nonuniform spatial or frequency sampling. These custom implementations form the computational backbone in fields as diverse as mesoscopic simulation, spectroscopic data analysis, photoacoustic tomography, convolutional neural operators, and large-scale parallel computing.

1. Mathematical Principles Underlying Non-Uniform FFT

The NUFFT generalizes the DFT by enabling fast computation of the transformation

$$F_k = \sum_{j=1}^{N} f(x_j) \exp(-i 2\pi v_k x_j),\quad k=1,\ldots,M$$

where $\{x_j\}$ are generally nonequispaced sample points, and $\{v_k\}$ are frequency grid nodes—which themselves may also be nonuniform (Type 3 NUFFTs). Standard DFT and FFT algorithms require $x_j$ to be equispaced and $v_k$ to lie on an integer frequency grid, mapping naturally to the circulant algebra underlying the FFT. In the NUFFT, explicit algorithms are developed to circumvent this restriction, typically by decomposing the kernel or convolution into manageable sub-problems.

Algorithmic approaches commonly encountered include:

• Gridding or Interpolation-Based Methods: The data are spread onto a regular grid using a localized window function (e.g., Gaussian, Kaiser–Bessel, exponential of semicircle), followed by standard FFT and deconvolution to correct for kernel effects (Barnett et al., 2018, 2208.00049).
• Frame Theoretic and Convolutional Gridding: The approximation of the inverse transform is put on a rigorous footing by the application of frame theory, leading to optimal density compensation weights for interpolation and nonuniform sampling (Gelb et al., 2014).
• Pruned and Fused Kernel Methods: Architecturally optimized implementations, such as TurboFNO, integrate frequency truncation, zero-padding, and spectral filtering directly into the FFT kernel itself, fusing FFT, GEMM, and iFFT stages for end-to-end efficiency in applications like Fourier neural operators (Wu et al., 16 Apr 2025).
• Multipole Expansions: For large-scale problems with strong nonuniformity, the FFT kernel is factorized, and the “slow-varying” term is accelerated by fast multipole summation tailored for periodic domains (Gumerov et al., 2016).

The selection and parameterization of window functions directly govern the trade-off between computational efficiency and numerical error, with exponential decay of aliasing error as a function of kernel width—a property exploited for user-tunable accuracy (Barnett et al., 2018, 2208.00049).

2. Algorithmic Structure and Types of Custom Implementations

Custom NUFFT implementations are classified according to the mapping between input and output grids:

NUFFT Type	Input grid	Output grid
1	Nonuniform	Uniform
2	Uniform	Nonuniform
3	Nonuniform	Nonuniform
4/5	Inverse of 1/2	Inverse of 1/2

In practice, efficient implementation involves decomposition into the following stages:

• Spreading/Gridding: Map nonuniform data to a fine uniform grid using a window function, minimizing support to control computational footprint (Barnett et al., 2018, Shih et al., 2021).
• FFT or iFFT: Perform standard (batch) FFTs on the upsampled grid; this is the primary source of acceleration.
• Deconvolution/Correction: Compensate the transform for the window’s effect in Fourier space, either by analytic division (if the window is well-behaved) or via efficient quadrature (Barnett et al., 2018).
• Adjoint/Interpolation: For Type 2/Type 3 transforms, interpolate the uniform FFT back to irregular points, reusing the window function.

Variants further integrate nonuniformity in both data and frequency, as required in image reconstruction (e.g., MRI), scalar diffraction, photoacoustic imaging, and PDE operator learning.

Recent research has shown that the classical iterative solvers for inverse problems (Type 4/5) can be replaced with explicit non-iterative convolutional and Lagrange interpolation-based methods, dramatically reducing computational overhead (Selva, 2016).

In massively parallel and GPU contexts, custom strategies such as cache-aware blocked spreading, shared-memory binning, and swizzled access to shared memory banks eliminate memory bottlenecks, leading to dramatic throughput improvements, especially in the high-dimensional setting (Shih et al., 2021, Wu et al., 16 Apr 2025).

3. Parameter Selection and Error Control

A key theme in custom NUFFT design is parameter optimization for accuracy and speed:

• Window Parameters: The kernel width and type are set according to the desired accuracy; for exponential of semicircle and Kaiser–Bessel kernels, window width $w = \lceil|\log_{10}\epsilon|\rceil + 1$ ensures uniform-norm error $\lesssim \epsilon$ for a target tolerance $\epsilon$ (Barnett et al., 2018, 2208.00049).
• Oversampling Factor: Choosing the oversampling $\sigma$ (typically $\sim 2$) balances aliasing and kernel support size.
• Truncation and Pruning: For kernel expansions or butterfly stages in the FFT (especially in fused GPU kernels), pruning of high-frequency bands that are not used in subsequent computation can yield significant computational savings, as observed in TurboFNO, where only a prescribed fraction (e.g., 25%) of FFT outputs are produced and processed (Wu et al., 16 Apr 2025).

Specific applications—such as Ewald summation in DPD (Wang et al., 2013)—require additional parameters: Ewald splitting parameter $\alpha$, reciprocal and real space cutoffs, and charge smearing corrections. These are derived analytically to balance convergence and ensure that both energy and force calculations are accurate within a prescribed threshold.

Robust frame-theoretic algorithms further permit principled, global error control by optimizing density compensation factors via least-squares or pseudoinverse solutions, supported by proven convergence rates (Gelb et al., 2014).

4. Numerical Efficiency and Scaling Properties

Custom NUFFT implementations achieve quasi-linear complexity, typically $O(N \log N)$ for $N$ sample points, in contrast to the $O(N^{3/2})$ or $O(N^2)$ scaling of direct computation for nonuniform grids.

Empirical benchmarks in both CPU and GPU settings support these scaling claims:

Implementation	Key Features	Scaling	Notable Performance
FINUFFT (Barnett et al., 2018)	ES kernel, load balance	$O(N\log N)$	Up to 8–10× faster in 3D imaging
cuFINUFFT (Shih et al., 2021)	GPU sharing, blocking	$O(N\log N)$	Throughput $>10^9$ pts/s on V100
TurboFNO (Wu et al., 16 Apr 2025)	Fused FFT-GEMM-iFFT	$O(N)$ (with pruned I/O)	Up to 150% speedup over cuFFT+cuBLAS
FMM-based NUFFT (Gumerov et al., 2016)	Multipole expansion	$O(N)$ at machine precision	Rapid evaluation for large N

Performance is sustained even for highly clustered or irregular sample distributions by partitioning data into blocks, optimizing for cache and memory access patterns, and, when relevant, overlapping communication with computation (in MPI settings (Balty et al., 2022)). Block partitioning and bin-based gridding in shared memory contribute to eliminating write collisions and optimizing bandwidth utilization, as documented for cuFINUFFT and TurboFNO.

5. Representative Applications in Scientific Computing

Custom non-uniform FFT implementations underpin a diverse set of application domains:

• Mesoscopic Simulation: In ENUF-DPD, nonuniform FFT computes contributions to Ewald-split electrostatics with O(N log N) scaling, enabling simulations of polyelectrolyte conformations and charged dendrimer–membrane systems with high fidelity (Wang et al., 2013).
• Photoacoustic Tomography: NEDNER-NUFFT reconstructs images from arbitrarily placed sensors, with dramatic improvements in axial and lateral resolution compared to interpolation-based FFT methods, and natural exclusion of limited-view artifacts by optimal sensor placement (equi-steradian/equiangular) (Schmid et al., 2015, Schmid et al., 2015).
• Fourier Transform Spectrometry: Direct application of NUFFT avoids spectral distortions or amplitude loss under undersampling, and suppresses aliasing present in traditional interpolated FFT approaches, with superior computational and noise properties (Wen et al., 2022).
• Time Series Periodogram Estimation: Lomb–Scargle periodograms recast as Type 1 NUFFT yield orders-of-magnitude improvements in speed and accuracy over extirpolation approaches, with seamless GPU acceleration and drop-in integration in established analysis pipelines (Garrison et al., 12 Sep 2024).
• Fourier Neural Operators: TurboFNO fuses FFT, GEMM, filtering and iFFT in spectral learning operators, eliminating redundant memory steps and focusing computation on spectral bands actually used by the neural model (Wu et al., 16 Apr 2025).
• Finite Element and Boundary Integral Evaluation: Efficient treatment of singular convolution kernels (e.g., single-layer potentials) on unstructured grids via NFFT enables O(N log N) evaluation in micromagnetics and beyond (Exl et al., 2013).

6. Architectural and Software Design Considerations

Custom NUFFT implementations exploit hardware and domain-specific structure:

• Multithreading and SIMD Vectorization: Kernel computations in FINUFFT and NFFT.jl are vectorized and thread-aware, enabling strong and weak scaling on CPU architectures, and efficient use of modern high-level languages with type and dimension genericity (Barnett et al., 2018, 2208.00049).
• Cache-Aware, Load-Balanced Spreading: Data are sorted into spatial bins to improve memory access, and, for GPU codes, blocks in shared memory minimize global atomic conflicts (Shih et al., 2021).
• Shared Memory Swizzling Patterns: For GPU kernels in TurboFNO, two-tiered swizzling of shared memory bank access eliminates conflicts between FFT output, GEMM operation, and iFFT, ensuring 100% bandwidth utilization (Wu et al., 16 Apr 2025).
• Distributed and Parallel Infrastructure: Libraries such as FLUPS provide flexible data layout (cell-centered vs. node-centered), tolerant communication patterns (blocking, non-blocking, and custom MPI datatypes), and rigorous validation for massively parallel FFT and, by extension, custom NUFFT implementations (Balty et al., 2022).

Open source codebases with BSD-style licenses, modularity, and language bindings across C/C++/Julia/MATLAB/Python/Fortran (e.g., FINUFFT, NFFT.jl, cuFINUFFT) facilitate integration into broader scientific workflows and further customization to application-specific requirements.

7. Limitations, Challenges, and Ongoing Directions

Significant challenges in custom NUFFT implementations include:

• Error Control in Highly Oscillatory and High-D Grid Regimes: The approximation error from windowing and kernel truncation can saturate due to floating point limitations, and expansions must be tailored to the problem’s spectral properties. Parameter tuning (e.g., truncation order, window width) remains nontrivial for new domains (Gumerov et al., 2016, 2208.00049).
• Inverse and Multidimensional Extensions: Explicit, non-iterative type 4/5 algorithms are currently restricted to one-dimensional settings; multidimensional generalization requires new mathematical tools beyond the reach of univariate Lagrange interpolation (Selva, 2016).
• Efficiency for Highly Clustered Nonuniform Data: While recent GPU and parallel implementations mitigate the computational cost, performance still degrades or becomes memory-bound for pathological (e.g., singularly clustered) sampling distributions without careful load balancing (Shih et al., 2021, Barnett et al., 2018).
• Integration with Domain-Specific Pipelines and PDE Solvers: The translation from theoretical advances to robust, high-performance application pipelines remains ongoing, especially for domains requiring nonstandard boundary conditions, convolutions with highly singular kernels, or explicit coupling to adaptive mesh refinement (Exl et al., 2013, Schmid et al., 2015, Balty et al., 2022).

A plausible implication is that ongoing research will focus on making NUFFT routines ever more architecture- and problem-specific, integrating hardware-aware optimizations (e.g., kernel fusion, dataflow alignment, memory swizzling), advanced density compensation, and opportunistic exploitation of task parallelism, to meet the demands of increasingly heterogeneous and high-dimensional scientific workloads.

---

In summary, custom non-uniform FFT implementations constitute a mathematically sophisticated, algorithmically diverse, and application-driven frontier in computational science. The progress in this field is characterized by the synergy of mathematical rigor (e.g., frame theory, multipole expansions), hardware-specific optimizations (e.g., kernel fusion, shared memory tuning), and validation in challenging real-world applications from materials simulation to machine learning and imaging. The current literature demonstrates both broad utility across domains and a trend toward truly plug-and-play, scalable, and high-accuracy solutions adaptable to the next generation of scientific and engineering challenges.