PaperlayerView exact revision on arXiv

Source-linked AI summary

A parallel non-uniform fast Fourier transform library based on an "exponential of semicircle" kernel

Alex H. Barnett, Jeremy F. Magland, Ludvig af Klinteberg

arXiv:1808.06736v2math.NAcs.MS

TL;DR

Efficient type 1 and type 2 NUFFT computation must handle fine-grid spreading, FFT evaluation, and kernel deconvolution. FINUFFT uses the exponential-of-semicircle kernel with numerical quadrature and parallel implementation strategies; its error convergence approaches Kaiser–Bessel rates, while general type 1 error bounds are unavailable when M > N.

  • Problem

    Type 1 and type 2 NUFFTs require fine-grid spreading, FFT computation, and deconvolution, including kernel Fourier-transform evaluation.

  • Method

    FINUFFT uses the exponential-of-semicircle kernel, numerical quadrature for its Fourier transform, phase winding, and adaptive blocking for parallel spreading.

  • Results

    The ES kernel's aliasing error converges almost exponentially with kernel width and approaches the Kaiser–Bessel exponential rate; a factor 2.30 gives slightly lower typical error than factors nearer 2.356.

  • Takeaways & Limitations

    The ES kernel is simpler and faster to evaluate than the Kaiser–Bessel pair while achieving essentially identical error, and reported efficiency gains also accelerate inverse-NUFFT methods.

  • Takeaways & Limitations

    When M > N, no general type 1 relative-error bound exists because the task's relative condition number may be infinite.

Abstract

from arXiv · show

The nonuniform fast Fourier transform (NUFFT) generalizes the FFT to off-grid data. Its many applications include image reconstruction, data analysis, and the numerical solution of differential equations. We present FINUFFT, an efficient parallel library for type 1 (nonuiform to uniform), type 2 (uniform to nonuniform), or type 3 (nonuniform to nonuniform) transforms, in dimensions 1, 2, or 3. It uses minimal RAM, requires no precomputation or plan steps, and has a simple interface to several languages. We perform the expensive spreading/interpolation between nonuniform points and the fine grid via a simple new kernel---the `exponential of semicircle' $e^{β\sqrt{1-x^2}}$ in $x\in[-1,1]$---in a cache-aware load-balanced multithreaded implementation. The deconvolution step requires the Fourier transform of the kernel, for which we propose efficient numerical quadrature. For types 1 and 2, rigorous error bounds asymptotic in the kernel width approach the fastest known exponential rate, namely that of the Kaiser--Bessel kernel. We benchmark against several popular CPU-based libraries, showing favorable speed and memory footprint, especially in three dimensions when high accuracy and/or clustered point distributions are desired.

1. Introduction.

FINUFFT addresses the computational cost of NUFFT exponential sums with a general-purpose parallel library and a new exponential-of-semicircle spreading kernel. The library targets types 1, 2, and 3 in dimensions 1–3, combining near-linear NUFFT algorithms with favorable accuracy, speed, and memory characteristics.

  • FINUFFT implements type 1, type 2, and type 3 transforms for one, two, or three spatial dimensions.Type 1 maps nonuniform data to uniform Fourier modes, type 2 evaluates Fourier series at arbitrary points, and type 3 handles nonuniform frequencies.
  • NUFFT algorithms approximate nonuniform exponential sums in close to linear time rather than the naive O(NM) cost.The transforms support user-specified relative tolerance ε.
  • The standard fine-grid approach spreads or interpolates with a width-w kernel, applies an upsampled FFT, and performs Fourier-space deconvolution, costing O(M| log ε|^d + N log N).Spreading costs O(w^dM), while the FFT costs O(N log N).
  • With σ = 2, Gaussian and B-spline kernels yield approximately 0.5w correct digits, whereas the referenced fastest-known rate exceeds 0.9w correct digits and reaches 12 digits at w = 13.The ES and Kaiser–Bessel kernels have comparable Fourier-transform shapes and outperform the Gaussian’s convergence rate in the figure comparison.
  • FINUFFT avoids precomputation, keeps RAM overhead small, and uses parallelization that adapts to nonuniform point distributions.For medium and large 2D and 3D problems, its speed is competitive with state-of-the-art CPU codes and can be 10× faster at high accuracies.
  • The ES kernel achieves an exponential error convergence rate arbitrarily close to the Kaiser–Bessel rate while being simpler and faster to evaluate.Its Fourier transform is computed by numerical quadrature because no analytic transform is known.

2. Use of the FINUFFT library.

FINUFFT provides callable interfaces for computing NUFFT transforms in C, Fortran, MATLAB, Python, and Julia. Its no-plan design simplifies use but incurs overhead for repeated small problems, motivating additional planned interfaces.

  • The 1D type 1 interface accepts nonuniform points, complex strengths, transform options, and a requested relative precision, writing outputs to a preallocated array.The sign parameter controls the sign of the imaginary unit, and options can enable debugging or bounds checks.
  • FINUFFT exposes similar interfaces from C and Fortran and callable interfaces for MATLAB, Python, and Julia.In high-level interfaces, the number of nonuniform points can be inferred from input sizes.
  • Repeated small problems incur approximately 100 microseconds per thread of FFTW planning overhead because the interface has no plan stage.The affected regime is N and M of order 10^4 or less; multi-input or plan-stage interfaces were being added for this case.
  • Type 1 uses the standard three-step procedure, type 2 reverses those steps, and type 3 combines types 1 and 2.

3. Algorithms.

The algorithms use fine-grid spreading, FFT-based processing, and kernel correction to compute NUFFT types 1 and 2, while type 3 combines type 1 and type 2 operations with additional coordinate handling. Kernel Fourier transforms are evaluated by numerical quadrature, with parameter choices tied to requested precision.

  • Parameter choices: For σ = 2, FINUFFT sets w = ⌈log10 1/ε⌉ + 1 and β = 2.30w, with the first formula targeting the desired accurate digits.A low-upsampling option σ = 5/4 is also available but was not tested.
  • Type 1: The rescaled kernel satisfies ψ(x) = φ(x/α) and ψ̂(k) = α φ̂(αk), enabling correction factors from sampled kernel Fourier transforms.In multiple dimensions, the correction is separable and requires d one-dimensional evaluations followed by an outer product.
  • Type 1: Type 1 spreads nonuniform sources to a fine grid, applies an FFT, then truncates and deconvolves the central Fourier modes.The fine grid uses n = σ^dN points, and spreading uses a periodized kernel of width w.
  • Kernel Fourier transform: Gauss–Legendre quadrature evaluates the kernel Fourier transform without an analytic formula; p ≥ 1.5w + 2 gives errors below ε over |k| ≤ N/2.Because the kernel is real and even, only positive quadrature nodes are needed, and phase winding reduces the cosine work.
  • Type 2: Type 2 reverses type 1: coefficients are pre-corrected and zero-padded, an FFT is performed, and interpolation evaluates the Fourier series at nonuniform targets.Interpolation is the adjoint of spreading and uses wd kernel evaluations per target point.
  • Type 3: Type 3 uses spreading and a type 2 NUFFT, but nonuniform frequency targets require dpN cosine evaluations because no phase-winding shortcut applies.Despite this extra cost, the step consumes only a small fraction of total computation time; for σ = 2, rdyn ≈ e^0.057β and remains at most about 8 for ε ≥ 10^-15.

4. Error analysis and parameter choices.

The analysis derives aliasing-error bounds for the ES kernel, compares their asymptotic rate with Kaiser–Bessel bounds, and uses empirical behavior to justify practical parameter choices. It also explains why observed relative ℓ2 errors can differ from pessimistic worst-case bounds and identifies conditioning and rounding-error limits.

  • Aliasing-error bounds: The type 1 and type 2 aliasing error is represented through a matrix whose entries are controlled by aliased Fourier-transform tails of the spreading kernel.The resulting uniform bound ε∞ leads to ℓ1–ℓ∞ output-error bounds for both transform types.
  • Theoretical convergence: The ES-kernel aliasing bound converges almost exponentially with kernel width w, and its rate approaches the Kaiser–Bessel rate as γ approaches 1 from below.The ES kernel’s Fourier transform has the same below-cutoff exponential form asymptotically, up to algebraic factors.
  • Parameter choices: For σ = 2, the practical factor 2.30 corresponds to γ ≈ 0.976 and gives slightly lower typical error than choosing γ closer to or equal to 1.The factor 2.30 is close to the γ = 1 value 2.356 and is consistent with a factor 2.34 reported for an optimized Kaiser–Bessel kernel.
  • Parameter choices: The kernel width is selected approximately as w ≈ |log ε∞|/(πγ), after dropping the algebraic prefactor in the theoretical bound.The parameter expression is interpreted using the requested tolerance ε as ε∞, with γ ≈ 0.976 for the default setting.
  • Limitations: For M > N, the type 1 task can have infinite relative condition number because the smallest singular value of its transform matrix is zero.Consequently, no general relative-error bound for all nontrivial inputs exists in this regime.
  • Empirical error: Empirical relative ℓ2 error usually tracks the requested tolerance ε, because typical input norms and phase cancellations improve substantially over uniform worst-case bounds.The match holds for the tested routines and random inputs except when rounding error dominates; users are advised to verify convergence for their own data.

5. Implementation issues.

FINUFFT addresses cache use, parallel load balancing, kernel evaluation speed, and numerical conditioning in its implementation. Its spreading scheme adapts to nonuniform point distributions, while sorting and lookup-table approximations improve locality and throughput.

  • Numerical conditioning: Condition numbers can be huge even when M ≤ N, with a random half-interval example yielding a minimum singular value typically of 10^-15.The example uses M=80 and N=100.
  • Cache-aware execution: Sorting nonuniform points to preserve RAM locality typically speeds spreading and interpolation by a factor of 2–10, including sorting time.Each point accesses a block extending ±w/2 grid points in each dimension.
  • Parallel load balancing: Interpolation parallelizes over nonuniform points because each thread reads a fine-grid block and writes a distinct output value.The computation uses weighted sums with tensor-product 1D kernels.
  • Parallel load balancing: Spreading uses sorted point subproblems assigned to threads, avoiding write collisions while adaptively keeping threads busy across nonuniform distributions.Subproblems contain up to 10^4 points, and the scheme requires additional RAM of order the fine-grid size.
  • Kernel evaluation: Kernel evaluation throughput varies substantially with CPU, compiler, and kernel width, while future SIMD support makes the preferred ES-kernel choice uncertain.A Xeon E5-2643 example ranges from 40M to 50–200M evaluations/sec/thread depending on compiler.
  • Kernel evaluation: Piecewise polynomial lookup tables enable SIMD-friendly kernel evaluation, with polynomial degree p = w + 3 sufficient for the stated approximation accuracy.The approach evaluates ordinates within equal-width kernel intervals using centered polynomials.

6. Performance tests.

FINUFFT is benchmarked across dimensions, accuracies, point distributions, thread counts, and competing CPU-based NUFFT libraries. It shows strong parallel performance and favorable speed and memory overhead, especially for large or clustered 3D problems.

  • Parallel scaling: 90% weak-scaling efficiency is achieved through 12 physical-core threads for 3D type 1 and type 2 transforms at ε = 10^-12.Strong-scaling efficiency at 12 threads is 62–74%.
  • Parallel scaling: At ε = 10^-3 and 24 threads, both 3D transforms remain 5–7× faster than a single core.The authors attribute lower efficiency in this case to memory bandwidth because the narrower kernel touches only 4^3 = 64 fine-grid points.
  • 1D comparisons: In 1D, FINUFFT outperforms most tested libraries single-threaded and is 1.5–2× faster than NFFT without precomputation for multithreaded type 1.For type 2, FINUFFT and NFFT with precomputation have similar speed, excluding NFFT’s precomputation time and higher RAM overhead.
  • 2D comparisons: In 2D, FINUFFT is 2–5× faster than CMCL and 4–8× faster than NFFT when single-threaded without precomputation.With NFFT precomputation, multithreaded type 2 speeds are similar, while FINUFFT is 2× faster for type 1 at high accuracy.
  • 3D comparisons: In 3D, FINUFFT is at least 2× faster than other codes for random points except MIRT at 1-digit accuracy, while CMCL is 4–50× slower.For ordered-access spherical quadrature points, CMCL is 2–10× slower, suggesting point ordering affects its performance.
  • 3D comparisons: For large 3D spherical-quadrature type 1 problems, FINUFFT is 8–10× faster than NFFT with precomputation at ε < 10^-6 and at least 10× faster than BART.In a matched large task, FINUFFT is 16× faster than NFFT and 8.6× faster than NFFT with precomputation when precomputation is excluded.
  • Speed and memory: NFFT with precomputation has around 8× FINUFFT’s RAM overhead, while FINUFFT avoids precomputation through on-the-fly kernel evaluation.The library’s conclusions emphasize order-of-magnitude speedups and lower RAM overhead for certain clustered 3D problems.

7. Conclusion.

The conclusion identifies future benchmarking and interface work, and notes a kernel limitation for applications requiring spatial derivatives. A smooth variant addresses endpoint behavior but adds computational cost.

  • Future work: Future work includes benchmarking type 3 transforms, completing an interface for repeated small problems, and developing a carefully benchmarked general-purpose GPU NUFFT.Implementation of the repeated-small-problem interface and GPU code is in progress.
  • Kernel limitations: The spreading kernel has unbounded derivatives, limiting its suitability when spatial derivatives are required in particle-mesh Ewald applications.The unbounded derivatives arise from an inverse square-root singularity at the endpoints.
  • Kernel limitations: A smooth exponentially-close kernel variant avoids the endpoint singularity but requires one extra reciprocal or an approximation.The variant is smooth up to the endpoints.
  • Acknowledgments: The acknowledgments credit discussions with contributors from numerical analysis, imaging, and computational science.The paper also notes that the Flatiron Institute is a division of the Simons Foundation.
Loading 1808.06736v2…