Nonuniform fast Fourier transforms

Specific Domains Numerics
1

FastTransforms.jl v0.2.2 includes nonuniform fast Fourier transforms thanks to Alex Townsend:

  • nufft1 assumes uniform samples and noninteger frequencies;
  • nufft2 assumes nonuniform samples and integer frequencies;
  • nufft3 ( = nufft) assumes nonuniform samples and noninteger frequencies;
  • inufft1 inverts an nufft1; and,
  • inufft2 inverts an nufft2.

Basic usage is like so:

julia> Pkg.add("FastTransforms")

julia> using FastTransforms

julia> n = 10^4;

julia> c = complex(rand(n));

julia> ω = collect(0:n-1) + rand(n);

julia> nufft1(c, ω, eps());

julia> p1 = plan_nufft1(ω, eps());

julia> @time p1*c;
  0.002383 seconds (6 allocations: 156.484 KiB)

The new algorithms are described in the following pre-print:

D. Ruiz—Antolín and A. Townsend. A nonuniform fast Fourier transform based on low rank approximation, arXiv:1701.04492, 2017.

7 Likes
2

@MikaelSlevinsky: Are you aware of https://github.com/tknopp/NFFT.jl? It seems that nufft just allows for 1D NFFTs whereas NFFT.jl allows for d-dimensional transforms.

I am pretty interested in a performance comparison of both approaches.

3

And from an interface point of view: Isn’t nufft2 the adjoint of nufft2 (modulo conjugation). The NFFT is quite often used as an operator when solving linear systems of equations and there you usually want that pair. In my (internal) software on MRI reconstruction we overload A_mul_B! and Ac_mul_B! for this purpose.

4

Hi Tobias, there is some support for 2D nufft’s of types I-I and II-II, but the multi-dimensional support is not fully fleshed out. You are right, nufft1 and nufft2 are adjoints, modulo conjugation.

I’m interested in a performance comparison, too, and the first step is actually getting a release of this approach.

There are four levels to the 1D interface:

  • the functions nufft1/2/3(args...) plan & execute;
  • the functions plan_nufft1/2/3(args...) construct plans;
  • if p is an NUFFTPlan, then p*c applies the plan; and,
  • A_mul_B!(out, p, in) applies the plan to the in array in-place on the out array.

For a fair comparison, methods should be compared at the same tolerance ϵ.

Note that a hierarchical approach to the Chebyshev–Legendre transform accelerated the execution phase (of any previous methods in FastTransforms.jl) by a factor of ~20. With some effort, Dutt & Rokhlin’s hierarchical nonuniform FFT method could be implemented based on the HierarchicalMatrices.jl framework.

1 Like
5

I just saw that the paper actually makes a comparison with the NFFT.jl, which is great.

I wonder if it would be good to share an interface. NFFT.jl could have different implementations that could be switchable.

6

Yes, in Figure 2.2 we did a quick comparison to NFFT.jl for the 1D NUFFT of type-I:

1 Like
7

Maybe a meta-package a la DifferentialEquations.jl is in order? Transforms.jl is available it looks like…there could even be an org JuliaTransforms.

8

No, in this case I do not really see a need for that. It should be possible to have a single NFFT package that allows to switch different algorithms in its backend.
FastTransforms.jl is more high level and can then integrate the low level packages (where NFFT is just one example)

9

I think that’s what I meant: a package that wraps the different algorithms for the different transforms, with a consistent naming scheme and consistent way to specify algorithms.

10

Talking about a common interface maybe this package could be interesting: AbstractOperators.jl.

This offers a common interface that allows to easily combine linear and nonlinear transformations and use them in inverse problems (through the package StructuredOptimization).

In case you find this interesting feel free to contribute!

11

FYI: There’s also GitHub - ludvigak/FINUFFT.jl: Julia interface to the nonuniform FFT library FINUFFT that wasn’t available during the above discussion (a wrapper for flatironinstitute’s C++ code, or I guess its C interface), nor the paper on it:

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. […] 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.

Also interesting: GitHub - gwater/HexFFT.jl: Fast Fourier transform on hexagonal grids using Birdsong and Rummelt's algorithm “implements Birdsong and Rummelt’s 2016 algorithm for fast Fourier transforms on hexagonal lattices.”

The current implementation is ca. 10 times slower than Base.fft() for comparable rectangular grids. However significant optimization should be possible by preallocating and reusing temporary arrays.