Type 3 NUFFT: Arbitrary Grid Fourier Transform

• Type 3 NUFFT is a nonuniform-to-nonuniform discrete Fourier transform that handles arbitrary spatial and frequency grids for flexible spectral analysis.
• It employs fast algorithms using gridding, oversampling, and kernel-based interpolation to significantly reduce the computational complexity of the direct transform.
• Practical implementations use low-rank approximations and hierarchical solvers to achieve quasi-linear complexity, benefiting applications in imaging, MRI, and numerical simulations.

The Type 3 Nonuniform Fast Fourier Transform (NUFFT)—also referred to as the "nonuniform-to-nonuniform" transform—computes the discrete Fourier transform (DFT) between a set of arbitrary nonuniform spatial (or temporal) input positions and another set of nonuniform frequency points. Distinguished from Type 1 (nonuniform to uniform frequency grid) and Type 2 (uniform to nonuniform frequency grid) NUFFT, Type 3 is the most general and versatile, enabling the transformation of data sampled on arbitrary grids to arbitrary target frequency locations. It is central to a wide range of applications requiring flexible sampling and spectral analysis and has motivated substantial algorithmic research due to the computational complexity of the direct (quadratic-time) approach.

1. Formal Definition and Mathematical Structure

Let $x_j \in \mathbb{R}^d$ for $j=1,\ldots,M$ be the nonuniform input points, $f(x_j)\in\mathbb{C}$ the data at those points, and $\omega_k \in \mathbb{R}^d$ for $k=1,\ldots,K$ the nonuniform target frequencies. The Type 3 NUFFT computes

$$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$

In one-dimensional scenarios, this reduces to $f_j = \sum_{k=1}^{N} u_k\,e^{2\pi i x_j \omega_k}$ for given nonuniform $\{ x_j \}$ and $\{ \omega_k \}$, encapsulating the core algebraic structure (Li et al., 3 Dec 2025, Sultan et al., 9 Jan 2025).

The corresponding transformation matrix $A$ with entries $j=1,\ldots,M$0 is typically full-rank with no exploitable Toeplitz or circulant structure, making naive direct computation an $j=1,\ldots,M$1 operation. The inverse problem, often overdetermined (i.e., $j=1,\ldots,M$2), is commonly solved in a least-squares sense.

2. Fast Algorithms and Kernel-Based Approximations

To achieve computational tractability, fast Type 3 NUFFT algorithms combine clever interpolation (gridding), oversampling, and windowing techniques with uniform FFTs:

1. Gridding/Spreading: Input data is spread onto a uniform oversampled grid using a localized kernel $j=1,\ldots,M$3 (e.g., Kaiser–Bessel or exponential-of-semicircle functions).
2. Uniform FFT: A standard FFT of the oversampled grid is performed; this leverages highly optimized FFT implementations for uniform data.
3. Degridding/Regridding: The transform is evaluated at nonuniform target frequencies using interpolation with the same kernel, together with an explicit deconvolution by the kernel's Fourier transform (Sultan et al., 9 Jan 2025, Barnett et al., 2018).

Parameter selection—particularly the oversampling ratio $j=1,\ldots,M$4 and kernel width $j=1,\ldots,M$5—controls accuracy and performance. For $j=1,\ldots,M$6-dimensional data, the total computational cost is $j=1,\ldots,M$7 for balanced regimes ($j=1,\ldots,M$8), with the error controlled by $j=1,\ldots,M$9 for an appropriate choice of kernel $f(x_j)\in\mathbb{C}$0.

Libraries such as FINUFFT implement the spreading and degridding steps using the exponential-of-semicircle (ES) kernel. ES offers compact support and exponential alias suppression, facilitating minimal RAM usage and high throughput for large-scale 1D, 2D, and 3D problems. The underlying error analysis is precise, with kernel support $f(x_j)\in\mathbb{C}$1 for target accuracy $f(x_j)\in\mathbb{C}$2 (Barnett et al., 2018).

3. Low-Rank and Structured Matrix Decompositions

Another approach exploits numerical low-rank structure in the Type 3 NUDFT matrix. Key steps:

• Matrix Factorization: Each entry can be expressed as a product of a smooth low-rank "quotient" and a uniform DFT entry.
• Low-Rank Approximation: The smooth factor is approximated via a bivariate Chebyshev expansion, yielding a rank-$f(x_j)\in\mathbb{C}$3 separation (Ruiz-Antolin et al., 2017).
• Implementation: The transform is evaluated as a sum of $f(x_j)\in\mathbb{C}$4 diagonally-scaled FFTs, with an overall cost of

$f(x_j)\in\mathbb{C}$5

The approximation error is bounded entrywise by $f(x_j)\in\mathbb{C}$6, with total $f(x_j)\in\mathbb{C}$7-error $f(x_j)\in\mathbb{C}$8.

This approach is robust in higher dimensions and can be extended to 2D by dimension-by-dimension low-rank approximations, leading to a sum of $f(x_j)\in\mathbb{C}$9 2D FFTs for required accuracy.

4. Direct and Inverse Type 3 NUDFT Solvers

Superfast direct solvers for the Type 3 inverse NUDFT (INUDFT) have been proposed based on hierarchical matrix techniques:

• Matrix Factorization: The NUDFT-III matrix $\omega_k \in \mathbb{R}^d$0 is decomposed approximately as $\omega_k \in \mathbb{R}^d$1, where $\omega_k \in \mathbb{R}^d$2 is a Type II NUDFT matrix and $\omega_k \in \mathbb{R}^d$3 an explicit matrix embodying the residual structure. Both $\omega_k \in \mathbb{R}^d$4 and $\omega_k \in \mathbb{R}^d$5 are amenable to hierarchical semiseparable (HSS) approximation (Li et al., 3 Dec 2025).
• Forward and Inverse Application: Forward application leverages fast Type II NUDFT algorithms and HSS matrix multiplication. The inversion proceeds efficiently due to quasi-linear complexity of HSS-based direct solvers, with the step $\omega_k \in \mathbb{R}^d$6 requiring only one type-II direct solve and one HSS solve.

Given $\omega_k \in \mathbb{R}^d$7 and for i.i.d. uniform samples, Frobenius-norm error bounds for the approximation $\omega_k \in \mathbb{R}^d$8 and $\omega_k \in \mathbb{R}^d$9 are sharp, and in practice the truncation parameter $k=1,\ldots,K$0 suffices for small errors.

Numerical performance exhibits $k=1,\ldots,K$1 construction and factorization, and $k=1,\ldots,K$2 per-solve runtime. As a preconditioner for conjugate gradients on the normal equations, the method dramatically accelerates convergence—reducing iteration counts from $k=1,\ldots,K$3 (unpreconditioned) to $k=1,\ldots,K$4–$k=1,\ldots,K$5 (preconditioned) (Li et al., 3 Dec 2025).

5. Error Estimates, Parameter Choices, and Parallel Performance

Rigorous error analysis underpins kernel-based and low-rank approaches:

• Kernel Interpolation Error: For a kernel of half-width $k=1,\ldots,K$6 and oversampling $k=1,\ldots,K$7, the error decays exponentially as $k=1,\ldots,K$8, with $k=1,\ldots,K$9 and $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$0 depending on the kernel (Sultan et al., 9 Jan 2025, Barnett et al., 2018).
• Low-Rank Approximation Error: The Chebyshev expansion rank $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$1 is chosen so $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$2, where $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$3 is a slowly growing function of $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$4 (Ruiz-Antolin et al., 2017).
• Hierarchical Compression Error: HSS approximation achieves quasi-linear complexity with the rank $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$5 for target accuracy $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$6 (Li et al., 3 Dec 2025).

Optimized libraries exploit cache and memory access patterns, multithreading, and on-the-fly kernel evaluation. For example, FINUFFT achieves 80–90% parallel efficiency up to one thread per physical core, with no plan/precompute stage and minimal RAM escalation (Barnett et al., 2018).

6. Applications and Practical Implementations

Type 3 NUFFT is pivotal for scientific computing tasks where both data and reconstruction loci are aperiodic or non-Cartesian:

• Non-Line-of-Sight (NLOS) Imaging: Type 3 NUFFT enables flexible reconstruction and sampling strategies in phasor-field-based NLOS imaging, accommodating arbitrary, irregular sensor and reconstruction grids. This flexibility is crucial for compressed measurements and scalable 3D reconstruction (Sultan et al., 9 Jan 2025).
• Signal Processing and Medical Imaging: Arbitrary k-space trajectories in MRI and advanced image reconstruction techniques necessitate Type 3 NUFFT algorithms for efficient and accurate resampling.
• Numerical PDE Solvers: The requirement to map between spatially and spectrally nonuniform discretizations arises in high-fidelity scientific simulations.

Contemporary open-source implementations such as FINUFFT (Barnett et al., 2018) provide highly efficient, high-accuracy routines for Type 3 transforms in dimensions 1–3, with demonstrated throughput and robustness for both scattered and clustered data.

7. Summary Table: Algorithmic Approaches for Type 3 NUFFT

Approach	Main Reference	Complexity
Kernel-based (gridding)	(Sultan et al., 9 Jan 2025, Barnett et al., 2018)	$$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$7
Low-rank/Chebyshev	(Ruiz-Antolin et al., 2017)	$$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$8
HSS direct solver	(Li et al., 3 Dec 2025)	Setup: $$F(\omega_k) = \sum_{j=1}^M f(x_j) e^{-2\pi i \langle x_j,\,\omega_k\rangle},\qquad k=1,\ldots,K.$$9,<br> Solve: $f_j = \sum_{k=1}^{N} u_k\,e^{2\pi i x_j \omega_k}$0

Each approach supports high-precision transforms for arbitrary nonuniform grids, with different balances between setup time, memory, and runtime, depending on application context and underlying data structure.

---

References:

• (Li et al., 3 Dec 2025) Li & Liu, "A Superfast Direct Solver for Type-III Inverse Nonuniform Discrete Fourier Transform"
• (Sultan et al., 9 Jan 2025) Sultan et al., "Optimized Sampling for Non-Line-of-Sight Imaging Using Modified Fast Fourier Transforms"
• (Barnett et al., 2018) Barnett et al., "A parallel non-uniform fast Fourier transform library based on an 'exponential of semicircle' kernel"
• (Ruiz-Antolin et al., 2017) Ruiz-Antolín & Townsend, "A nonuniform fast Fourier transform based on low rank approximation"