Efficient multivariate inversion of the NFFT

Various applications such as MRI, solution of PDEs, etc. need to perform an inverse nonequispaced fast Fourier transform (iNFFT), i. e., compute Fourier coefficients from given nonequispaced data. In contrast to iterative solvers we introduce a direct method for this inversion in the overdetermined setting. For this purpose, we use the matrix representation of the NFFT. Modifying one of the matrix factors leads to an optimization problem, which can be solved in a precomputational step using normal equations. Therefore, we are able to compute an inverse NFFT up to a certain accuracy by means of a modified adjoint NFFT and therefore preserve its arithmetic complexity.


Nonequispaced fast Fourier transform
The nonequispaced fast Fourier transform (NFFT), cf. [1][2][3], is a fast algorithm to evaluate a trigonometric polynomial this problem can be written as f = Af . The fast algorithm can also be formulated in matrix-vector notation and we receive the approximation A ≈ BF D, where D is a diagonal matrix, F is a truncated d-variate Fourier matrix and B is sparse.
The adjoint problem can be written as h = A * f with the adjoint matrix A * . Hence, a fast algorithm for this problem follows by using the adjoint matrix-vector factorization A * ≈ D * F * B * .

Inversion of the NFFT
Similar to [4] we face the following two problems.
(2) Solve This will be solved by an inverse adjoint NFFT.
Considering (2.1) we seek to find a matrix X such that XA ≈ I |I M | , i. e., the identity of size |I M |. From the equispaced case we know A * A = N I |I M | for |I M | ≤ N . Therefore, we aim to modify the matrix B such that D * F * B * A ≈ I |I M | holds. For this purpose, the entries of the matrix B should be calculated such that its sparse structure with at most (2m + 1) d entries per row and consequently the arithmetic complexity of the algorithms is preserved. A matrix B satisfying this property we call (2m + 1 ) d -sparse. In other words, we consider the optimization problem where M σ = σM with the d-dimensional oversampling factor σ ≥ 1 and the Frobenius norm is denoted by · F . Suppose A * B F D ≈ I |I M | , i. e., A * B is a pseudoinverse of F D. Since F * F = |I M σ | I |I M | , a pseudoinverse is given by Here we already rewrote the norm based on the definition of the Frobenius norm in terms of the nonzeros of the columns ofB, namedb l ∈ R |I M σ ,m (l)| , where we defined the index set Additionally, we set H l ∈ C |I M |×|I M σ ,m (l)| as the appropriate cutout of A * and f l are the columns of F * . Thus, if the matrix H l has full rank the solution is given by b opt Having these vectors b opt l we compose the modified matrix B opt , observing that b opt l only consist of the nonzero entries of B opt . Then the approximation of the Fourier coefficients is given byf ≈f = D * F * B * opt f . In other words, this approach yields an inverse NFFT by modifying the adjoint NFFT.
Remark 2.1 Considering (2.2) we recognize that this problem can also be solved by achieving an approximation of the form A * BF D ≈ I |I M | . In other words, our approach provides an inverse NFFT as well as an inverse adjoint NFFT.

Numerical Example
Example 3.1 We have a look at the reconstruction of the Shepp-Logan phantom, cf. Figure 1 (a). For the phantom P (M ) of size M we setf = vec(P (M )) and for the reconstructed coefficientsf we setP := vec(f ) −1 .
The reconstruction of the phantom of size 1024 × 1024 is presented in Figure 1 including a detailed view of the 832nd row of this reconstruction. For the purpose of comparison we also added the result of reconstructing using the original adjoint NFFT. It can be seen that by our optimization procedure we achieve a very good reconstruction.