## ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 28 June 2023
Sec. Mathematics of Computation and Data Science

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

• Faculty of Mathematics, Chemnitz University of Technology, Chemnitz, Germany

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.

## 1. 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

with given Fourier coefficients , $k\in {\mathcal{I}}_{M}$, at nonequispaced points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, j = 1,…, N, N∈ℕ, where ${\mathcal{I}}_{M}:={ℤ}^{d}\cap {\left[-\frac{M}{2},\frac{M}{2}\right)}^{d}$ with $|{\mathcal{I}}_{M}|={M}^{d}$. In case we are given equispaced points xj and $|{\mathcal{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. [1, 2] and solution of PDEs, cf. [3] need to perform an inverse nonequispaced fast Fourier transform (iNFFT), i.e., compute the Fourier coefficients from given function evaluations f(xj) of the trigonometric polynomial (1.1). Hence, we are interested in an inversion also for nonequispaced data.

In general, the number N of points xj is independent of the number $|{\mathcal{I}}_{M}|$ of Fourier coefficients and hence the nonequispaced Fourier matrix

which we would have to invert, is rectangular in most cases. Considering the corresponding linear system $A\stackrel{^}{f}=f$ with $f:={\left(f\left({x}_{j}\right)\right)}_{j=1}^{N}$ and $\stackrel{^}{f}:={\left({\stackrel{^}{f}}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$, this can either be overdetermined, if $|{\mathcal{I}}_{M}|\le N$, or underdetermined, if $|{\mathcal{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 [48] 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.

### 1.1. Iterative inversion methods

We start surveying the iterative methods. For the one-dimensional setting d = 1 with $|{\mathcal{I}}_{M}|=N$ an algorithm was published in Ruiz-Antolin and Townsend [9], 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 $|{\mathcal{I}}_{M}|\le N$ can be found in Feichtinger et al. [4], where the Toeplitz structure of the matrix product A*WA with a diagonal matrix $W:=diag{\left({w}_{j}\right)}_{j=1}^{N}$ of Voronoi weights is used to compute the solution iteratively by means of the CG algorithm.

For higher dimensional problems with d ≥ 1, several approaches compute a least squares approximation to the linear system $A\stackrel{^}{f}=f$. In the overdetermined case $|{\mathcal{I}}_{M}|\le N$, the given data can typically only be approximated up to a residual $r:=A\stackrel{^}{f}-f$. Therefore, the weighted least squares problem

$Minimizef^∈ℂ|IM|∑j=1Nwj|∑k∈IMf^ke2πikxj-f(xj)|2$

is considered, which is equivalent in solving the weighted normal equations of the first kind ${A}^{*}WA\stackrel{^}{f}={A}^{*}Wf$ with the diagonal matrix $W:=diag{\left({w}_{j}\right)}_{j=1}^{N}$ of weights in the time domain. In [1012], these normal equations are solved iteratively by means of CG using the NFFT to realize fast matrix-vector multiplications involving A, whereas in Pruessmann and Wayer [13] a fast convolution is used.

In the consistent underdetermined case $|{\mathcal{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

$Minimizef^∈ℂ|IM|∑k∈IM|f^k|2w^k subject to Af^=f.$

This interpolation problem is equivalent to the weighted normal equations of the second kind AŴA*y = f, $\stackrel{^}{f}=Ŵ{A}^{*}y$ with the diagonal matrix $Ŵ:=diag{\left({ŵ}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$ of weights in the frequency domain. In Kunis and Potts [14], the CG method was used in connection with the NFFT to iteratively compute a solution to this problem, see also [15, Section 7.6.2].

### 1.2. Regularization methods

Moreover, there also exist several regularization techniques for the multidimensional setting d ≥ 1. For example, [1618] all solve the 1-regularized problem

$Minimizef^∈ℂ|IM|12∥Af^-f∥22+λ∥Lmf^∥1$

with regularization parameter λ > 0 and the m-th order polynomial annihilation operator ${\mathcal{L}}^{m}\in {ℝ}^{N×|{\mathcal{I}}_{M}|}$ as sparsifying transform, see Archibald et al. [19]. Based on this, weighted p-schemes

$Minimizef^∈ℂ|IM|12∥Af^-f∥22+1p∥WLmf^∥pp$

were introduced in [2023], which are designed to reduce the penalty at locations where ${\mathcal{L}}^{m}\stackrel{^}{f}$ is nonzero. For instance, Churchill et al. [24] and Scarnati and Gelb [25] each state a two-step method, that first uses edge detection to create a mask, i.e., a weighting that which indicates where nonzero 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.

### 1.3. Direct inversion methods

In contrast to these iterated procedures, there are also so-called direct methods that do not require multiple iteration steps. Already in Dutt and Rokhlin [26] a direct method was explained for the setting d = 1 and $|{\mathcal{I}}_{M}|=N$ that 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 [27], or utilize fast summation in Kircheis and Potts [28] for the fast evaluation of occurring sums, see also [15, Section 7.6.1].

In the overdetermined setting $|{\mathcal{I}}_{M}|\le 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 Heinig and Rost [29], can be used to solve the normal equations ${A}^{*}A\stackrel{^}{f}={A}^{*}f$ exactly by analogy with Averbuch et al. [72]. 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 Gelb and Song [30], 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 Gelb and Song [30] is based only on optimizing a diagonal matrix [the matrix D defined in (2.13)], whereas in Kircheis and Potts [28] 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 3B, such that the inversion involves only one-dimensional FFTs and interpolations. On the one hand, in Averbuch et al. [31], a least squares solution is computed iteratively by using a fast multiplication technique with the matrix A, which can be derived in the case of the linogram grid. On the other hand, [32] 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 to the special case of the linogram grid, see Figure 3B, or the polar grid by another resampling, cf. Figure 3A. Since we are interested in more generally applicable methods, a brief introduction to direct inversion for general sampling patterns can be found in Kircheis and Potts [33].

### 1.4. Current approach

In this article, we focus on direct inversion methods that are applicable to general sampling patterns and present new methods for the computation of an iNFFT. Note that the direct method in the context of the linear system $A\stackrel{^}{f}=f$ means, that for a fixed set of points xj, j = 1,…, N, the reconstruction of $\stackrel{^}{f}$ 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 these precomputations might be rather costly, they need to be done only once for a given set of points xj, j = 1,…, N. Therefore, direct methods are especially beneficial in the case of fixed points for several measurement vectors f.

For this reason, the current paper is concerned with two different approaches of this kind. First, we consider the very well-known approach of so-called sampling density compensation, which can be written as $\stackrel{^}{f}\approx {A}^{*}Wf$ with a diagonal matrix $W:=diag{\left({w}_{j}\right)}_{j=1}^{N}$ of weights. The already mentioned precomputations then consist of computing suitable weights wj, while the actual reconstruction step includes only one adjoint NFFT applied to the scaled coefficient vector Wf. In this article, we examine several existing approaches for computing the weights wj 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 Kircheis and Potts [33]. Therefore, the idea is using the matrix representation ABFD of the NFFT to modify one of the matrix factors, such that an inversion is given by $\stackrel{^}{f}\approx {D}^{*}{F}^{*}{B}_{\text{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.

### 1.5. Outline of this paper

This article is organized as follows. In Section 2, we introduce the already mentioned algorithm, the NFFT, as well as its adjoint version, the adjoint NFFT. Second, 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 and 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. Afterward, we examine another direct inversion method in Section 4, where we aimed 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.

## 2. Nonequispaced fast Fourier transform

Let

$𝕋d: =ℝd\ℤd≅[-12,12)d={x∈ℝd:-12≤xt<12,t=1,…,d}$

denote the d-dimensional torus with d ∈ ℕ. For M := (M,…, M)T, M ∈ 2ℕ, we define the multi-index set

$IM: =ℤd∩[-M2,M2)d={k∈ℤd:-M2≤kt

with cardinality $|{\mathcal{I}}_{M}|={M}^{d}$. The inner product of two vectors shall be defined as usual as kx := k1x1 + ⋯ + kdxd. and the componentwise product as $x\odot y:={\left({x}_{1}{y}_{1},\dots ,{x}_{d}{y}_{d}\right)}^{T}$. In addition, we define the all ones vector ${\text{1}}_{d}:={\left(1,\dots ,1\right)}^{T}$ and the reciprocal of a vector x with nonzero components shall be given by ${x}^{-1}:={\left({x}_{1}^{-1},\dots ,{x}_{d}^{-1}\right)}^{T}$.

We consider the Hilbert space ${L}_{2}\left({𝕋}^{d}\right)$ of all 1-periodic, complex-valued functions, which possess the orthonormal basis {e2πikx : k ∈ ℤd}. It is known that every function $f\in {L}_{2}\left({T}^{d}\right)$ is uniquely representable in the form

with the coefficients

where the sum in (2.1) converges to f in the ${L}_{2}\left({𝕋}^{d}\right)$-norm, cf. [15, Thm. 4.5]. A series of the form (2.1) is called Fourier series with the Fourier coefficients (2.2). Numerically, the Fourier coefficients (2.2) are approximated on the uniform grid $\left\{{M}^{-1}\odot \ell ,\ell \in {\mathcal{I}}_{M}\right\}$ by the trapezoidal rule for numerical integration as

which is an acceptable approximation for $k\in {\mathcal{I}}_{M}$, see e. g., [15, 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

with ${\stackrel{^}{f}}_{k}\approx {c}_{k}\left(f\right)$, $k\in {\mathcal{I}}_{M}$, can be realized by means of an inverse fast Fourier transform (iFFT), which is the same algorithm except for some reordering and scaling, cf. [15, Lem. 3.17].

Now, suppose we are given nonequispaced points ${x}_{j}\in {𝕋}^{d}$, j = 1,…, N, instead. Then, we consider the computation of the sums

for given ${\stackrel{^}{f}}_{k}\in ℂ$, $k\in {\mathcal{I}}_{M}$, as well as the adjoint problem of computing

for given fj ∈ ℂ, j = 1,…, N. By defining the nonequispaced Fourier matrix

as well as the vectors $f:={\left({f}_{j}\right)}_{j=1}^{N}$, $\stackrel{^}{f}:={\left({\stackrel{^}{f}}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$ and $h:={\left({h}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$, the computation of sums of the form (2.5) and (2.6) can be written as $f=A\stackrel{^}{f}$ and h = A*f, where ${A}^{*}:={\overline{A}}^{T}$ denotes the adjoint matrix of A.

Since the naive computation of (2.5) and (2.6) is of complexity $\mathcal{O}\left(N·|{\mathcal{I}}_{M}|\right)$, a fast approximate algorithm, the so-called nonequispaced fast Fourier transform (NFFT), is briefly explained below. For more information see [3438] or [15, p. 377–381].

### 2.1. NFFT

We first restrict our attention to the problem (2.5), which is equivalent to the evaluation of a trigonometric polynomial

with given ${\stackrel{^}{f}}_{k}\in ℂ$, $k\in {\mathcal{I}}_{M}$, at given nonequispaced points ${x}_{j}\in {𝕋}^{d}$, j = 1,…, N. Let $w\in {L}_{2}\left({ℝ}^{d}\right)\cap {L}_{1}\left({ℝ}^{d}\right)$ be a so-called window function, which is well localized in space and frequency domain. Now, we define the 1-periodic function $\stackrel{~}{w}\left(x\right):=\sum _{r\in {ℤ}^{d}}w\left(x+r\right)$ with absolute convergent Fourier series. As a consequence, the Fourier coefficients of the periodization $\stackrel{~}{w}$ have the form

$ck(w~)=∫𝕋dw~(x)e-2πikxdx=∫ℝdw(x)e-2πikxdx=:w^(k).$

For a given oversampling factor σ ≥ 1, we define 2ℕ ∋ Mσ:= 2⌈⌈σM⌉/2⌉ as well as Mσ:= Mσ·1d, and approximate f by a linear combination translates of the periodized window function, i.e.,

where g ∈ ℂ, $\ell \in {\mathcal{I}}_{{M}_{\sigma }}$, are coefficients to be determined such that (2.9) yields a good approximation. By means of the convolution theorem (see [15, Lem. 4.1]), the approximant ${s}_{1}\in {L}_{2}\left({𝕋}^{d}\right)$ in (2.9) can be represented as

where the discrete Fourier transform of the coefficients g is defined by

Comparing (2.5) and (2.10) then yields

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 ${\left[-\text{}m\text{}\text{}/\text{}\text{}{M}_{\sigma }\text{},\text{}m\text{}\text{}/\text{}\text{}{M}_{\sigma }\text{}\right]}^{d}$ with truncation parameter mMσ. In this case, w can be approximated by the compactly supported function

$wm(x): ={w(x):x∈[-mMσ,mMσ]d,0:otherwise.$

Thereby, we approximate s1 by the short sums

where the index set

contains at most (2m+1)d entries for each fixed xj. Thus, the obtained algorithm can be summarized as follows.

ALGORITHM 2.1

Algorithm 2.1. NFFT.

Remark 2.2. Suitable window functions can be found, e.g., in [34, 35, 3741].

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

• the diagonal matrix

• the truncated Fourier matrix

• 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 ABFD of (2.7), cf. [15, p. 383]. In other words, the NFFT uses the approximation

$e2πikxj≈1|IMσ|·w^(k)∑ℓ∈IMσ,m(xj)e2πik(Mσ-1⊙ℓ)w~m(xj-Mσ-1⊙ℓ).$

Remark 2.3. It has to be pointed out that because of consistency the factor $|{\mathcal{I}}_{{M}_{\sigma }}{|}^{-1}$ is here not located in the matrix F as usual but in the matrix D.

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

The algorithms presented in this section (Algorithms 2.1 and 2.4) are part of the NFFT software [42]. For algorithmic details we refer to Keiner et al. [38].

## 3. Direct inversion using density compensation

Having introduced the fast methods for nonequispaced data, we remark that various applications such as MRI and solution of PDEs are interested in the inverse problem, i.e., instead of the evaluation of (2.5) the aim was to compute the Fourier coefficients ${\stackrel{^}{f}}_{k}$, $k\in {\mathcal{I}}_{M}$, from given nonequispaced data f(xj), 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}=\frac{1}{n}j\in {T}^{d}$, $j\in {\mathcal{I}}_{n}$, with n := n·1d and $|{\mathcal{I}}_{n}|=N$, the nonequispaced Fourier matrix $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7) turns into the equispaced Fourier matrix $F\in {ℂ}^{|{\mathcal{I}}_{{M}_{\sigma }}|×|{\mathcal{I}}_{M}|}$ from (2.14) with $|{\mathcal{I}}_{{M}_{\sigma }}|=N$. Thereby, it results from the geometric sum formula that

as well as

and $|{\mathcal{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}\in {𝕋}^{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 coefficients $\stackrel{^}{f}={\left({\stackrel{^}{f}}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$ from given function values $f={\left(f\left({x}_{j}\right)\right)}_{j=1}^{N}$. This problem is referred to as inverse NDFT (iNDFT) and an efficient solver shall be called inverse NFFT (iNFFT).

(2) Solve the linear system

i.e., reconstruct the coefficients $f={\left({f}_{j}\right)}_{j=1}^{N}$ from given data $h={\left({h}_{k}\right)}_{k\in {\mathcal{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 $|{\mathcal{I}}_{M}|$ and N are independent, such that the nonequispaced Fourier matrix $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7) is generally rectangular.

At first, we restrict our attention to the 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 procedure, with so-called direct methods. In the setting of the problem (3.4) we hereby mean methods, where for a fixed set of points xj, j = 1,…, N, the reconstruction of $\stackrel{^}{f}$ 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 these precomputations might be rather costly, they need to be done only once for a given set of points xj, 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}\in {𝕋}^{d}$, j = 1,…, N. Thereby, the Fourier coefficients (2.2) are approximated by a general quadrature rule using quadrature weights wj ∈ ℂ, j = 1,…, N, which are needed for sampling density compensation due to the nonequispaced sampling. Thus, for a trigonometric polynomial (2.8), we have

Using the nonequispaced Fourier matrix $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7), the diagonal matrix of weights $W:=\text{diag}{\left({w}_{j}\right)}_{j=1}^{N}\in {ℂ}^{N×N}$ as well as the vector ${h}^{\text{w}}:={\left({h}_{k}^{\text{w}}\right)}_{k\in {\mathcal{I}}_{M}}$, the nonequispaced quadrature rule (3.6) can be written as $\stackrel{^}{f}\approx {h}^{\text{w}}:={A}^{*}Wf$. 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 by

with the matrices $D\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{M}|}$, $F\in {ℂ}^{|{\mathcal{I}}_{{M}_{\sigma }}|×|{\mathcal{I}}_{M}|}$, and $B\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ defined in (2.13), (2.14), and (2.15). In other words, for density compensation methods the already mentioned precomputations consists of computing the quadrature weights wj ∈ ℂ, j = 1,…, N, while the actual reconstruction step includes only one adjoint NFFT (see Algorithm 2.4) applied to the scaled measurement vector Wf.

The aim of all density compensation techniques was then to choose appropriate weights wj ∈ ℂ, 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 wj.

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., [43]. 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. First, in Section 3.1, we introduce density compensation factors wj, 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 (3.7) for bandlimited functions $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ as well. 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 to each other as well as to the method introduced in Section 3.1.

Remark 3.1. Recapitulating, we have a closer look at some possible interpretation perspectives on the reconstruction (3.7).

(i) If we define g := Wf, i.e., each entry of f is scaled with respect to the points xj, j = 1,…, N, the approximation (3.7) can be written as $\stackrel{^}{f}\approx {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 matrix $\stackrel{~}{B}:={W}^{*}B$, i.e., scaling the rows of B with respect to the points xj, j = 1, …, N, the approximation (3.7) can be written as $\stackrel{^}{f}\approx {D}^{*}{F}^{*}{\stackrel{~}{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\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ in (2.15), as it shall be done in Section 4. We remark that density compensation methods allow only N degrees of freedom.

### 3.1. Exact quadrature weights for trigonometric polynomials

Similar to Gräf et al. [44], we aimed to introduce density compensation factors wj, 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 first examine certain properties that arise from (3.6) being exact.

Theorem 3.2. Let a polynomial degree M ∈ (2ℕ)d, nonequispaced points ${x}_{j}\in {T}^{d}$, j = 1, …, N, and quadrature weights wj ∈ ℂ be given. Then an exact reconstruction formula (3.6) for trigonometric polynomials (2.8) with maximum degree M satisfying

implies the following equivalent statements.

is exact for all trigonometric polynomials (2.8) with maximum degree M.

(ii) The linear system of equations

is fulfilled with the matrix $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7) and $w:={\left({w}_{j}\right)}_{j=1}^{N}$.

Proof: (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}_{k}^{\text{w}}$ we proceed with

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

This together with the property (i) and (3.11) leads to

$f^0=∑k∈IMf^k∑j=1Nwje2πikxj$

and thus to assertion (3.10).

(ii) ⇒ (i): Combining (3.11), (3.10), and (3.13) yields the assertion via

$∫𝕋df(x)dx=∑k∈IMf^k·δ0,k=∑k∈IMf^k∑j=1Nwje2πikxj=∑j=1Nwjf(xj).$

Remark 3.3. Comparable results can also be found in literature. A fundamental theorem in numerical integration, see [45], states that for any integral ${\int }_{{𝕋}^{d}}f\left(x\right)dx$ there exists an exact quadrature rule (3.9), i.e., optimal points ${x}_{j}\in {𝕋}^{d}$ and weights wj ∈ ℂ, j = 1,…, N, such that (3.9) is fulfilled. In Gröchenig [46, Lemma 2.6], it was shown that for given points ${x}_{j}\in {𝕋}^{d}$, j = 1,…, N, certain quadrature weights wj can be stated by means of frame theoretical considerations that lead to an exact quadrature rule (3.9) by definition. Moreover, it was shown (cf. [46, 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 [47]. According to Huybrechs [47, Sec. 2.1], these quadrature weights wj, j = 1,…, N, can be found by solving a linear system of equations Φw = v, where Φk, j = ϕk(xj) and ${v}_{k}={\int }_{{𝕋}^{d}}{\varphi }_{k}\left(x\right)dx$ for a given set of basis functions ${\left\{{\varphi }_{k}\right\}}_{k\in {\mathcal{I}}_{M}}$. In our setting, we have ${\varphi }_{k}\left(x\right)={\text{e}}^{2\pi \text{i}kx}$ and, therefore,

$Φ=(e2πikxj)k,j=AT and vk=∫𝕋d1·e2πikxdx=δ0,k,$

i.e., the same linear system of equations as in (3.10). We remark that both Gröchenig [46] and Huybrechs [47] 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 guarantees (3.6) being exact for all trigonometric polynomials (2.8) with maximum degree M.

Corollary 3.4. The two statements (i) and (ii) in Theorem 3.2 are not equivalent to property (3.8), since (3.10) does not imply an exact reconstruction in (3.6).

However, an augmented variant of (3.10), namely,

yields an exact reconstruction in (3.6) for trigonometric polynomials (2.8) with maximum degree M. In addition, (3.14) implies the matrix equation ${A}^{*}WA={I}_{|{\mathcal{I}}_{M}|}$ with $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7) and the identity matrix ${I}_{|{\mathcal{I}}_{M}|}$ of size $|{\mathcal{I}}_{M}|$.

Proof: Utilizing definitions (3.6) and (2.8), we have

Since (3.10) only holds for $k,\ell \in {\mathcal{I}}_{M}$ with $\left(\ell -k\right)\in {\mathcal{I}}_{M}$, this implies

$hkw=f^k+∑ℓ∈IM(ℓ-k)∉IMf^ℓ∑j=1Nwje2πi(ℓ-k)xj, k∈IM,$

where for all $k\in {\mathcal{I}}_{M}\\left\{0\right\}$ there exists an $\ell \in {\mathcal{I}}_{M}$ with $\left(\ell -k\right)\in {\mathcal{I}}_{2M}\{\mathcal{I}}_{M}$.

As $\left(\ell -k\right)\in {\mathcal{I}}_{2M}$ for $k,\ell \in {\mathcal{I}}_{M}$, the augmented variant (3.14) yields

$hkw=∑ℓ∈IMf^ℓ∑j=1Nwje2πi(ℓ-k)xj=∑ℓ∈IMf^k·δ0,k=f^k, k∈IM.$

Moreover, since δ(k), 0 = δk,, the condition (3.14) implies

$δk,ℓ=∑j=1Nwje2πi(ℓ-k)xj=∑j=1Ne-2πikxj(wje2πiℓxj), k,ℓ∈IM.$

In matrix-vector notation, this can be written as ${A}^{*}WA={I}_{|{\mathcal{I}}_{M}|}$ with $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7) and the identity matrix ${I}_{|{\mathcal{I}}_{M}|}$ of size $|{\mathcal{I}}_{M}|$. We remark that this matrix equation immediately shows that we have an exact reconstruction of the form (3.8), since if ${A}^{*}WA={I}_{|{\mathcal{I}}_{M}|}$ is fulfilled, (3.4) implies that $\stackrel{^}{f}={A}^{*}WA\stackrel{^}{f}={A}^{*}Wf$.

Remark 3.5. Let $f\in {L}_{2}\left({𝕋}^{d}\right)$ be an arbitrary 1-periodic function (2.1). Then (3.14) yields

i.e., for a function $f\in {L}_{2}\left({𝕋}^{d}\right)$ we only have a good approximation in case the coefficients c(f) are small for $\ell \notin {\mathcal{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

So far, we have seen in Corollary 3.4 that an exact solution $w={\left({w}_{j}\right)}_{j=1}^{N}$ 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 aimed to use this condition (3.14) to numerically find optimal density compensation factors wj ∈ ℂ, j = 1,…, N.

Having a closer look at the condition (3.14) we recognize that it can be written as the linear system of equations with the matrix ${A}_{|{\mathcal{I}}_{2M}|}\in {ℂ}^{N×|{\mathcal{I}}_{2M}|}$, cf. (2.7), and right side ${e}_{0}:={\left({\delta }_{0,k}\right)}_{k\in {\mathcal{I}}_{2M}}$. We remark that in contrast to $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ we now deal with the enlarged matrix ${A}_{|{\mathcal{I}}_{2M}|}\in {ℂ}^{N×|{\mathcal{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 needs to be solved only once for fixed points ${x}_{j}\in {𝕋}^{d}$, j = 1,…, N. Its solution w can then be used to efficiently approximate $\stackrel{^}{f}$ for multiple measurement vectors f, whereas iterative methods for (3.4) need to solve $A\stackrel{^}{f}=f$ each time.

As already mentioned in [47, Sec. 3.1] an exact solution to (3.14) can only be found if , i.e., in case is an underdetermined system of equations. By [47, Lem. 3.1] this system has at least one solution, which is why we may choose the one with a minimal 2-norm. If $\text{rank}\left({\mathbit{A}}_{|{\mathcal{I}}_{2M}|}\right)=|{\mathcal{I}}_{2M}|$, then the system is consistent and the unique solution is given by the normal equations of the 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 ${\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}^{T}$ and $\overline{{\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}}$ are needed. Since fast multiplication with ${\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}^{T}$ and $\overline{{\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}}$ 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) is of complexity $\mathcal{O}\left(|{I}_{\mathbit{2}\mathbit{M}}|\mathrm{log}\left(|{I}_{\mathbit{2}\mathbit{M}}|\right)+N\right)$, where

$|I2M|=(2M)d=2dMd=2d|IM|.$

Thus, to receive exact quadrature weights wj ∈ ℂ, j = 1,…, N, via (3.15) we need to satisfy the full rank condition $\text{rank}\left({\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}\right)=|{\mathcal{I}}_{2M}|$. In case of a low-rank matrix ${A}_{|{\mathcal{I}}_{2M}|}\in {ℂ}^{N×|{\mathcal{I}}_{2M}|}$ for $|{ℐ}_{\mathbit{2}\mathbit{M}}|\le N$, we may still use (3.15) to obtain a least squares approximation to (3.14).

#### 3.1.2. Practical computation in the overdetermined setting $|{ℐ}_{\mathbit{2}\mathbit{M}}|>N$

In the setting $|{ℐ}_{\mathbit{2}\mathbit{M}}|>N$, we cannot expect to find an exact solution w to (3.14), since we have to deal with an overdetermined system possessing more conditions than variables. However, we still aimed to numerically find optimal density compensation factors wj∈ℂ, j = 1,…, N, by considering a least squares approximation to (3.14) that minimizes . In [60, Thm. 1.1.2] it was shown that every least squares solution satisfies the normal equations of the first kind

By means of definitions of ${\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}\in {\text{ℂ}}^{N×|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}$, cf. (2.7), and ${e}_{0}={\left({\delta }_{0,k}\right)}_{k\in {\mathcal{I}}_{2M}}$ we simplify the right-hand side via

$A|I2M|¯e0=(∑k∈I2Mδ0,ke-2πikxj)j=1N=1N.$

Since fast multiplication with ${\mathbit{A}}_{|{I}_{\mathbit{2}\mathbit{M}}|}^{T}$ and $\overline{{\mathbit{A}}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}}$ 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 $\mathcal{O}\left(|{\mathcal{I}}_{2M}|\mathrm{log}\left(|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|\right)+N\right)$ arithmetic operations. Note that the solution to (3.16) is only unique if the full rank condition $\text{rank}\left({A}_{|{\mathcal{I}}_{\mathbit{2}\mathbit{M}}|}\right)=N$ is satisfied, cf. [60, p. 7]. We remark that the computed weight matrix W = diag(w) can further be used in an iterative procedure as in Plonka et al. [[15, Alg. 7.27] to improve the approximation of $\stackrel{^}{f}$.

The previous considerations lead to the following algorithms.

ALGORITHM 3.6

Algorithm 3.6. Computation of the optimal density compensation factors.

ALGORITHM 3.7

Algorithm 3.7. iNFFT – density compensation approach.

### 3.2. Bandlimited functions

In some numerical examples, such as in MRI, we are concerned with bandlimited functions $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ instead of trigonometric polynomials $f\in {L}_{2}\left({𝕋}^{d}\right)$ in (2.8), cf. [1]. In the following we show that it is reasonable to consider the inversion problem (3.4) as well as the density compensation via (3.7) for bandlimited functions as well.

To this end, let $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ be a bandlimited function with bandwidth M, i.e., its (continuous) Fourier transform

is supported on ${\left[-\frac{M}{2},\frac{M}{2}\right)}^{d}$. Utilizing this fact, we have and thus by the Fourier inversion theorem [15, Thm. 2.10] the inverse Fourier transform of f can be written as

Analogous to (2.3), the approximation using equispaced quadrature points $k\in {\mathcal{I}}_{M}$ yields

such that evaluation at the given nonequispaced points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, j = 1,…, N, leads to

$f(xj)≈∑k∈IMf^(k)e2πikxj.$

By means of the definition (2.7) of the matrix $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ this can be written as $f\approx A\stackrel{^}{f}$, where we used the notation $\stackrel{^}{f}:={\left(\stackrel{^}{f}\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}$ in this setting. Thus, also for bandlimited functions f its evaluations at points xj 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 transform . However, in practical applications, such as MRI, this is only a hypothetical case, since f cannot be sampled on whole ℝd, cf. [1]. 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 ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{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\in {\mathcal{I}}_{M}$ can be approximated via

$f^(k)≈h~(k): =∑j=1Nwjf(xj)e-2πikxj, k∈IM.$

This is to say, equispaced samples of the Fourier transform of a bandlimited function may be approximated in the same form ${\left(\stackrel{^}{f}\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}=\stackrel{^}{f}\approx {h}^{\text{w}}:={A}^{*}Wf$ as in (3.6), where we used the notation ${h}^{\text{w}}:={\left(\stackrel{~}{h}\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}$ in this setting. Moreover, we extend this approximation onto the whole interval ${\left[-\frac{M}{2},\frac{M}{2}\right)}^{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 aimed to find a numerical method for computing suitable weights wj ∈ ℂ, 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.

Theorem 3.8. Let a bandwidth M ∈ ℕd, nonequispaced points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$ as well as quadrature weights wj ∈ ℂ, j = 1,…, N, be given. Then an exact reconstruction formula (3.21) for bandlimited functions $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ with bandwidth M, i.e.,

$∫ℝdf(x)dx=∑j=1Nwjf(xj)$

is exact for all bandlimited functions $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ with bandwidth M.

Proof: By (3.17) the assumption (3.22) can be written as

Especially, for v = 0 evaluation of (3.23) yields the assertion

$∫ℝdf(x)dx=f^(0)=∑j=1Nwjf(xj).$

However, in contrast to Theorem 3.2, this Theorem 3.8 does not yield an explicit condition for computing suitable weights wj∈ℂ, j = 1,…, N. To derive a numerical procedure anyway, we generalize the notion of an exact reconstruction $\stackrel{~}{h}$ of f and have a look at the theory of tempered distributions. To this end, let $\mathcal{S}$(ℝd) be the Schwartz space of rapidly decaying functions, cf. [15, Sec. 4.2.1]. The tempered Dirac distribution δ shall be defined by $〈\delta ,\phi 〉:={\int }_{{ℝ}^{d}}\phi \left(v\right)\delta \left(v\right)dv=\phi \left(0\right)$ for all φ ∈ $\mathcal{S}$(ℝd), cf. [15, Ex. 4.36]. For a slowly increasing function f:ℝd → ℂ satisfying $|f\left(x\right)|\le c{\left(1+||x|{|}_{2}\right)}^{n}$ almost everywhere with c > 0 and n ∈ ℕ0, the induced distribution Tf shall be defined by $〈{T}_{f},\phi 〉:=\underset{{ℝ}^{d}}{\int }\phi \left(x\right)f\left(x\right)dx$ for all φ ∈ $\mathcal{S}$(ℝd). For a detailed introduction to the topic we refer to [15, Sec. 4.2.1 and Sec. 4.3].

Then the following property can be shown.

Theorem 3.9. Let nonequispaced points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, j = 1, …, N, and quadrature weights wj ∈ ℂ be given. Furthermore, let Tf be the distribution induced by some bandlimited function $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ with bandwidth M. Then

with

implies

with the function $\stackrel{~}{h}$ defined in (3.21).

Proof: Using the definition of the function $\stackrel{~}{h}$ in (3.21) as well as the fact that the inversion formula (3.18) holds for all x ∈ ℝd, we have

Hence, by (3.24), this implies

Considering the property (3.26), we remark that this indeed states an exact reconstruction 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 aimed to use this result to compute suitable weights wj ∈ ℂ, j = 1,…, N, for bandlimited functions as well. To this end, suppose we have (3.14), i.e.,

$∑j=1Nwje2πikxj=δ0,k, k∈I2M.$

Then this yields

Having a look at Theorem 3.9, an exact reconstruction (3.26) is implied by (3.24), i.e.,

$φ(0)=〈δ,φ〉=〈Tξ,φ〉=∫ℝdφ(v)∑j=1Nwje2πivxjdv, φ∈S(ℝd).$

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\in {\mathcal{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 by

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 wj ∈ ℂ, j = 1,…, N, can also be derived by means of a quadrature formula applied to the property $〈\delta ,\stackrel{^}{\phi }〉=〈{T}_{\stackrel{~}{\xi }},\stackrel{^}{\phi }〉$ with $\stackrel{~}{\xi }$ 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 ∈ ℕd, nonequispaced points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, and quadrature weights wj ∈ ℂ, j = 1, …, N, be given. Then the following two statements are equivalent.

(i) For all φ ∈ $\mathcal{S}$(ℝd) we have $〈\delta ,\stackrel{^}{\phi }〉=〈{T}_{\stackrel{~}{\xi }},\stackrel{^}{\phi }〉$ with $\stackrel{~}{\xi }$ defined in (3.28).

(ii) We have 〈1, φ〉 = 〈Tψ, φ〉 for all φ ∈ $\mathcal{S}$(ℝd), where

$ψ(x): =∑j=1Nwj·|I2M|sinc(2Mπ(xj-x)), x∈ℝd,$

with the d-variate sinc function $\text{sinc}\left(x\right):={\prod }_{t=1}^{d}\text{sinc}\left({x}_{t}\right)$ and

$sinc(x): ={sinxx x∈ℝ\{0},1 x=0.$

Proof: By the definition $〈\stackrel{^}{T},\phi 〉=〈T,\stackrel{^}{\phi }〉$ of the Fourier transform of a tempered distribution T$\mathcal{S}$′(ℝd) we have $〈1,\phi 〉=〈\delta ,\stackrel{^}{\phi }〉$, cf. [15, 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

$∫[-M,M)de2πiv(xj-x)dv=|I2M|sinc(2Mπ(xj-x)),$

such that (3.29) shows the equality $〈{T}_{\stackrel{~}{\xi }},\stackrel{^}{\phi }〉=〈{T}_{\psi },\phi 〉$. 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 wj ∈ ℂ, j = 1,…, N. To this end, we have a closer look at

for all φ ∈ $\mathcal{S}$(ℝ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 $|{\mathcal{I}}_{2M}|$ of equispaced quadrature points ${y}_{\ell }:={\left(2M\right)}^{-1}\odot \ell$, $\ell \in {\mathcal{I}}_{2M}$, as in (3.27), i.e., we consider

$∑ℓ∈I2Mφ(yℓ)=∑ℓ∈I2Mφ(yℓ)∑j=1Nwj·|I2M|sinc(2Mπ(xj-yℓ)).$

In order that this applies for all φ ∈ $\mathcal{S}$(ℝd), we need to satisfy

i.e., one could also compute weights wj∈ℂ, 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.

Remark 3.12. Since we derived discretizations out of both statements of Theorem 3.11, we examine if also the two linear systems (3.14) and (3.31) are related. Considering the statements in Theorem 3.11 we notice that in some sense they are the Fourier transformed versions of each other. To this end, we need to Fourier transform one of the linear systems for better comparability. More precisely, we apply an iFFT of length $|{\mathcal{I}}_{2M}|$, cf. (2.4), to both sides of equation (3.31). Since the left side transforms to

$∑ℓ∈I2M1·e2πikyℓ=|I2M|·δ0,k, k∈I2M,$

we obtain the transformed system

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 ${\text{e}}^{2\pi \text{i}k{x}_{j}}$ and $\sum _{\ell \in {\mathcal{I}}_{2M}}\text{sinc}\left(2M\pi \left({x}_{j}-{y}_{\ell }\right)\right){\text{e}}^{2\pi \text{i}k{y}_{\ell }}$.

For this purpose, we consider the function f(t) = e2πitx, t ∈ [−M, M)d, for fixed x ∈ ℂd. By means of $\stackrel{~}{f}\left(t\right):=\sum _{k\in {ℤ}^{d}}f\left(t+2Mk\right)$ we extend it into a (2M)-periodic function. This periodized version then possesses the Fourier coefficients

cf. (2.2), i. e., the Fourier expansion of $\stackrel{~}{f}\left(t\right)$, t ∈ [−M, M)d, for fixed x is given by

$e2πitx=∑ℓ∈ℤde2πityℓsinc(2Mπ(x-yℓ)), x∈ℂd,$

cf. (2.1). Since $\stackrel{~}{f}\left(t\right)$ is continuous and piecewise differentiable, this Fourier series converges absolutely and uniformly, cf. [71, Ex. 1.22]. Thereby, we may consider the point evaluations at x = xj, j = 1,…, N, and $t=k\in {\mathcal{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}_{\ell }\in {ℝ}^{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 [4850], states that any bandlimited function $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ with maximum bandwidth M can be recovered from its uniform samples f(L−1), ∈ ℤd, with LM, L: = L·1d, 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 [5155] considered the regularized Shannon sampling formula with localized sampling

instead. Here, ${\phi }_{m}:{ℝ}^{d}\to \left[0,1\right]$ is a compactly supported window function with truncation parameter m ∈ ℕ\{1}, such that for φm with small support the direct evaluation of (3.35) is efficient, see [55] for the relation to the NFFT window functions.

On the other hand, in the one-dimensional setting a fast sinc transform was introduced in [56], which is based on the Clenshaw-Curtis quadrature

$sinc(Lπ(x-ℓL))≈∑k=0nwke-πiL(x-ℓL)zk$

using Chebyshev points ${z}_{k}=cos\left(\frac{k\pi }{n}\right)\in \left[-1,1\right]$, k = 0,…, n, and corresponding Clenshaw-Curtis weights wk > 0. Thereby, sums of the form

with uniform truncation parameter T ∈ 2ℕ, 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 [34, 57, 58] or [15, p. 394–397].

### 3.3. 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\in {L}_{2}\left({𝕋}^{d}\right)\cap C\left({𝕋}^{d}\right)$, and bandlimited functions $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ as well.

Theorem 3.14. Let p, q ∈ {1, 2, ∞} with $\frac{1}{p}+\frac{1}{q}=1$. For given d, N ∈ ℕ, M ∈ (2ℕ)d, and nonequispaced points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, j = 1,…, N, let $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ be the nonequispaced Fourier matrix in (2.7). Furthermore, assume we can compute density compensation factors $W=diag{\left({w}_{j}\right)}_{j=1}^{N}\in {ℂ}^{N×N}$ by means of Algorithm 3.6, such that

with small εk ∈ ℝ for all $k\in {\mathcal{I}}_{2M}$.

Then there exists an ε ≥ 0 such that the corresponding density compensation procedure with $W=\text{diag}{\left({w}_{j}\right)}_{j=1}^{N}$ satisfies the following error bounds.

(i) For any trigonometric polynomial $f\in {L}_{2}\left({T}^{d}\right)$ of degree M given in (2.8), we have

where $\stackrel{^}{f}:={\left({\stackrel{^}{f}}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$ are the coefficients given in (2.8).

(ii) For any 1-periodic function $f\in {L}_{2}\left({T}^{d}\right)\cap C\left({T}^{d}\right)$, we have

where $\stackrel{^}{f}:={\left({c}_{k}\left(f\right)\right)}_{k\in {\mathcal{I}}_{M}}$ are the first $|{\mathcal{I}}_{M}|$ coefficients given in (2.1) and pM is the best approximating trigonometric polynomial of degree M of f.

(iii) For any bandlimited function $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ with bandwidth M, we have

where $\stackrel{^}{f}:={\left(\stackrel{^}{f}\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}$ are the integer evaluations of (3.17) and Q in (3.45) is the pointwise quadrature error of the equispaced quadrature rule (3.19).

Proof: We start with some general considerations that are independent of the function f. By (3.36) we can find $\epsilon :=\underset{k\in {\mathcal{I}}_{M}}{max}|{\epsilon }_{k}|\ge 0$, such that |εk| ≤ ε, $k\in {\mathcal{I}}_{2M}$, and thereby

$|∑j=1Nwje2πikxj-δ0,k|≤ε, k∈I2M.$

Then for all $k,\ell \in {\mathcal{I}}_{M}$ with $\left(\ell -k\right)\in {\mathcal{I}}_{2M}$ this yields

$|[Er]k,ℓ|=|∑j=1Nwje2πi(ℓ-k)xj-δk,ℓ|≤ε,$

where ${E}_{\text{r}}:={A}^{*}WA-{I}_{|{\mathcal{I}}_{M}|}$. Hence, we have

and

Considering the approximation error of (3.6), it can be estimated by

Using $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ from (2.7) as well as $W=\text{diag}{\left({w}_{j}\right)}_{j=1}^{N}=diag\left(w\right)$, we have

$∥A*W∥∞=maxk∈IM∑j=1N|wj|·|e-2πikxj|≤∑j=1N|wj|·maxk∈IM1=∥w∥1,$

and

Hence, from (3.40) to (3.43) and ||·||2 ≤ ||·||F it follows that

for p ∈ {1, 2, ∞} with $\frac{1}{p}+\frac{1}{q}=1$. Now it merely remains to estimate $|A\stackrel{^}{f}-f{|}_{p}$ for the specific choice of f.

(i): Since a trigonometric polynomial (2.8) of degree M satisfies $A\stackrel{^}{f}=f$, the second error term in (3.44) vanishes and we obtain the assertion (3.37).

(ii): When considering a general 1-periodic function $f\in {L}_{2}\left({T}^{d}\right)\cap C\left({T}^{d}\right)$ in (2.1) we have

$r|[Af^-f]j|=|f(xj)-∑k∈IMck(f)e2πikxj|≤maxx∈𝕋d|∑k∈ℤd\IMck(f)e2πikx|=‖f-pM‖C(𝕋d),j=1,…,N,$

with the best approximating trigonometric polynomial pM of degree M of f. Thus, this yields and by (3.44) the assertion (3.38).

(iii): For a bandlimited function $f\in {L}_{1}\left({ℝ}^{d}\right)\cap {C}_{0}\left({ℝ}^{d}\right)$ with bandwidth M we may use the notation $\stackrel{^}{f}:={\left(\stackrel{^}{f}\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}$ as well as the inverse Fourier transform (3.18) to estimate

of the uniform quadrature rule (3.19). For detailed investigations of quadrature errors for bandlimited functions we refer to Kircheis et al. [56] and Gopal and Rokhlin [59]. Hence, we obtain $‖A\stackrel{^}{f}-f{‖}_{p}\le {N}^{1/p}||Q|{|}_{C\left({T}^{d}\right)}$ 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*WA being equal to identity ${I}_{|{\mathcal{I}}_{M}|}$. The following theorem shows that the error of the reconstruction (3.6) also affects the condition of the matrix A*WA.

Theorem 3.15. Let $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ from (2.7), $W=\text{diag}{\left({w}_{j}\right)}_{j=1}^{N}$ and ε ≥ 0 be given as in Theorem 3.14. If additionally $\epsilon |{\mathcal{I}}_{M}|<1$ is fulfilled, then we have

for the condition number ${\kappa }_{2}\left(X\right):=||X|{|}_{2}||{X}^{-1}|{|}_{2}$.

Proof: To estimate the condition number ${\kappa }_{2}\left({A}^{*}WA\right)$ we need to determine the norms $‖{A}^{*}WA{‖}_{2}$ and $‖{\left({A}^{*}WA\right)}^{-1}{‖}_{2}$. By (3.36) it is known that ${A}^{*}WA={I}_{|{\mathcal{I}}_{M}|}+\mathcal{E}$, where $\mathcal{E}:={\left({\epsilon }_{\ell -k}\right)}_{\ell ,k\in {\mathcal{I}}_{M}}$, and, therefore, we have

Moreover, it is known by the theory of Neumann series, cf. [70, Thm. 4.20], that if $‖{I}_{|{\mathcal{I}}_{M}|}-T{‖}_{2}<1$ holds for a matrix $T\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{M}|}$, then T is invertible and its inverse is given by

$T-1=∑n=0∞(I|IM|-T)n.$

Using this property for T = A*WA, we have

in case that $‖{I}_{|{\mathcal{I}}_{M}|}-{A}^{*}WA{‖}_{2}=||\mathcal{E}|{|}_{2}<1$. Hence, by (3.47) and (3.48) we obtain

In addition, we know that |εk| ≤ ε, $k\in {\mathcal{I}}_{2M}$, with some ε > 0 and, therefore,

In other words, the correctness of (3.48) is ensured if $\epsilon |{\mathcal{I}}_{M}|<1$. Since the spectral norm is a sub-multiplicative norm, (3.50) also implies $||{\mathcal{E}}^{n}|{|}_{2}\le ||\mathcal{E}|{|}_{2}^{n}\le {\left(\epsilon |{\mathcal{I}}_{M}|\right)}^{n}$. Consequently, we have

Thus, combining (3.49), (3.50), and (3.51) yields the assertion (3.46).

### 3.4. 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\in {ℝ}^{N×|{\mathcal{I}}_{M}|}$ from (3.52), the weight matrix $W=\text{diag}{\left({w}_{j}\right)}_{j=1}^{N}$ as well as the vectors $f={\left(f\left({x}_{j}\right)\right)}_{j=1}^{N}$ and $\stackrel{~}{f}={\left(f\left({M}^{-1}\odot \ell \right)\right)}_{l\in {\mathcal{I}}_{M}}$, the evaluation of (3.53) at equispaced points M−1, $\ell \in {\mathcal{I}}_{M}$, can be denoted as $\stackrel{~}{f}\approx {C}^{*}Wf$. Using the equispaced quadrature rule in (2.3), we find that evaluations of (3.17) at the uniform grid points $k\in {\mathcal{I}}_{M}$ can be approximated by (3.20) by means of a simple FFT. In matrix-vector notation, this can be written as $\stackrel{^}{f}\approx {\stackrel{~}{D}}^{*}{F}_{|{\mathcal{I}}_{M}|}^{*}\stackrel{~}{f}$ where $\stackrel{^}{f}={\left(\stackrel{^}{f}\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}$, ${F}_{|{\mathcal{I}}_{M}|}:={\left({\text{e}}^{2\pi \text{i}k\left({M}^{-1}\odot \ell \right)}\right)}_{\ell ,k\in {\mathcal{I}}_{M}}$, cf. (2.14), and $\stackrel{~}{D}:=\frac{1}{|{\mathcal{I}}_{M}|}{I}_{|{\mathcal{I}}_{M}|}$. Thus, all in all one obtains an approximation of the form $\stackrel{^}{f}\approx {\stackrel{~}{D}}^{*}{F}_{|{\mathcal{I}}_{M}|}^{*}{C}^{*}Wf$.

Here some of these approaches, cf. [2], 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.

#### 3.4.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-A\stackrel{^}{f}{‖}_{2}$. It is known (e.g., [60, 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 wj 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

where ||·||F denotes the Frobenius norm of a matrix. It was shown in Sedarat and Nishimura [61] that the solution to this least squares problem can be computed as

However, since a singular value decomposition is necessary for the calculation in (3.56), we obtain a high complexity of $\mathcal{O}\left({N}^{2}|{\mathcal{I}}_{M}|+|{\mathcal{I}}_{M}{|}^{3}\right)$. Therefore, we study some more sophisticated least squares approaches in the following.

#### 3.4.2. Density compensation using weighted normal equations of the first kind

It is known that every least squares solution to (3.4) satisfies the weighted normal equations of the first kind ${A}^{*}WA\stackrel{^}{f}={A}^{*}Wf$, see e.g., [60, Thm. 1.1.2]. As already mentioned in Corollary 3.4, we have an exact reconstruction formula (3.6) for all trigonometric polynomials (2.8) of degree M, if ${A}^{*}WA={I}_{|{\mathcal{I}}_{M}|}$ is fulfilled. Thus, we aimed to compute optional weights wj, j = 1,…, N, by considering the optimization problem

Analogous to Sedarat and Nishimura [61] 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 Er:= E·A, where E:= A*WA is the error matrix in (3.55).

In Rosenfeld et al. [62], it was shown that a solution W = diag(w) to (3.57) can be obtained by solving Sw = b with

However, since Sw = b is not separable for single wj, j = 1,…, N, computing these weights is of complexity $\mathcal{O}\left({N}^{3}\right)$. This is why the authors in Sedarat and Nishimura [61] restricted themselves to a maximal image size of 64 × 64 pixels, which corresponds to setting M = 64.

#### 3.4.3. Density compensation using weighted normal equations of the second kind

Another approach for density compensation factors is based on the weighted normal equations of the second kind

We recognize that by (3.59) we are given an exact approximation $\stackrel{^}{f}={A}^{*}Wf$ of the Fourier coefficients in (3.6) in case y = f, and thereby $A{A}^{*}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 El:= 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 Pipe and Menon [63]. Another version using a sinc transform evaluated at pointwise differences of the nonequispaced points instead of (3.52) was studied in [64, 65], where it was claimed that this approach coincides with the one in [63]. However, we remark that due to the sampling theorem of Shannon-Whittaker-Kotelnikov, see (3.34), applied to the function f(x) = sinc((xjx)), i.e.,

and its evaluation at x = xh, h = 1,…, N, this claim only holds asymptotically for $|{\mathcal{I}}_{M}|\to \infty$ 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 Greengard et al. [65] utilizes an approximation of the form fHWf, where the matrix H is defined as the system matrix of (3.4) evaluated at pointwise differences of the nonequispaced points, i.e.,

Since by (2.7) we have H = AA*, minimizing the approximation error

$‖H*Wf-f‖22=‖AA*Wf-f‖22≤‖AA*W-IN‖F2·‖f‖22,$

leads to the optimization problem (3.60) as well.

It was shown in [63] that the minimizer of (3.60) is given by

Since for fixed j the computation of ${\left[A{A}^{*}\right]}_{j,h}$, h = 1,…, N, is of complexity $\mathcal{O}\left(N|{\mathcal{I}}_{M}|\right)$, the weights (3.62) can be computed in $\mathcal{O}\left({N}^{2}|{\mathcal{I}}_{M}|\right)$ arithmetic operations. However, due to the explicit representation (3.61) the computation of ${\left[A{A}^{*}\right]}_{j,h}$, h = 1,…, N, for fixed j can be accelerated by means of the NFFT (see Algorithm 2.1). Then, this step takes $\mathcal{O}\left(|{\mathcal{I}}_{M}|log\left(|{\mathcal{I}}_{M}|\right)+N\right)$ arithmetic operations and the overall complexity is given by $\mathcal{O}\left(N·|{\mathcal{I}}_{M}|log\left(|{\mathcal{I}}_{M}|\right)+{N}^{2}\right)$.

As mentioned in Pipe and Menon [63] 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 $A{A}^{*}W={I}_{N}$ as

By means of (2.7) this can be written as $A{A}^{*}w={0}_{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 $A{A}^{*}w={0}_{N}$ can be computed iteratively with arithmetic complexity $\mathcal{O}\left(|{\mathcal{I}}_{M}|log\left(|{\mathcal{I}}_{M}|\right)+N\right)$.

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 ∈ ℂN, i.e., by ${A}^{*}=\overline{{A}^{T}}$ we have ${\left({\delta }_{0,k}\right)}_{k\in {\mathcal{I}}_{M}}=\overline{{A}^{T}w}={A}^{*}\overline{w}$. Then multiplication with $A\in {ℂ}^{N×|{\mathcal{I}}_{M}|}$ in (2.7) yields

$AA*w¯=A·(δ0,k)k∈IM=(∑k∈IMδ0,k·e2πikxj)j=1N=1N.$

In other words, an exact solution w to the linear system (3.10) implies that the conjugate complex weights $\overline{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 reconstruction in (3.6) for trigonometric polynomials (2.8) with maximum degree M.

## 4. 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\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ 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. [33]. To this end, we first 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 aimed 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\approx {I}_{|{\mathcal{I}}_{M}|}$, since then we can simply compute an approximation of the Fourier coefficients by means of $Xf=XA\stackrel{^}{f}\approx \stackrel{^}{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 $|{\mathcal{I}}_{M}|\le 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\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{M}|}$, $F\in {ℂ}^{|{\mathcal{I}}_{{M}_{\sigma }}|×|{\mathcal{I}}_{M}|}$, and $B\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ defined in (2.13), (2.14), and (2.15). Combining both ingredients we aimed for an approximation of the form ${D}^{*}{F}^{*}{B}^{*}A\approx {I}_{|{\mathcal{I}}_{M}|}$. To achieve an approximation like this, we aimed 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}^{*}WA={I}_{|{\mathcal{I}}_{M}|}$. Since the reconstruction shall be realized efficiently by means of an adjoint NFFT, one rather studies ${D}^{*}{F}^{*}{B}^{*}WA\approx {I}_{|{\mathcal{I}}_{M}|}$. Using the definition $\stackrel{~}{B}:={W}^{*}B$ as in Remark 3.1, we end up with an approximation of the form ${D}^{*}{F}^{*}{\stackrel{~}{B}}^{*}A\approx {I}_{|{\mathcal{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.

Let $\stackrel{~}{B}$ denote such a modified matrix. By defining $\stackrel{~}{h}:={D}^{*}{F}^{*}{\stackrel{~}{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 Nieslony and Steidl [66]. 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 ${F}^{*}F=|{\mathcal{I}}_{{M}_{\sigma }}|{I}_{|{\mathcal{I}}_{M}|}$ and $D\in {ℝ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{M}|}$ is diagonal by (2.13), we have $\frac{1}{|{\mathcal{I}}_{{M}_{\sigma }}|}{D}^{-1}{F}^{*}FD={I}_{|{\mathcal{I}}_{M}|}$. 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 ∈ ℝk×n and the definition of the Euclidean norm of a vector y ∈ ℝn, we obtain for zj being the columns of Z ∈ ℝk×n that

Since we aimed 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 of $\stackrel{~}{B}$ 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\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$. Thus, we have

where ${\stackrel{~}{b}}_{\ell }\in {ℝ}^{|{\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)|}$ denotes the vectors of the nonzeros of each column of $\stackrel{~}{B}$,

are the corresponding submatrices of ${A}^{*}\in {ℂ}^{|{\mathcal{I}}_{M}|×N}$, cf. (2.7), and ${f}_{\ell }\in {ℂ}^{|{\mathcal{I}}_{M}|}$ are the columns of ${F}^{*}\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{{M}_{\sigma }}|}$, cf. (2.14). Hence, we receive the equivalent optimization problems

Thus, if the matrix (4.9) has full column rank, the solution of the least squares problem (4.10) can be computed by means of the pseudoinverse ${H}_{\ell }^{†}$ as

Having these vectors ${b}_{\ell }^{\text{opt}}$ we compose the optimized matrix Bopt, observing that ${b}_{\ell }^{\text{opt}}$ only consist of the nonzero entries of Bopt. Then, the approximation of the Fourier coefficients is given by

In other words, this approach yields an inverse NFFT by modifying the adjoint NFFT.

Remark 4.2. To achieve an efficient algorithm we now have a closer look at the computation scheme (4.11). We start with the computation of the matrix ${H}_{\ell }^{*}{H}_{\ell }$. By introducing the d-dimensional Dirichlet kernel

the matrix ${H}_{\ell }^{*}{H}_{\ell }$ in (4.11) can explicitly be stated via

i.e., for given index set ${\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)$ the matrix ${H}_{\ell }^{*}{H}_{\ell }$ can be determined in $\mathcal{O}\left(|{\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)|\right)$ operations. Considering the right-hand side of (4.11), by definitions (4.9), (2.13), and (2.14), we have

Thus, since $\frac{1}{|{\mathcal{I}}_{{M}_{\sigma }}|}{D}^{-1}=\text{diag}{\left(ŵ\left(k\right)\right)}_{k\in {\mathcal{I}}_{M}}$, the computation of v involves neither multiplication with nor division by the (possibly) huge number $|{\mathcal{I}}_{{M}_{\sigma }}|$ and is, therefore, numerically stable.

This leads to the following algorithm.

ALGORITHM 4.3

Algorithm 4.3. Optimization of the sparse matrix B.

Note that a general statement about the dimensions of ${H}_{\ell }\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)|}$ is not possible, since the size of the set ${\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)$ 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\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ in Figure 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\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ contains at most (2m+1)d entries, each column contains $\frac{N}{|{\mathcal{I}}_{{M}_{\sigma }}|}{\left(2m+1\right)}^{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 $|{\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)|$ is a small constant compared to $|{\mathcal{I}}_{M}|$, such that Algorithm 4.3 ends up with total arithmetic costs of approximately $\mathcal{O}\left(|{\mathcal{I}}_{M}|\right)$.

FIGURE 1

Figure 1. Nonzero entries of the matrix $B\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ for several choices of the points ${x}_{j}\in {𝕋}^{d}$, j = 1,…, N, with d = 1, Mσ = M = 16, N = 2M, and m = 2. (A) Equispaced points. (B) Jittered points. (C) Chebyshev points. (D) Random points.

In conclusion, our approach for an inverse NFFT can be summarized as follows.

ALGORITHM 4.4

Algorithm 4.4. iNFFT – optimization approach.

Theorem 4.5. Let ${B}_{\text{opt}}\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ be the optimized matrix computed by means of Algorithm 4.3 and let ${h}_{\text{opt}}={D}^{*}{F}^{*}{B}_{\text{opt}}^{*}f\in {ℂ}^{|{\mathcal{I}}_{M}|}$ be the corresponding approximation of $\stackrel{^}{f}$ computed by means of Algorithm 4.4. Further assume that each column ${b}_{\ell }^{\text{opt}}\in {ℝ}^{|{\mathcal{I}}_{{M}_{\sigma },m}|}$ of ${B}_{\text{opt}}\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ as solution to (4.10) possesses a small residual

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

$‖A*Bopt-1‖IMσ|D-1F*|F2=∑ℓ∈IMσ‖Hℓbℓopt-1|IMσ|D-1fℓ‖22,$

where ${b}_{\ell }^{\text{opt}}\in {ℝ}^{|{\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)|}$ are the nonzeros of the columns of Bopt, ${H}_{\ell }\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)|}$ in (4.9) are the corresponding submatrices of ${A}^{*}\in {ℂ}^{|{\mathcal{I}}_{M}|×N}$, cf. (2.7), and ${f}_{\ell }\in {ℂ}^{|{\mathcal{I}}_{M}|}$ are the columns of ${F}^{*}\in {ℂ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{{M}_{\sigma }}|}$, cf. (2.14). Since ${b}_{\ell }^{\text{opt}}\in {ℝ}^{|{\mathcal{I}}_{{M}_{\sigma },m}|}$ as solutions to the least squares problems (4.10) satisfy (4.15), we can find $\epsilon :=\underset{\ell \in {\mathcal{I}}_{{M}_{\sigma }}}{max}{\epsilon }_{\ell }\ge 0$, such that ε ≤ ε, $\ell \in {\mathcal{I}}_{{M}_{\sigma }}$, and, therefore,

$∑ℓ∈IMσ‖Hℓbℓopt-1|IMσ|D-1fℓ‖22≤∑ℓ∈IMσεℓ≤ε|IMσ|.$

Therefore, we may write (4.19) as

Thus, it remains to estimate the Frobenius norm $‖FD{‖}_{F}^{2}$. By the definitions of the Frobenius norm and the trace tr(Z) of a matrix Z, it is clear that $||Z|{|}_{F}^{2}=\text{tr}\left({Z}^{*}Z\right)$. Since by (2.14) we have that ${F}^{*}F=|{\mathcal{I}}_{{M}_{\sigma }}|{I}_{|{\mathcal{I}}_{M}|}$, this yields

Using the definition (2.13) of the diagonal matrix $D\in {ℝ}^{|{\mathcal{I}}_{M}|×|{\mathcal{I}}_{M}|}$, we obtain

Then combination of (4.20), (4.21), and (4.22) implies

$‖D*F*Bopt*A-I|IM|‖F2≤ε∑k∈IM1w^(k)2,$

such that (4.18) yields the assertion (4.16).

Since it is known that 0 ≤ ŵ(k) ≤ 1, $k\in {\mathcal{I}}_{M}$, for suitable window functions of the NFFT, cf. [67], we have

$1≤1w^(k)≤1w^(k)2$

and, therefore.

$∑k∈IM1w^(k)2≥∑k∈IM1=|IM|.$

Hence, the smallest constant is achieved in (4.16) when ŵ(k) = 1, $k\in {\mathcal{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 $|{\mathcal{I}}_{M}|\le 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 the second kind

$A*Ay=h, f=Ay.$

Incorporating the matrix decomposition of the NFFT, cf. Section 2.1, we recognize that a modification of the matrix $B\in {ℝ}^{N×|{\mathcal{I}}_{{M}_{\sigma }}|}$ such that ${A}^{*}BFD\approx {I}_{|{\mathcal{I}}_{M}|}$ implies yh and hence fBFDh. 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.

## 5. Numerics

In conclusion, we have a look at some numerical examples. In addition to 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.,

t = 1,…, N1, j = 1,…, N2, 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 2A. In addition, we also consider the random grid ${x}_{t,j}:=\frac{1}{2}{\left({\eta }_{1},{\eta }_{2}\right)}^{T},$ see Figure 2B.

(ii) Moreover, we examine grids of polar kind, as mentioned in Fenn et al. [68]. For R, T ∈ 2ℕ, the points of the polar grid are given by a signed radius ${r}_{j}:=\frac{j}{R}\in \left[-\frac{1}{2},\frac{1}{2}\right)$ and an angle ${\theta }_{t}:=\frac{\pi t}{T}\in \left[-\frac{\pi }{2},\frac{\pi }{2}\right)$ as

Since it is known that the inversion problem is ill-conditioned for this grid we consider a modification, the modified polar grid

i.e., we added more concentric circles and excluded the points outside the unit square, see Figure 3A. Another sampling scheme that 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 3B. There we distinguish two sets of points, i.e.,

(iii) Another modification of these polar grids was introduced in Helou et al. [69], 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

The golden angle linogram grid is given by

$xt,j: ={(2j+12R,2j+12Rtan(θt-π4))T:θt∈[0,π2)(-2j+12Rcot(θt-π4),2j+12R)T:θt∈[-π2,0)}, j∈IR,$

with θt in (5.5), as illustrated in Figure 3C.

FIGURE 2

Figure 2. Exemplary randomized grids of size N1 = N2 = 12. (A) Jittered grid. (B) Random grid.

FIGURE 3

Figure 3. Polar grids of size R = 12 and T = 2R. (A) Polar (blue) and modified polar (red) grid. (B) Linogram/pseudo-polar grid. (C) Golden angle linogram grid.

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. First, 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 coefficients , $k\in {\mathcal{I}}_{M}$. In this test, we consider several d ∈ {1, 2, 3}. For the corresponding function evaluations of (2.8) at given points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, j = 1,…, N, we test how well these Fourier coefficients can be approximated. More precisely, we consider the estimate ${\stackrel{~}{h}}^{\text{w}}={D}^{*}{F}^{*}{B}^{*}Wf$, cf. (3.7), with the matrix $W=diag{\left({w}_{j}\right)}_{j=1}^{N}$ of density compensation factors computed by means of Algorithm 3.6, i.e., by (3.15), in case , or by (3.16), if , and compute the relative errors

By (3.37), it is known that

$‖f^-A*Wf‖p‖f^‖p≤|IM|ε, p∈{2,∞},$

with the residual , cf. (3.36).

In our experiment, we use random points xj with ${N}_{t}={2}^{9-d}$, t = 1,…, d, cf. Figure 2B, and, for several problem sizes M = M·1d, M = 2c with c = 1,…, 11−d, we choose random Fourier coefficients , $k\in {\mathcal{I}}_{M}$. Afterward, 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 4. It can clearly be seen that $|{\mathcal{I}}_{2M}|, i.e., as long as $M<\frac{{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 $|{\mathcal{I}}_{2M}|>N$ the least squares approximation via (3.16) does not yield good results anymore.

FIGURE 4

Figure 4. Relative errors (5.6) of the reconstruction of the Fourier coefficients of a trigonometric polynomial (2.8) with given , $k\in {\mathcal{I}}_{M}$, computed via the density compensation method from Algorithm 3.7, for random grids with ${N}_{t}={2}^{9-d}$, t = 1,…, d, and M = M·1d, M = 2c with c = 1,…, 11−d. (A) d = 1. (B) d = 2. (C) d = 3.

Example 5.3. To study the quality of our optimization method in Section 4, we consider the Frobenius norms

where B denotes the original matrix from the NFFT in (2.15) and Bopt is the optimized matrix generated by Algorithm 4.3. For the original matrix B we utilize the common B-Spline window function

with the centered B-Spline of order 2m, cf. [15, p. 388]. The optimized matrix Bopt shall be computed by means of the B-Spline (5.8) as well as the Dirichlet window function (4.17), which is the optimal window by Theorem 4.5.

Due to memory limitations in the computation of the Frobenius norms (5.7), we have to settle for very small problems, which however show the functionality of Algorithm 4.3. For this reason, we consider d = 2 and choose M = (12, 12)T as well as ${N}_{1}={N}_{2}=R={2}^{\mu }$, μ ∈ {2,…, 7}, and T = 2R for the grids mentioned in Remark 5.1. In other words, we test Algorithm 4.3 in the underdetermined setting $|{\mathcal{I}}_{M}|>N$ as well as for the overdetermined setting $|{\mathcal{I}}_{M}|\le N$.

Having a look at the results for the grids in Remark 5.1, it becomes apparent that they separate into two groups. Figure 5A 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 of these grids are depicted in Figure 5B exemplarily for the modified polar grid (5.3). Moreover, it can be seen that our optimization procedure in Algorithm 4.3 is the most effective in the overdetermined setting $|{\mathcal{I}}_{M}|\le N$.

FIGURE 5

Figure 5. Frobenius norms (5.7) of the original matrix B (violet) and the optimized matrix Bopt generated by Algorithm 4.3 using the B-Spline wB (orange) as well as the Dirichlet window wD (cyan) with R = 2μ, μ ∈ {2,…, 7}, and T = 2R as well as M = (12, 12)T, m ∈ {2, 4}, and σ ∈ {1, 2}. (A) Polar grid (5.2). (B) Modified polar grid (5.3).

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}_{\ell }^{*}{H}_{\ell }$, 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 3A, there are no grid points in the corners of the unit square. Therefore, some of the index sets ${\mathcal{I}}_{{M}_{\sigma },m}\left(\ell \right)$, 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 wj 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 Averbuch et al. [32] and Kircheis and Potts [33] we have a look at the reconstruction of the Shepp-Logan phantom, see Figure 6A. Here, we treat the phantom data as given Fourier coefficients $\stackrel{^}{f}:={\left({\stackrel{^}{f}}_{k}\right)}_{k\in {\mathcal{I}}_{M}}$ of a trigonometric polynomial (2.8). For given points ${\mathbit{x}}_{j}\in {\left[-\frac{1}{2},\frac{1}{2}\right)}^{d}$, 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.

FIGURE 6

Figure 6. Reconstruction of the Shepp-Logan phantom of size 1024 × 1024 (top) via the density compensation method from Section 3.1 using Voronoi weights and Algorithm 3.7 with weights computed by (3.16) compared to Algorithm 4.4 for the linogram grid (5.4) of size R = M = 1, 024, T = 2R; as well as a detailed view of the 832nd row each (bottom). (A) Original phantom. (B) Voronoi weights with e2=5.3040e-01. (C) Algorithm 3.7 with e2=5.0585e-01. (D) Algorithm 4.4 with e2=2.2737e-03.

In a first experiment, we test the inversion methods from Sections 3 and 4 as in Averbuch et al. [32] for increasing input sizes. To this end, we choose M = (M, M)T, M = 2c with c = 3,…, 10, and linogram grids (5.4) of size R = 2M, T = 2R, i.e., we consider the setting $|{\mathcal{I}}_{2M}|. 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 Bopt, as well as the time needed for the reconstruction, i.e., the corresponding adjoint NFFT, see Algorithms 3.7, 4.4. Moreover, for the reconstruction $\stackrel{~}{h}\in \left\{{\stackrel{~}{h}}^{\text{w}},{h}_{\text{opt}}\right\}$, cf. (3.7) and (4.12), we consider the relative errors

$e2: =‖h~-f^‖2‖f^$