Fast and direct inversion methods for the multivariate nonequispaced fast Fourier transform

The well-known discrete Fourier transform (DFT) can easily be generalized to arbitrary nodes in the spatial domain. The fast procedure for this generalization is referred to as nonequispaced fast Fourier transform (NFFT). Various applications such as MRI and solution of PDEs are interested in the inverse problem, i.e., computing Fourier coefficients from given nonequispaced data. In this article, we survey different kinds of approaches to tackle this problem. In contrast to iterative procedures, where multiple iteration steps are needed for computing a solution, we focus especially on so-called direct inversion methods. We review density compensation techniques and introduce a new scheme that leads to an exact reconstruction for trigonometric polynomials. In addition, we consider a matrix optimization approach using Frobenius norm minimization to obtain an inverse NFFT.


Introduction
The NFFT, short hand for nonequispaced fast Fourier transform or nonuniform fast Fourier transform (NUFFT), respectively, is a fast algorithm to evaluate a trigonometric polynomial f (x) = given equispaced points x j and |I M | = N , this evaluation can be realized by means of the well-known fast Fourier transform (FFT); an algorithm that is invertible. However, various applications such as magnetic resonance imaging (MRI), cf. [15,19], solution of PDEs, cf. [21], etc., need to perform an inverse nonequispaced fast Fourier transform (iNFFT), i. e., compute the Fourier coefficientsf k from given function evaluations f (x j ) of the trigonometric polynomial (1.1). Hence, we are interested in an inversion also for nonequispaced data.
In general, the number N of points x j is independent from the number |I M | of Fourier coefficientsf k and hence the nonequispaced Fourier matrix which we would have to invert, is rectangular in most cases. Considering the corresponding linear system Af = f with f := (f (x j )) N j=1 andf := (f k ) k∈I M , this can either be overdetermined, if |I M | ≤ N , or underdetermined, if |I M | > N . Generally, this forces us to look for a pseudoinverse solution. Moreover, we also require that the nonequispaced Fourier matrix A has full rank. Eigenvalue estimates in [22,6,9,45,43] indeed confirm that this condition is satisfied for sufficiently nice sampling sets.
In literature a variety of approaches for an inverse NFFT (iNFFT) can be found. This is why we give a short overview.
Iterative inversion methods We start surveying the iterative methods. For the onedimensional setting d = 1 with |I M | = N an algorithm was published in [60], which is specially designed for jittered equispaced points and is based on the conjugate gradient (CG) method in connection with low rank approximation. An approach for the overdetermined case |I M | ≤ N can be found in [22], where the Toeplitz structure of the matrix product A * W A with a diagonal matrix W := diag(w j ) N j=1 of Voronoi weights is used to compute the solution iteratively by means of the CG algorithm.
For higher dimensional problems with d ≥ 1 there are several approaches that compute a least squares approximation to the linear system Af = f . In the overdetermined case |I M | ≤ N , the given data can typically only be approximated up to a residual r := Af − f . Therefore, the weighted least squares problem is considered, which is equivalent to solving the weighted normal equations of first kind A * W Af = A * W f with the diagonal matrix W := diag(w j ) N j=1 of weights in time domain. In [69,24,41] these normal equations are solved iteratively by means of CG using the NFFT to realize fast matrix-vector multiplications involving A, whereas in [56] a fast convolution is used.
In the consistent underdetermined case |I M | > N the data can be interpolated exactly and therefore one can choose a specific solution, e. g. the one that solves the constrained minimization problem This interpolation problem is equivalent to the weighted normal equations of second kind AŴ A * y = f ,f =Ŵ A * y with the diagonal matrixŴ := diag(ŵ k ) k∈I M of weights in frequency domain. In [44] the CG method was used in connection with the NFFT to iteratively compute a solution to this problem, see also [53,Section 7.6.2].
Regularization methods Moreover, there also exist several regularization techniques for the multidimensional setting d ≥ 1. For example, [65,71, 1] all solve the 1 -regularized problem with regularization parameter λ > 0 and the m-th order polynomial annihilation operator L m ∈ R N ×|I M | as sparsifying transform, see [2]. Based on this, weighted p -schemes were introduced in [10,11,14,48], which are designed to reduce the penalty at locations where L mf is nonzero. For instance, [13,61] each state a two step method, that firstly uses edge detection to create a mask, i. e., a weighting matrix which indicates where non-zero entries are expected in the TV domain, and then targets weighted 2 -norm TV regularization appropriately to smooth regions of the function in a second minimization step.
Direct inversion methods In contrast to these iterated procedures, there are also socalled direct methods that do not require multiple iteration steps. Already in [18] a direct method was explained for the setting d = 1 and |I M | = N which uses Lagrange interpolation in combination with fast multiple methods. Based on this, further methods were deduced for the same setting, which also use Lagrange interpolation, but additionally incorporate an imaginary shift in [63], or utilize fast summation in [37] for the fast evaluation of occurring sums, see also [53,Section 7.6.1].
In the overdetermined setting |I M | ≤ N another approach for computing an inverse NFFT can be obtained by using the fact that A * A is of Toeplitz structure. To this end, the Gohberg-Semencul formula, see [32], can be used to solve the normal equations A * Af = A * f exactly by analogy with [5]. Here the computation of the components of the Gohberg-Semencul formula can be viewed as a precomputational step. In addition, also a frame-theoretical approach is known from [26], which provides a link between the adjoint NFFT and frame approximation, and could therefore be seen as a way to invert the NFFT. Note that the method in [26] is based only on optimizing a diagonal matrix (the matrix D defined in (2.13)), whereas in [37] similar ideas were used to modify a sparse matrix (the matrix B defined in (2.15)).
For the multidimensional setting d > 1 several methods have been developed that are tailored to the special structure of the linogram or pseudo-polar grid, respectively, see Figure 5.2b, such that the inversion involves only one-dimensional FFTs and interpolations. On the one hand, in [4] a least squares solution is computed iteratively by using a fast multiplication technique with the matrix A, which can be derived in case of the linogram grid. On the other hand, [3] is based on a fast resampling strategy, where a first step resamples the linogram grid onto a Cartesian grid, and the second phase recovers the image from these Cartesian samples. However, these techniques are exclusively applicable for the special case of the linogram grid, see Figure 5.2b, or the polar grid by another resampling, cf. Figure 5.2a. Since we are interested in more generally applicable methods, a brief introduction to direct inversion for general sampling patterns can be found in [38].
Current approach In this paper we focus on direct inversion methods that are applicable for general sampling patterns and present new methods for the computation of an iNFFT. Note that direct method in the context of the linear system Af = f means, that for a fixed set of points x j , j = 1, . . . , N , the reconstruction off from given f can be realized with the same number of arithmetic operations as a single application of an adjoint NFFT (see Algorithm 2.4). To achieve this, a certain precomputational step is compulsory, since the adjoint NFFT does not yield an inversion of the NFFT per se, cf. (3.3). Although this precomputations might be rather costly, they need to be done only once for a given set of points x j , j = 1, . . . , N . Therefore, direct methods are especially beneficial in case of fixed points for several measurement vectors f .
For this reason, the current paper is concerned with two different approaches of this kind. Firstly, we consider the very well known approach of so-called sampling density compensation, which can be written asf ≈ A * W f with a diagonal matrix W := diag(w j ) N j=1 of weights. The already mentioned precomputations then consist of computing suitable weights w j , while the actual reconstruction step includes only one adjoint NFFT applied to the scaled coefficient vector W f . In this paper we examine several existing approaches for computing the weights w j and introduce a new method, such that the iNFFT is exact for all trigonometric polynomials (1.1) of degree M .
As a second part, we reconsider and enhance our approach introduced in [38]. Here the idea is using the matrix representation A ≈ BF D of the NFFT to modify one of the matrix factors, such that an inversion is given byf ≈ D * F * B * opt f . Then the precomputational step consists of optimizing the matrix B, while the actual reconstruction step includes only one modified adjoint NFFT applied to the coefficient vector f .
Outline of this paper This paper is organized as follows. In Section 2 we introduce the already mentioned algorithm, the NFFT, as well as its adjoint version, the adjoint NFFT. Secondly, in Section 3 we introduce the inversion problem and deal with direct methods using so-called sampling density compensation. We start our investigations with trigonometric polynomials in Section 3.1. Here the main formula (3.14) yields exact reconstruction for all trigonometric polynomials of degree M . We also discuss practical computation schemes for the overdetermined as well as the underdetermined setting. Subsequently, in Section 3.2 we go on to bandlimited functions and show that the same numerical procedures as in Section 3.1 can be used in this setting as well. Section 3.3 then summarizes the previous findings by presenting a general error bound on density compensation factors computed by means of (3.14) in Theorem 3.14. In addition, this also yields an estimate on the condition number of a specific matrix product, as shown in Theorem 3.15. In Section 3.4 we have a look at certain approaches from literature and their connection among each other as well as to the method presented in Section 3.1. Afterwards, we examine another direct inversion method in Section 4, where we aim to modify the adjoint NFFT by means of matrix optimization such that this yields an iNFFT. Finally, in Section 5 we show some numerical examples to investigate the accuracy of our approaches.

Nonequispaced fast Fourier transform
The inner product of two vectors shall be defined as usual as kx := k 1 x 1 + · · · + k d x d . Additionally, we define the componentwise product as x y := (x 1 y 1 , . . . , x d y d ) T , the all ones vector 1 d := (1, . . . , 1) T and the reciprocal of a vector x with nonzero components shall be given by We consider the Hilbert space L 2 (T d ) of all 1-periodic, complex-valued functions, which possesses the orthonormal basis {e 2πikx : k ∈ Z d }. It it known that every function f ∈ L 2 (T d ) is uniquely representable in the form with the coefficients which is an acceptable approximation for k ∈ I M , see e. g. [53, p. 214]. The fast evaluation of (2.3) can then be realized by means of the fast Fourier transform (FFT). Moreover, it is know that this transformation is invertible and that the inverse problem of computing can be realized by means of an inverse fast Fourier transform (iFFT), which is basically the same algorithm except for some reordering and scaling, cf. [53,Lem. 3.17]. Now suppose we are given nonequispaced points x j ∈ T d , j = 1, . . . , N , instead. Then, we consider the computation of the sums for givenf k ∈ C, k ∈ I M , as well as the adjoint problem of computing for given f j ∈ C, j = 1, . . . , N . By defining the nonequispaced Fourier matrix as well as the vectors f := (f j ) N j=1 ,f := (f k ) k∈I M and h := (h k ) k∈I M , the computation of sums of the form (2.5) and (2.6) can be written as f = Af and h = A * f , where A * := A T denotes the adjoint matrix of A.

The NFFT
We firstly restrict our attention to problem (2.5), which is equivalent to the evaluation of a trigonometric polynomial with givenf k ∈ C, k ∈ I M , at given nonequispaced points x j ∈ T d , j = 1, . . . , N . Let w ∈ L 2 (R d ) ∩ L 1 (R d ) be a so-called window function, which is well localized in space and frequency domain. Now we define the 1-periodic functionw(x) := r∈Z d w(x + r) with absolute convergent Fourier series. As a consequence, the Fourier coefficients of the periodizationw have the form For a given oversampling factor σ ≥ 1, we define 2N M σ := 2 σM /2 as well as M σ := M σ · 1 d , and approximate f by a linear combination of translates of the periodized window function, i. e., where g ∈ C, ∈ I M σ , are coefficients to be determined such that (2.9) yields a good approximation. By means of the convolution theorem (see [53,Lem. 4.1]), the approximant s 1 ∈ L 2 (T d ) in (2.9) can be represented as where the discrete Fourier transform of the coefficients g is defined bŷ Comparing (2.5) and (2.10) then yieldŝ Consequently, the coefficients g in (2.9) can be obtained by inverting (2.11), i. e., by the application of an iFFT. Furthermore, we assume that w is well localized such that it is small outside the square [− m /Mσ, m /Mσ] d with truncation parameter m M σ . In this case, w can be approximated by the compactly supported function : otherwise.
Thereby, we approximate s 1 by the short sums where the index set 12) contains at most (2m + 1) d entries for each fixed x j . Thus, the obtained algorithm can be summarized as follows.

Setĝ
Next we give the matrix-vector representation of the NFFT. To this end, we define • the diagonal matrix

13)
• the truncated Fourier matrix 14) • and the sparse matrix where by definition (2.12) each row of B contains at most (2m + 1) d nonzeros. In doing so, the NFFT in Algorithm 2.1 can be formulated in matrix-vector notation such that we receive the approximation A ≈ BF D of (2.7), cf. [53, p. 383]. In other words, the NFFT uses the approximation 3. It has to be pointed out that because of consistency the factor |I M σ | −1 is here not located in the matrix F as usual but in the matrix D.

The adjoint NFFT
Now we proceed with the adjoint problem (2.6). As already seen, this can be written as h = A * f with the adjoint matrix A * of (2.7). Thus, using the matrices (2.13), (2.14) and (2.15) we receive the approximation A * ≈ D * F * B * , such that a fast algorithm for the adjoint problem can be denoted as follows.

Algorithm 2.4 (adjoint NFFT).
For d, N ∈ N let x j ∈ T d , j = 1, . . . , N , be given points as well as f j ∈ C given coefficients. Further, we are given the oversampling factor σ ≥ 1, 2N M σ := 2 σM /2 , M σ := M σ · 1 d , as well as the window function w, the truncated function w m with truncation parameter m M σ , and their 1-periodic versionsw andw m .
1. Compute the sparse sums The algorithms presented in this section (Algorithms 2.1 and 2.4) are part of the NFFT software [35]. For algorithmic details we refer to [36].

Direct inversion using density compensation
Having introduced the fast methods for nonequispaced data, we remark that various applications such as MRI, solution of PDEs, etc. are interested in the inverse problem, i. e., instead of the evaluation of (2.5) the aim is computing the Fourier coefficientsf k , k ∈ I M , from given nonequispaced data f (x j ), j = 1, . . . , N . Therefore, this section shall be dedicated to this task.
To clarify the major dissimilarity between equispaced and nonequispaced data, we start considering the equispaced case. When evaluating at the points x j = 1 n j ∈ T d , j ∈ I n , with n := n · 1 d and |I n | = N , the nonequispaced Fourier matrix A ∈ C N ×|I M | in (2.7) turns into the equispaced Fourier matrix F ∈ C |I M σ |×|I M | from (2.14) with |I M σ | = N . Thereby, it results from the geometric sum formula that as well as and |I M | is divisible by N . Thus, in the equispaced setting a one-sided inverse is given by the (scaled) adjoint matrix. However, when considering arbitrary points x j ∈ T d , j = 1, . . . , N , this property is lost, i. e., for the nonequispaced Fourier matrix we have Because of this, more effort is needed in the nonequispaced setting.
In general, we face the following two problems.
(1) Solve the linear system i. e., reconstruct the Fourier coefficientsf = (f k ) k∈I M from given function values f = (f (x j )) N j=1 . This problem is referred to as inverse NDFT (iNDFT) and an efficient solver shall be called inverse NFFT (iNFFT).
(2) Solve the linear system (3.5) i. e., reconstruct the coefficients f = (f j ) N j=1 from given data h = (h k ) k∈I M . This problem is referred to as inverse adjoint NDFT (iNDFT*) and an efficient solver shall be called inverse adjoint NFFT (iNFFT*).
Note that in both problems the numbers |I M | and N are independent, such that the nonequispaced Fourier matrix A ∈ C N ×|I M | in (2.7) is generally rectangular.
At first, we restrict our attention to problem (3.4). When considering iterative inversion procedures as those mentioned in the introduction, these methods require multiple iteration steps by definition. Therefore, multiple matrix vector multiplications with the system matrix A, or rather multiple applications of the NFFT (see Algorithm 2.1), are needed to compute a solution. To reduce the computational effort, we now proceed, in contrast to this iterated procedures, with so-called direct methods. In the setting of problem (3.4) we hereby mean methods, where for a fixed set of points x j , j = 1, . . . , N , the reconstruction off from given f can be realized with the same number of arithmetic operations as a single application of an adjoint NFFT (see Algorithm 2.4). To achieve this, a certain precomputational step is compulsory, since the adjoint NFFT does not yield an inversion of the NFFT per se, see (3.3). Although this precomputations might be rather costly, they need to be done only once for a given set of points x j , j = 1, . . . , N . In fact, the actual reconstruction step is very efficient. Therefore, direct methods are especially beneficial in case we are given fixed points for several measurement vectors f .
In this section we focus on a direct inversion method for solving problem (3.4) that utilizes so-called sampling density compensation. To this end, we consider the integral (2.2) and introduce a corresponding quadrature formula. In contrast to the already known equispaced approximation (2.3) we now assume given arbitrary, nonequispaced points x j ∈ T d , j = 1, . . . , N . Thereby, the Fourier coefficients (2.2) are approximated by a general quadrature rule using quadrature weights w j ∈ C, j = 1, . . . , N , which are needed for sampling density compensation due to the nonequispaced sampling. Thus, for a trigonometric polynomial (2.8) we havê Using the nonequispaced Fourier matrix A ∈ C N ×|I M | in (2.7), the diagonal matrix of weights W := diag(w j ) N j=1 ∈ C N ×N as well as the vector h w := (h w k ) k∈I M , the nonequispaced quadrature rule (3.6) can be written asf ≈ h w := A * W f . For achieving a fast computation method we make use of the approximation of the adjoint NFFT, cf. Section 2.2, i. e., the final approximation is given bŷ with the matrices D ∈ C |I M |×|I M | , F ∈ C |I M σ |×|I M | and B ∈ R N ×|I M σ | defined in (2.13), (2.14) and (2.15). In other words, for density compensation methods the already mentioned precomputations consist of computing the quadrature weights w j ∈ C, j = 1, . . . , N , while the actual reconstruction step includes only one adjoint NFFT (see Algorithm 2.4) applied to the scaled measurement vector W f . The aim of all density compensation techniques is then to choose appropriate weights w j ∈ C, j = 1, . . . , N , such that the underlying quadrature (3.6) is preferably exact. In the following we have a look at the specific choice of the so-called density compensation factors w j .
An intuitive approach for density compensation is based on geometry, where each sample is considered as representative of a certain surrounding area, as in numerical integration. The weights for each sample can be obtained for instance by constructing a Voronoi diagram and calculating the area of each cell, see e. g. [58]. This approach of Voronoi weights is well-known and widely used in practice. However, it does not necessarily yield a good approximation (3.7), which is why we examine some more sophisticated approaches in the remainder of this section.
To this end, this section is organized as follows. Firstly, in Section 3.1 we introduce density compensation factors w j , j = 1, . . . , N, that lead to an exact reconstruction formula (3.6) for all trigonometric polynomials (2.8) of degree M . In addition to the theoretical results, we also discuss methods for the numerical computation. Secondly, in Section 3.2 we show that it is reasonable to consider the inversion problem (3.4) and density compensation via Subsequently, we summarize our previous findings by presenting a general error bound on density compensation factors in Section 3.3. Finally, in Section 3.4 we reconsider certain approaches from literature and illustrate their connection among each other as well as to the method introduced in Section 3.1. (i) If we define g := W f , i. e., each entry of f is scaled with respect to the points x j , j = 1, . . . , N , the approximation (3.7) can be written asf ≈ D * F * B * g. As mentioned before, this coincides with an ordinary adjoint NFFT applied to a modified coefficient vector g.
(ii) By defining the matrixB := W * B, i. e., scaling the rows of B with respect to the points x j , j = 1, . . . , N , the approximation (3.7) can be written asf ≈ D * F * B * f . In this sense, density compensation can also be seen as a modification of the adjoint NFFT and its application to the original coefficient vector.
Note that (i) is the common viewpoint. However, we keep (ii) in mind, since this allows treating density compensation methods as an optimization of the sparse matrix B ∈ R N ×|I M σ | in (2.15), as it shall be done in Section 4. We remark that density compensation methods allow only N degrees of freedom.

Exact quadrature weights for trigonometric polynomials
Similar to [28], we aim to introduce density compensation factors w j , j = 1, . . . , N, that lead to an exact reconstruction formula (3.6) for all trigonometric polynomials (2.8) of degree M . To this end, we firstly examine certain properties that arise from (3.6) being exact.
. . , N, and quadrature weights w j ∈ C be given. Then an exact reconstruction formula (3.6) for trigonometric polynomials (2.8) with maximum degree M satisfyinĝ implies the following equivalent statements.
(i) The quadrature rule is exact for all trigonometric polynomials (2.8) with maximum degree M .
(ii) The linear system of equations (3.8) ⇒ (i): By inserting the definition (2.8) of a trigonometric polynomial of degree M into the integral considered in (3.9) we have with the Kronecker delta δ 0,k . Now using the property (3.8) as well as definition (3.6) of h w k we proceed witĥ such that (3.11) combined with (3.12) yields the assertion (3.9). (i) ⇒ (ii): Inserting the definition (2.8) of a trigonometric polynomial of degree M into the right-hand side of (3.9) implies (3.13) This together with the property (i) and (3.11) leads tô w j e 2πikx j and thus to assertion (3.10).
(ii) ⇒ (i): Combining (3.11), (3.10) and (3.13) yields the assertion via Comparable results can also be found in literature. A fundamental theorem in numerical integration, see [70], states that for any integral T d f (x) dx there exists an exact quadrature rule (3.9), i. e., optimal points x j ∈ T d and weights w j ∈ C, j = 1, . . . , N , such that (3.9) is fulfilled. In [31, Lemma 2.6] it was shown that for given points x j ∈ T d , j = 1, . . . , N, certain quadrature weights w j can be stated by means of frame theoretical considerations which lead to an exact quadrature rule (3.9) by definition. Moreover, it was shown (cf. [31,Lemma 3.6]) that these weights are the ones with minimal (weighted) 2 -norm, which are already known under the name "least squares quadrature", see [34]. According to [34,Sec. 2.1] these quadrature weights w j , j = 1, . . . , N, can be found by solving a linear system of equations dx for a given set of basis functions {φ k } k∈I M . In our setting we have φ k (x) = e 2πikx and therefore i. e., the same linear system of equations as in (3.10). We remark that both [31] and [34] state the results in the case d = 1, a generalization to d > 1, however, is straight-forward.
By means of Theorem 3.2 we can now give a condition that guaranties (3.6) being exact for all trigonometric polynomials (2.8) with maximum degree M . However, an augmented variant of (3.10), namely

14)
yields an exact reconstructionf k = h w k in (3.6) for trigonometric polynomials (2.8) with maximum degree M . Additionally, (3.14) implies the matrix equation Proof. Utilizing definitions (3.6) and (2.8) we have Since (3.10) only holds for k, ∈ I M with ( − k) ∈ I M , this implies where for all k ∈ I M \ {0} there exists an ∈ I M with ( − k) ∈ I 2M \ I M .
As ( − k) ∈ I 2M for k, ∈ I M , the augmented variant (3.14) yields Moreover, since δ ( −k),0 = δ k, , the condition (3.14) implies In matrix-vector notation this can be written as We remark that this matrix equation immediately shows that we have an exact reconstruction of the form (3.8), since if be an arbitrary 1-periodic function (2.1). Then (3.14) yields i. e., for a function f ∈ L 2 (T d ) we only have a good approximation in case the coefficients c (f ) are small for / ∈ I M , whereas this reconstruction can only be exact for f being a trigonometric polynomial (2.8).
3.1.1 Practical computation in the underdetermined setting |I 2M | ≤ N So far, we have seen in Corollary 3.4 that an exact solution w = (w j ) N j=1 to the linear system (3.14) leads to an exact reconstruction formula (3.6) for all trigonometric polynomials (2.8) with maximum degree M . Therefore, we aim to use this condition (3.14) to numerically find optimal density compensation factors w j ∈ C, j = 1, . . . , N .
Having a closer look at the condition (3.14) we recognize that it can be written as the , and right side e 0 := (δ 0,k ) k∈I 2M . We remark that in contrast to A ∈ C N ×|I M | we now deal with the enlarged matrix A |I 2M | ∈ C N ×|I 2M | , such that single matrix operations are more costly. Nevertheless, Corollary 3.4 yields a direct inversion method for (3.4), where the system A T |I 2M | w = e 0 needs to be solved only once for fixed points x j ∈ T d , j = 1, . . . , N . Its solution w can then be used to efficiently approximatef for multiple measurement vectors f , whereas iterative methods for (3.4) need to solve Af = f each time.
As already mentioned in [34,Sec. 3.1] an exact solution to (3.14) can only be found if |I 2M | ≤ N , i. e., in case A T |I 2M | w = e 0 is an underdetermined system of equations. By [34,Lem. 3.1] this system has at least one solution, which is why we may choose the one with minimal 2 -norm. If rank(A |I 2M | ) = |I 2M |, then the system A T |I 2M | w = e 0 is consistent and the unique solution is given by the normal equations of second kind More precisely, we may compute the vector v using an iterative procedure such as the CG algorithm, such that only matrix multiplications with A T |I 2M | and A |I 2M | are needed. Since fast multiplication with A T |I 2M | and A |I 2M | can easily be realized by means of an adjoint NFFT (see Algorithm 2.4) and an NFFT (see Algorithm 2.1), respectively, computing the solution w to (3.15 Thus, in order to receive exact quadrature weights w j ∈ C, j = 1, . . . , N, via (3.15) we need to satisfy the full rank condition rank(A |I 2M | ) = |I 2M |. In case of a low rank matrix A |I 2M | ∈ C N ×|I 2M | for |I 2M | ≤ N , we may still use (3.15) to obtain a least squares approximation to (3.14).

Practical computation in the overdetermined setting |I 2M | > N
In the setting |I 2M | > N , we cannot expect finding an exact solution w to (3.14), since we have to deal with an overdetermined system possessing more conditions than variables. However, we still aim to numerically find optimal density compensation factors w j ∈ C, j = 1, . . . , N , by considering a least squares approximation to (3.14) that mini- , and e 0 = (δ 0,k ) k∈I 2M we simplify the right hand side via Since fast multiplication with A T |I 2M | and A |I 2M | can easily be realized by means of an adjoint NFFT (see Algorithm 2.4) and an NFFT (see Algorithm 2.1), respectively, the solution w to (3.16) can be computed iteratively by means of the CG algorithm in O(|I 2M | log(|I 2M |) + N ) arithmetic operations. Note that the solution to (3.16) is only unique if the full rank condition rank(A |I 2M | ) = N is satisfied, cf. [8, p. 7]. We remark that the computed weight matrix W = diag(w) can further be used in an iterative procedure as in [53,Alg. 7.27] to improve the approximation off .
The previous considerations lead to the following algorithms.
Algorithm 3.6 (Computation of the optimal density compensation factors).
For d, N ∈ N let x j ∈ T d , j = 1, . . . , N , as well as M ∈ 2N and M := M · 1 d be given.
Compute the solution w to (3.16) iteratively using the NFFT.
0. Precompute the weights matrix W using Algorithm 3.6.

Computeh
, by means of an adjoint NFFT.

Bandlimited functions
In some numerical examples, such as in MRI, we are concerned with bandlimited func- cf. [15]. In the following we show that it is reasonable to consider the inversion problem (3.4) as well as density compensation via (3.7) for bandlimited functions as well.
To this end, let f ∈ L 1 (R d ) ∩ C 0 (R d ) be a bandlimited function with bandwidth M , i. e., its (continuous) Fourier transform Analogous to (2.3), the approximation using equispaced quadrature points k ∈ I M yields 19) such that evaluation at the given nonequispaced points By means of the definition (2.7) of the matrix A ∈ C N ×|I M | this can be written as f ≈ Af , where we used the notationf := (f (k)) k∈I M in this setting. Thus, also for bandlimited functions f its evaluations at points x j can be approximated in the form (2.5), such that it is reasonable to consider the inversion problem (3.4) for bandlimited functions as well. Considering (3.17) we are given an exact formula for the evaluation of the Fourier transformf . However, in practical applications, such as MRI, this is only a hypothetical case, since f cannot be sampled on whole R d , cf. [15]. Due to a limited coverage of space by the acquisition, the function f is typically only known on a bounded domain, w.l.o.g. for x ∈ − 1 2 , 1 2 d . Thus, we have to deal with the approximation Using the nonequispaced quadrature rule in (3.6), we find that evaluation at uniform grid points k ∈ I M can be approximated viâ This is to say, equispaced samples of the Fourier transform of a bandlimited function may be approximated in the same form (f (k)) k∈I M =f ≈ h w := A * W f as in (3.6), where we used the notation h w := (h(k)) k∈I M in this setting. Moreover, we extend this approximation onto the whole interval − M 2 , M 2 d , i. e., we consider So all in all, we have seen that it is reasonable to study the inversion problem (3.4) and the associated density compensation via (3.7) for bandlimited functions as well. Analogous to Section 3.1 we now aim to find a numerical method for computing suitable weights w j ∈ C, j = 1, . . . , N , such that the reconstruction formula (3.21) is preferably exact. To this end, we have a closer look at (3.21) being exact and start with analogous considerations as in Theorem 3.2.
as quadrature weights w j ∈ C, j = 1, . . . , N, be given. Then an exact reconstruction formula implies that the quadrature rule Proof. By (3.17) the assumption (3.22) can be written as Especially, for v = 0 evaluation of (3.23) yields the assertion However, in contrast to Theorem 3.2, this Theorem 3.8 does not yield an explicit condition for computing suitable weights w j ∈ C, j = 1, . . . , N . To derive a numerical procedure anyway, we generalize the notion of an exact reconstructionh of f and have a look at the theory of tempered distributions. To this end, let S (R d ) be the Schwartz space of rapidly decaying functions, cf. [53, Sec. 4.2.1]. The tempered Dirac distribution δ shall be defined by δ, ϕ : [53,Ex. 4.36]. For a slowly increasing function f : R d → C satisfying |f (x)| ≤ c(1 + x 2 ) n almost everywhere with c > 0 and n ∈ N 0 , the induced distribution T f shall be defined by T f , ϕ : Then the following property can be shown. weights w j ∈ C be given. Further let T f be the distribution induced by some bandlimited
Proof. Using the definition of the functionh in (3.21) as well as the fact that the inversion formula (3.18) holds for all x ∈ R d , we have w j e 2πivx j dv du.
Hence, by (3.24) this implies Considering the property (3.26), we remark that this indeed states an exact reconstructionf =h in the sense of tempered distributions, as distinct from (3.22). Since it is known by Corollary 3.4 that the condition (3.14) yields an exact reconstruction for trigonometric polynomials, we aim to use this result to compute suitable weights w j ∈ C, j = 1, . . . , N , for bandlimited functions as well. To this end, suppose we have (3.14), i. e., N j=1 w j e 2πikx j = δ 0,k , k ∈ I 2M .
Then this yields (3.27) Having a look at Theorem 3.9, an exact reconstruction (3.26) is implied by (3.24), i. e., Thus, the property (3.27) that is fulfilled by (3.14) could be interpreted as an equispaced quadrature of (3.24) at integer frequencies k ∈ I 2M .
Remark 3.10. We remark that for deriving the quadrature rule (3.27) from (3.24), we implicitly truncate the integral bounds in (3.24) as i. e., instead of (3.25) we rather deal with a distribution induced bỹ

28)
However, an analogous implication as in Theorem 3.9 cannot be shown when using the function (3.28) instead of (3.25).
As seen in Remark 3.10, our numerical method for computing suitable weights w j ∈ C, j = 1, . . . , N , can also be derived by means of a quadrature formula applied to the property δ,φ = Tξ,φ withξ defined in (3.28). Having a closer look at this property, the following equivalent characterization can be shown.
Theorem 3.11. Let a bandwidth M ∈ N d , nonequispaced points x j ∈ − 1 2 , 1 2 d and quadrature weights w j ∈ C, j = 1, . . . , N, be given. Then the following two statements are equivalent.
with the d-variate sinc function sinc(x) := d t=1 sinc(x t ) and Proof. By the definition T , ϕ = T,φ of the Fourier transform of a tempered distribution T ∈ S (R d ) we have 1, ϕ = δ,φ , cf. [53,Ex. 4.46]. Moreover, the distribution induced by (3.28) can be rewritten using the Fourier transform (3.17) as The inner integral can be determined by such that (3.29) shows the equality Tξ,φ = T ψ , ϕ . Hence, the assertions (i) and (ii) are equivalent.
Thus, since the statements (i) and (ii) of Theorem 3.11 are equivalent, one could also consider the property (ii) for deriving a numerical method to compute suitable weights w j ∈ C, j = 1, . . . , N . To this end, we have a closer look at for all ϕ ∈ S (R d ). Due to the integrals on both sides of (3.30), we need to discretize twice and therefore use the same quadrature rule on both sides of (3.30). For better comparability to (3.27) we utilize the same number |I 2M | of equispaced quadrature points y := (2M ) −1 , ∈ I 2M , as in (3.27), i. e., we consider In order that this applies for all ϕ ∈ S (R d ), we need to satisfy i. e., one could also compute weights w j ∈ C, j = 1, . . . , N , as a least squares solution to the linear system of equations (3.31). Hence, it merely remains the comparison of the two computation schemes. 1 · e 2πiky = |I 2M | · δ 0,k , k ∈ I 2M , we obtain the transformed system sinc (2M π (x j − y )) e 2πiky , k ∈ I 2M . (3.32) Comparing this linear system of equations to (3.14), we recognize an identical structure. Hence, we have a closer look at the connection between the expressions e 2πikx j and ∈I 2M sinc (2M π (x j − y )) e 2πiky . For this purpose, we consider the function f (t) = e 2πitx , t ∈ [−M, M ) d , for fixed x ∈ C d . By means off (t) := k∈Z d f (t + 2M k) we extend it into a (2M )-periodic function. This periodized version then possesses the Fourier coefficients cf. (2.2), i. e., the Fourier expansion off (t), t ∈ [−M, M ) d , for fixed x is given by cf. (2.1). Sincef (t) is continuous and piecewise differentiable, this Fourier series converges absolutely and uniformly, cf. [49,Ex. 1.22]. Thereby, we may consider the point evaluations at x = x j , j = 1, . . . , N , and t = k ∈ I 2M , such that we obtain the representation Thus, we recognize that (3.32) is a truncated version of (3.33). In other words, this implies that the linear system (3.14) is equivalent to a discretization of (3.30) incorporating infinitely many points y ∈ R d in (3.31).
Remark 3.13. For bandlimited functions several fast evaluation methods including the sinc function are known. The classical sampling theorem of Shannon-Whittaker-Kotelnikov, see [72,64,42], states that any bandlimited function f ∈ L 1 (R d ) ∩ C 0 (R d ) with maximum bandwidth M can be recovered from its uniform samples f (L −1 ), ∈ Z d , with L ≥ M , L := L · 1 d , and we have Since the practical use of this sampling theorem is limited due to the infinite number of samples, which is impossible in practice, and the very slow decay of the sinc function, various authors such as [57,68,50,47,40] considered the regularized Shannon sampling formula with localized sampling instead. Here ϕ m : R d → [0, 1] is a compactly supported window function with truncation parameter m ∈ N \ {1}, such that for ϕ m with small support the direct evaluation of (3.35) is efficient, see [40] for the relation to the NFFT window functions. On the other hand, in the one-dimensional setting a fast sinc transform was introduced in [39], which is based on the Clenshaw-Curtis quadrature sinc Lπ x − L ≈ n k=0 w k e −πiL(x− L )z k using Chebyshev points z k = cos( kπ n ) ∈ [−1, 1], k = 0, . . . , n, and corresponding Clenshaw-Curtis weights w k > 0. Thereby, sums of the form with uniform truncation parameter T ∈ 2N, can efficiently be approximated by means of fast Fourier transforms. More precisely, for the term in brackets one may utilize an NFFT, cf. (2.5). Then the resulting outer sum can be computed using an NNFFT, also referred to as NFFT of type III, see [17,20,46]

General error bound
In this section we summarize our previous findings by presenting a general error bound on density compensation factors computed by means of (3.14), that applies to trigonometric polynomials, 1-periodic functions f ∈ L 2 (T d ) ∩ C(T d ) and bandlimited functions f ∈ L 1 (R d ) ∩ C 0 (R d ) as well. w j e 2πikx j = δ 0,k + ε k , k ∈ I 2M , (3.36) with small ε k ∈ R for all k ∈ I 2M . Then there exists an ε ≥ 0 such that the corresponding density compensation procedure with W = diag(w j ) N j=1 satisfies the following error bounds.
(i) For any trigonometric polynomial f ∈ L 2 (T d ) of degree M given in (2.8) we have
(ii) For any 1-periodic function

38)
wheref := (c k (f )) k∈I M are the first |I M | coefficients given in (2.1) and p M is the best approximating trigonometric polynomial of degree M of f .
Proof. We start with some general considerations that are independent of the function f . By (3.36) we can find ε := max k∈I M |ε k | ≥ 0, such that |ε k | ≤ ε, k ∈ I 2M , and thereby N j=1 w j e 2πikx j − δ 0,k ≤ ε, k ∈ I 2M .
Then for all k, ∈ I M with ( − k) ∈ I 2M this yields

40)
and Considering the approximation error of (3.6), it can be estimated by Using A ∈ C N ×|I M | from (2.7) as well as W = diag(w j ) N j=1 = diag(w) we have Hence, from (3.40) -(3.43) and · 2 ≤ · F it follows that for p ∈ {1, 2, ∞} with 1 p + 1 q = 1. Now it merely remains to estimate Af − f p for the specific choice of f .
(i): Since a trigonometric polynomial (2.8) of degree M satisfies Af = f , the second error term in (3.44) vanishes and we obtain the assertion (3.37).
(ii): When considering a general 1-periodic function with the pointwise quadrature error of the uniform quadrature rule (3.19). For detailed investigations of quadrature errors for bandlimited functions we refer to [39,27]. Hence, we obtain Af − f p ≤ N 1/p Q C(T d ) and by (3.44) the assertion (3.39).
By Corollary 3.4 it is known that in the setting of trigonometric polynomials there is a linkage between an exact reconstruction (3.8) and the matrix product A * W A being equal to identity I |I M | . The following theorem shows that the error of the reconstruction (3.6) also affects the condition of the matrix A * W A.
and ε ≥ 0 be given as in Theorem 3.14. If additionally ε |I M | < 1 is fulfilled, then we have for the condition number κ 2 (X) := X 2 X −1 2 .
Proof. To estimate the condition number κ 2 (A * W A) we need to determine the norms A * W A 2 and (A * W A) −1 2 . By (3.36) it is known that A * W A = I |I M | + E, where E := (ε −k ) ,k∈I M , and therefore we have (3.47) Moreover, it is known by the theory of Neumann series, cf. [67,Thm. 4.20], that if I |I M | − T 2 < 1 holds for a matrix T ∈ C |I M |×|I M | , then T is invertible and its inverse is given by Using this property for T = A * W A we have Additionally, we know that |ε k | ≤ ε, k ∈ I 2M , with some ε > 0 and therefore In other words, the correctness of (3.48) is ensured if ε |I M | < 1. Since the spectral norm is a sub-multiplicative norm, (3.50) also implies E n 2 ≤ E n 2 ≤ (ε |I M |) n . Consequently, we have

Connection to certain density compensation approaches from literature
In literature a variety of density compensation approaches can be found that are concerned with the setting of bandlimited functions and make use of a sinc transform instead of the Fourier transform (3.6). Namely, instead of directly using the quadrature (3.21) for reconstruction, in these methods it is inserted into the inverse Fourier transform (3.18), i. e., By using the sinc matrix C ∈ R N ×|I M | from (3.52), the weight matrix W = diag(w j ) N j=1 as well as the vectors f = (f (x j )) N j=1 andf = (f (M −1 )) l∈I M , the evaluation of (3.53) at equispaced points M −1 , ∈ I M , can be denoted asf ≈ C * W f . Using the equispaced quadrature rule in (2.3), we find that evaluationsf (k) of (3.17) at the uniform grid points k ∈ I M can be approximated by (3.20) by means of a simple FFT. In matrix-vector notation this can be written asf ≈D * F * |I M |f wheref = (f (k)) k∈I M , F |I M | := (e 2πik(M −1 ) ) , k∈I M , cf. (2.14), andD := 1 |I M | I |I M | . Thus, all in all one obtains an approximation of the formf ≈D * F * |I M | C * W f . Here some of these approaches, cf. [19], shall be reconsidered in the context of the Fourier transform (3.6). We especially focus on the connection of the approaches among each other as well as to our new method introduced in Section 3.1.

Density compensation using the pseudoinverse
Since (3.4) is in general not exactly solvable, we study the corresponding least squares problem, instead, i. e., we look for the approximant that minimizes the residual norm f − Af 2 . It is known (e. g. [8, p. 15]) that this problem always has the unique solution with the Moore-Penrose pseudoinverse A † . Comparing (3.54) to the density compensation approach (3.6), the weights w j should be chosen such that the matrix product A * W approximates the pseudoinverse A † as best as possible, i. e., we study the optimization problem Minimize where · F denotes the Frobenius norm of a matrix. It was shown in [62] that the solution to this least squares problem can be computed as However, since a singular value decomposition is necessary for the calculations in (3.56), we obtain a high complexity of O(N 2 |I M | + |I M | 3 ). Therefore, we study some more sophisticated least squares approaches in the following.

Density compensation using weighted normal equations of first kind
It is known, that every least squares solution to (3.4) Analogous to [62] this could also be derived from (3.55) by introducing a right-hand scaling in the domain of measured data and minimizing the Frobenius norm of the weighted error matrix E r := E · A, where E := A * W − A † is the error matrix in (3.55).
In [59] it was shown that a solution W = diag(w) to (3.57) can be obtained by solving However, since Sw = b is not separable for single w j , j = 1, . . . , N , computing these weights is of complexity O(N 3 ). This is why the authors in [62] restricted themselves to a maximal image size of 64 × 64 pixels, which corresponds to setting M = 64.

Density compensation using weighted normal equations of second kind
Another approach for density compensation factors is based on the weighted normal equations of second kind We recognize that by (3.59) we are given an exact approximationf = A * W f of the Fourier coefficients in (3.6) in case y = f , and thereby AA * W = I N . To this end, we consider the optimization problem As in Section 3.4.2, we remark that this optimization problem (3.60) could also be derived from (3.55) by introducing an additional left-hand scaling in the Fourier domain and minimizing the Frobenius norm of the weighted error matrix E l := A · E.
Remark 3.16. An analogous approach considering the sinc transform (3.52) instead of the Fourier transform (3.6) was already studied in [52]. Another version using a sinc transform evaluated at pointwise differences of the nonequispaced points instead of (3.52) was studied in [12,30], where it was claimed that this approach coincides with the one in [52]. However, we remark that due to the sampling theorem of Shannon-Whittaker-Kotelnikov, see (3.34), applied to the function f (x) = sinc(M π(x j − x)), i. e., and its evaluation at x = x h , h = 1, . . . , N , this claim only holds asymptotically for |I M | → ∞ in the setting of the sinc transform. In contrast, when using the Fourier transform (3.6) this equality can directly be seen. Then the analog to [30] utilizes an approximation of the form f ≈ HW f , where the matrix H is defined as the system matrix of (3.4) evaluated at pointwise differences of the nonequispaced points, i. e., . (3.61) Since by (2.7) we have H = AA * , minimizing the approximation error leads to the optimization problem (3.60) as well.
It was shown in [52] that the minimizer of (3.60) is given by As mentioned in [52] one could also consider a simplified version of the optimization problem (3.60) by reducing the number of conditions, e. g. by summing the columns on both sides of AA * W = I N as (3.63) By means of (2.7) this can be written as AA * w = 1 N . Since fast multiplication with A and A * can be realized using the NFFT (see Algorithm 2.1) and the adjoint NFFT (see Algorithm 2.4), respectively, a solution to the linear system of equations AA * w = 1 N can be computed iteratively with arithmetic complexity O(|I M | log(|I M |) + N ). Finally, we investigate the connection of this approach to our method introduced in Section 3.1. To this end, suppose the linear system (3.10) is fulfilled for given w ∈ C N , i. e., by A * = A T we have (δ 0,k ) k∈I M = A T w = A * w. Then multiplication with A ∈ C N ×|I M | in (2.7) yields In other words, an exact solution w to the linear system (3.10) implies that the conjugate complex weights w exactly solve the system (3.63). However, the reversal does not hold true and therefore (3.63) is not equivalent to (3.10). Moreover, we have seen in Corollary 3.4 that an augmented variant of (3.10), namely (3.14), is necessary to obtain an exact reconstructionf k = h w k in (3.6) for trigonometric polynomials (2.8) with maximum degree M .

Direct inversion using matrix optimization
As seen in Remark 3.1, the previously considered density compensation techniques can be regarded as an optimization of the sparse matrix B ∈ R N ×|I M σ | from the NFFT, cf. Section 2.1. Since density compensation allows only N degrees of freedom, this limitation shall now be softened, i. e., instead of searching for optimal scaling factors for the rows of B, we now study the optimization of each nonzero entry of the sparse matrix B, cf. [38]. To this end, we firstly have another look at the equispaced setting. It is known by (3.1) and (3.2), that for equispaced points and appropriately chosen parameters a one-sided inversion is given by composition of the Fourier matrix and its adjoint. Hence, we aim to use this result to find a good approximation of the inverse in the general setting.
Considering problem (3.4) we seek to find an appropriate matrix X such that we have XA ≈ I |I M | , since then we can simply compute an approximation of the Fourier coefficients by means of Xf = XAf ≈f . To find this left-inverse X, we utilize the fact that in the equispaced case it is known that (3.1) holds in the overdetermined setting |I M | ≤ N . In addition, we also incorporate the approximate factorization A * ≈ D * F * B * of the adjoint NFFT, cf. Section 2.2, with the matrices D ∈ C |I M |×|I M | , F ∈ C |I M σ |×|I M | and B ∈ R N ×|I M σ | defined in (2.13), (2.14) and (2.15). Combining both ingredients we aim for an approximation of the form D * F * B * A ≈ I |I M | . To achieve an approximation like this, we aim to modify the matrix B such that its sparse structure with at most (2m + 1) d entries per row and consequently the arithmetic complexity of its evaluation is preserved. A matrix satisfying this property we call (2m + 1) d -sparse.
Remark 4.1. We remark that this approach can also be deduced from the density compensation method in Section 3 as follows. By Corollary 3.4 it is known that an exact reconstruction needs to satisfy A * W A = I |I M | . Since the reconstruction shall be realized efficiently by means of an adjoint NFFT, one rather studies D * F * B * W A ≈ I |I M | . Using the definitionB := W * B as in Remark 3.1, we end up with an approximation of the form D * F * B * A ≈ I |I M | . Thus, optimizing each nonzero entry of the sparse matrix B using this approximation is the natural generalization of the density compensation method from Section 3.
LetB denote such a modified matrix. By definingh := D * F * B * f , we recognize that the minimization of the approximation error implies the optimization problem Note that a similar idea for the forward problem, i. e., the evaluation of (2.5), was already studied in [51]. By the definition of the Frobenius norm we have Z F = Z * F , such that (4.2) is equivalent to its adjoint Since it is known by (2.14) that Thus, due to the fact that the Frobenius norm is a submultiplicative norm, we have Hence, we consider the optimization problem Based on the definition of the Frobenius norm of a matrix Z ∈ R k×n and the definition of the Euclidean norm of a vector y ∈ R n , we obtain for z j being the columns of Z ∈ R k×n that Since we aim to preserve the property that B is a (2m + 1) d -sparse matrix, we rewrite the norm in (4.5) by (4.6) in terms of the columns ofB considering only the nonzero entries of each column. To this end, analogously to (2.12) we define the index set of the nonzero entries of the -th column of B ∈ R N ×|I M σ | . Thus, we have (4.14) Thus, since 1 |I M σ | D −1 = diag(ŵ(k)) k∈I M , the computation of v involves neither multiplication with nor division by the (possibly) huge number |I M σ | and is therefore numerically stable.
This leads to the following algorithm.  Compute the right side v via (4.14) .
Note that a general statement about the dimensions of H ∈ C |I M |×|I M σ ,m ( )| is not possible, since the size of the set I M σ ,m ( ) heavily depends on the distribution of the points. To visualize this circumstance, we depicted some exemplary patterns of the nonzero entries of the original matrix B ∈ R N ×|I M σ | in Figure 4.1. It can easily be seen that for all choices of the points each row contains the same number of nonzero entries, i. e., all index sets (2.12) are of the same size of maximum (2m + 1) d . However, when considering the columns instead, we recognize an evident mismatch in the number of nonzero entries. We remark that because of the fact that each row of B ∈ R N ×|I M σ | contains at most (2m + 1) d entries, each column contains N |I M σ | (2m + 1) d entries on average. A general statement about the maximum size of the index sets (4.7) cannot be made. Roughly speaking, the more irregular the distribution of the points is, the larger the index sets (4.7) can be. Nevertheless, in general |I M σ ,m ( )| is a small constant compared to |I M |, such that Algorithm 4.3 ends up with total arithmetic costs of approximately O(|I M |).
In conclusion, our approach for an inverse NFFT can be summarized as follows.
Then there exists an ε ≥ 0 such that Moreover, the (asymmetric) Dirichlet kernel is the optimal window function for the inverse NFFT in Algorithm 4.4.
Proof. As in (4.1) the approximation error can be estimated by Using the same arguments as for (4.3) and (4.4) we proceed with To estimate the first Frobenius norm in (4.19), we rewrite it analogously to (4.8) columnwise as . Since b opt ∈ R |I M σ ,m | as solutions to the least squares problems (4.10) satisfy (4.15), we can find ε := max ∈I M σ ε ≥ 0, such that ε ≤ ε, ∈ I M σ , and thereby Therefore, we may write (4.19) as Thus, it remains to estimate the Frobenius norm F D 2 F . By the definitions of the Frobenius norm and the trace tr(Z) of a matrix Z, it is clear that Z 2 F = tr(Z * Z). Since by (2.14) we have that F * F = |I M σ | I |I M | , this yields Using the definition (2.13) of the diagonal matrix D ∈ R |I M |×|I M | , we obtain Then combination of (4.20), (4.21) and (4.22) implies such that (4.18) yields the assertion (4.16).
Since it is known that 0 ≤ŵ(k) ≤ 1, k ∈ I M , for suitable window functions of the NFFT, cf. [55], we have Hence, the smallest constant is achieved in (4.16) whenŵ(k) = 1, k ∈ I M , i. e., the (asymmetric) Dirichlet kernel (4.17) is the optimal window function for the inverse NFFT in Algorithm 4.4.
Note that for trigonometric polynomials (2.8) the error bound of Theorem 4.5 with the optimal window function (4.17) is the same as the error bound (3.37) from Theorem 3.14.
Remark 4.6. Up to now, we only focused on the problem (3.4). Finally, considering the inverse adjoint NFFT in (3.5), we remark that this problem can also be solved by means of the optimization procedure in Algorithm 4.3. Assuming again |I M | ≤ N , this is the underdetermined setting for the adjoint problem (3.5). Therefore, the minimum norm solution of (3.5) is given by the normal equations of second kind Incorporating the matrix decomposition of the NFFT, cf. Section 2.1, we recognize that a modification of the matrix B ∈ R N ×|I M σ | such that A * BF D ≈ I |I M | implies y ≈ h and hence f ≈ BF Dh. Thus, the optimization problem (4.3) is also the one to consider for (3.5). In other words, our approach provides both, an inverse NFFT as well as an inverse adjoint NFFT.

Numerics
Concluding, we have a look at some numerical examples. Besides comparing the density compensation approach from Section 3 to the optimization approach from Section 4, for both trigonometric polynomials (2.8) and bandlimited functions, we also demonstrate the accuracy of these approaches.
Remark 5.1. At first we introduce some exemplary grids. For visualization we restrict ourselves to the two-dimensional setting d = 2.
(i) We start with a sampling scheme that is somehow "close" to the Cartesian grid, but also possesses a random part. To this end, we consider the two-dimensional Cartesian grid and add a two-dimensional perturbation, i. e., . . , N 1 , j = 1, . . . , N 2 , and η 1 , η 2 ∼ U (−1, 1), where U (−1, 1) denotes the uniform distribution on the interval (−1, 1). A visualization of this jittered grid can be found in Figure 5.1a. Additionally, we also consider the random grid x t,j := 1 2 (η 1 , η 2 ) T , see Figure 5.1b. (ii) Moreover, we examine grids of polar kind, as mentioned in [23]. For R, T ∈ 2N the points of the polar grid are given by a signed radius r j := j R ∈ − 1 2 , 1 2 and an angle θ t := πt T ∈ − π 2 , π 2 as Since it is known that the inversion problem is ill-conditioned for this grid we consider a modification, the modified polar grid x t,j := r j (cos θ t , sin θ t ) T , (j, t) T ∈ I √ 2R × I T , (5.3) i. e., we added more concentric circles and excluded the points outside the unit square, see Figure 5.2a. Another sampling scheme which is known to lead to more stable results than the polar grid is the linogram or pseudo-polar grid, where the points lie on concentric squares instead of concentric circles, see Figure 5.2b. There we distinguish two sets of points, i. e., (5.4) (iii) Another modification of these polar grids was introduced in [33], where the angles θ t are not chosen equidistantly but are obtained by golden angle increments. For the golden angle polar grid we only exchange the equispaced angles of the polar grid to , t = 0, . . . , T − 1. The golden angle linogram grid is given by with θ t in (5.5), as illustrated in Figure 5.2c.
Before comparing the different approaches from Sections 3 and 4, we study the quality of our methods for the grids mentioned in Remark 5.1. More specifically, in Example 5.2 we investigate the accuracy of the density compensation method from Algorithm 3.7 with the weights introduced in Section 3.1, and in Example 5.3 we check if the norm minimization targeted in Section 4 is successful.
Example 5.2. Firstly, we examine the quality of our density compensation method in Algorithm 3.7 for a trigonometric polynomial f as in (2.8) with given Fourier coefficientŝ f k ∈ [1, 10], k ∈ I M . In this test we consider several d ∈ {1, 2, 3}. For the corresponding function evaluations of (2.8) at given points x j ∈ − 1 2 , 1 2 d , j = 1, . . . , N , we test how well these Fourier coefficients can be approximated. More precisely, we consider the estimateh w = D * F * B * W f , cf. (3.7), with the matrix W = diag (w j ) N j=1 of density compensation factors computed by means of Algorithm 3.6, i. e., by (3.15), in case |I 2M | ≤ N , or by (3.16), if |I 2M | > N , and compute the relative errors By (3.37) it is known that with the residual ε = A T |I 2M | w − e 0 ∞ ≥ 0, cf. (3.36). In our experiment we use random points x j with N t = 2 9−d , t = 1, . . . , d, cf. Figure 5.1b, and, for several problem sizes M = M · 1 d , M = 2 c with c = 1, . . . , 11 − d, we choose random Fourier coefficientsf k ∈ [1, 10], k ∈ I M . Afterwards, we compute the evaluations of the trigonometric polynomial (2.8) by means of an NFFT and use the resulting vector f as input for the reconstruction. Due to the randomness we repeat this 10 times and then consider the maximum error over all runs. The corresponding results are displayed in Figure 5.3. It can clearly be seen that |I 2M | < N , i. e., as long as M < N 1 2 = 2 8−d , the weights computed by means of (3.15) lead to an exact reconstruction of the given Fourier coefficients. However, as soon as we are in the setting |I 2M | > N the least squares approximation via (3.16) does not yield good results anymore. Having a look at the results for the grids in Remark 5.1, it becomes apparent that they separate into two groups. Figure 5.4a displays the results for the polar grid (5.2), which are the same as for the golden angle polar grid, cf. (5.5). In these cases there is only a slight improvement by the optimization. However, for all other mentioned grids the minimization procedure in Algorithm 4.3 is very effective. The results for these grids are depicted in Figure 5.4b exemplarily for the modified polar grid (5.3). Moreover, it can be seen that our optimization procedure in Algorithm 4.3 is most effective in the overdetermined setting |I M | ≤ N . One reason for the different behavior of polar and modified polar grid could be the ill-posedness of the inversion problem for the polar grid, which becomes evident in huge condition numbers of H * H , whereas the problem for modified polar grids is well-posed. Another reason can be found in the optimization procedure itself. Having a closer look at the polar grid, see Figure 5.2a, there are no grid points in the corners of the unit square. Therefore, some of the index sets I M σ ,m ( ), cf. (4.7), are empty and no optimization can be done for the corresponding matrix columns. This could also cause the worsened minimization properties of the polar grid.
Next, we proceed with comparing the density compensation approach from Section 3 using the weights w j introduced in Section 3.1 to the optimization approach for modifying the matrix B from Section 4. To this end, we show an example concerning trigonometric polynomials (2.8) of degree M and a second one that deals with bandlimited functions of bandwidth M . Here we restrict ourselves to the two-dimensional setting d = 2 for better visualization of the results.
Example 5.4. Similar to [3,38] we have a look at the reconstruction of the Shepp-Logan phantom, see Figure 5.5a. Here we treat the phantom data as given Fourier coefficientŝ f := (f k ) k∈I M of a trigonometric polynomial (2.8). For given points x j ∈ − 1 2 , 1 2 2 , j = 1, . . . , N , we then compute the evaluations of the trigonometric polynomial (2.8) by means of an NFFT and use the resulting vector as input for the reconstruction.
In a first experiment, we test the inversion methods from Sections 3 and 4 as in [3] for increasing input sizes. To this end, we choose M = (M, M ) T , M = 2 c with c = 3, . . . , 10, and linogram grids (5.4) of size R = 2M , T = 2R, i. e., we consider the setting |I 2M | < N . For using Algorithm 4.4 we choose the oversampling factor σ = 1.0 and the truncation parameter m = 4. For each input size we measure the computation time of the precomputational steps, i. e., the computation of the weight matrix W or the computation of the optimized sparse matrix B opt , as well as the time needed for the reconstruction, i. e., the corresponding adjoint NFFT, see Algorithms 3.7 and 4.4. Moreover, for the reconstructionh ∈ {h w , h opt }, cf. (3.7) and (4.12), we consider the relative errors The corresponding results can be found in Table 5.1. We remark that since we are in the setting |I 2M | < N , the density compensation method in Algorithm 3.7 with weights computed by (3.15) indeed produces nearly exact results. Although, our optimization procedure from Algorithm 4.4 achieves small errors as well, this reconstruction is not as good as the one by means of our density compensation method. Note that in comparison to [3] our method in Algorithm 3.7 using density compensation produces errors of the same order, but is much more effective for solving several problems using the same points x j for different input values f . Since our precomputations have to be done only once in this setting, we strongly profit from the fact that we only need to perform an adjoint NFFT as reconstruction, which is very fast, whereas in [3] they would need to execute their whole routine each time again.
As a second experiment we aim to decrease the amount of overdetermination, i. e., we want to keep the size |I M | of the phantom, but reduce the number N of the points x j , j = 1, . . . , N . To this end, we now consider linogram grids (5.4) of the smaller size R = M , T = 2R, i. e., now we have |I 2M | > N . The reconstruction of the phantom of size 1024 × 1024 is presented in Figure 5.5 (top) including a detailed view of the 832nd row of this reconstruction (bottom). Despite the reconstruction via Algorithm 4.4 as well as the density compensation method using weights computed by means of (3.16), we also considered the result using Voronoi weights. For all approaches we added the corresponding relative errors (5.9) to Figure 5.5 as well.
Due to the fact that the exactness condition |I 2M | < N (cf. Section 3.1.1) is violated, it can be seen in Figure 5.5c that the density compensation method using weights computed by means of (3.16) does not yield an exact reconstruction is this setting. On the contrary, we recognize that our optimization method, see Figure 5.5d, achieves a huge improvement in comparison to the density compensation techniques in Figure 5.5b and 5.5c since no artifacts are visible. Presumably, this arises because there are more degrees of freedom  in the optimization of the matrix B from Section 4 than with the density compensation techniques from Section 3, cf. Remark 3.1. We remark that although the errors are not as small as in Table 5.1, by comparing Figure 5.5a and 5.5d it becomes apparent that the differences are not even visible anymore. Note that for this result the number N of points is ca. 4 times lower as for the results in depicted in Table 5.1, i. e., we only needed twice as much function values as Fourier coefficients, whereas e. g. in [3] they worked with a factor of more than 4.
Example 5.5. Finally, we examine the reconstruction properties for bandlimited functions f ∈ L 1 (R d ) ∩ C 0 (R d ) with maximum bandwidth M . To this end, we firstly specify a compactly supported functionf and consequently compute its inverse Fourier transform f , such that its samples f (x j ) for given x j ∈ − 1 2 , 1 2 2 , j = 1, . . . , N , can be used for the reconstruction of the samplesf (k), k ∈ I M . Here we consider the tensorized functionf (v) = g(v 1 ) · g(v 2 ), where g(v) is the one-dimensional triangular pulse g(v) : . Then for all b ∈ N with b ≤ M 2 the associated inverse Fourier transform (v) e 2πivx dv = b 2 sinc 2 (bπx) = b 2 sinc 2 (bπx 1 ) sinc 2 (bπx 2 ), x ∈ R 2 , is bandlimited with bandwidth M . In this case, we consider M = 64 and b = 24 as well as the jittered grid (5.1) of size N 1 = N 2 = 144, i. e., we study the setting |I 2M | ≤ N . Now the aim is comparing the different density compensation methods considered in Section 3 and the optimization approach from Section 4. More precisely, we consider the reconstruction using Voronoi weights, the weights computed via (3.58), the weights in easily be seen that Voronoi weights, see Figure 5.6a, and the weights in (3.62), see Figure 5.6d, do not yield a good reconstruction, as expected. The other three approaches produce nearly the same reconstruction error, which is also obtained by reconstruction on an equispaced grid and therefore is the best possible. In other words, in case of bandlimited functions the truncation error in (3.20) is dominating and thus reconstruction errors smaller than the ones shown in Figure 5.6 cannot be expected. Note that the comparatively small choice of M = 64 was made in order that the computation of the weights via (3.58), see Figure 5.6c, as well as the weights in (3.62), see Figure 5.6d, is affordable, cf. Section 3.4.2. In contrast, our new methods using Algorithm 3.7, see Figure 5.6b, or Algorithm 4.4, see Figure 5.6e, are much more effective and therefore better suited for the given problem.

Conclusion
In the present paper we considered several direct methods for computing an inverse NFFT, i. e., reconstructing the Fourier coefficientsf k , k ∈ I M , from given nonequispaced data f (x j ), j = 1, . . . , N . Being a direct method here means, that for a fixed set of points x j , j = 1, . . . , N , the reconstruction can be realized with the same number of arithmetic operations as a single application of an adjoint NFFT (see Algorithm 2.4). As we have seen in (3.3), a certain precomputational step is compulsory, since the adjoint NFFT does not yield an inversion by itself. Although this precomputations might be rather costly, they need to be done only once for a given set of points x j , j = 1, . . . , N . Therefore, direct methods are especially beneficial in case of fixed points. For this reason, we studied two different approaches of this kind and especially focused on methods for the multidimensional setting d ≥ 1 that are applicable for general sampling patterns.
Firstly, in Section 3 we examined the well known approach of sampling density compensation, where suitable weights w j ∈ C, j = 1, . . . , N, are precomputed, such that the reconstruction can be realized by means of an adjoint NFFT applied to the scaled data w j f (x j ). We started our investigations with trigonometric polynomials in Section 3.1. In Corollary 3.4 we introduced the main formula (3.14), that yields exact reconstruction for all trigonometric polynomials of degree M . In addition to this theoretical considerations, we also discussed practical computation schemes for the overdetermined as well as the underdetermined setting, as summarized in Algorithm 3.6. Afterwards, in Section 3.2 we studied the case of bandlimited functions, which often occurs in the context of MRI, and discussed that the same numerical procedures as in Section 3.1 can be used in this setting as well. In Section 3.3 we then summarized the previous findings by presenting a general error bound on density compensation factors computed by means of Algorithm 3.6 in Theorem 3.14. In addition, this also yields an estimate on the condition number of the matrix product A * W A, as shown in Theorem 3.15. In Section 3.4 we surveyed certain approaches from literature and commented on their connection among each other as well as to the method presented in Section 3.1.
Subsequently, in Section 4 we studied another direct inversion method, where the matrix representation A ≈ BF D of the NFFT is used to modify the sparse matrix B, such that a reconstruction is given byf ≈ D * F * B * opt f . In other words, the inversion is done by a modified adjoint NFFT, while the optimization of the matrix B can be realized in a precomputational step, see Algorithm 4.4.
Finally, in Section 5 we had a look at some numerical examples to investigate the accuracy of the previously introduced methods. We have seen that our approaches are best-suited for the overdetermined setting |I M | ≤ N and work for many different sampling patterns. More specifically, in the highly overdetermined case |I 2M | ≤ N we have theoretically proven as well as numerically verified in several examples that the density compensation technique in Algorithm 3.7 leads to an exact reconstruction for trigonometric polynomials. In case not that much data is available and we have to reduce the amount of overdetermination such that |I 2M | > N , we have shown that the optimization approach from Algorithm 4.4 is preferable, since the higher number of degrees of freedom in the optimization (see Remark 3.1) yields better results. In addition, also for the setting of bandlimited functions we demonstrated that our methods are much more efficient than existing ones.