$\displaystyle=\big{\|}\big{(}\boldsymbol{A}^{*}\boldsymbol{\tilde{B}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\big{)}\boldsymbol{F}\boldsymbol{D}\big{\|}_{\mathrm{F}}$
$\displaystyle\leq\big{\|}\boldsymbol{A}^{*}\boldsymbol{\tilde{B}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\big{\|}_{\mathrm{F}}\,\big{\|}\boldsymbol{F}\boldsymbol{D}\big{\|}_{\mathrm{F}}.$
(4.4)

Hence, we consider the optimization problem

$\underset{\boldsymbol{\tilde{B}}\in\mathbb{R}^{N\times{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}}\colon\boldsymbol{\tilde{B}}\,(2m+1)^{d}\text{-sparse
}}{\text{Minimize }}\
\big{\|}\boldsymbol{A}^{*}\boldsymbol{\tilde{B}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\big{\|}_{\mathrm{F}}^{2}.$
(4.5)

Based on the definition of the Frobenius norm of a matrix
$\boldsymbol{Z}\in\mathbb{R}^{k\times n}$ and the definition of the Euclidean
norm of a vector $\boldsymbol{y}\in\mathbb{R}^{n}$, we obtain for
$\boldsymbol{z}_{j}$ being the columns of
$\boldsymbol{Z}\in\mathbb{R}^{k\times n}$ that

$\|\boldsymbol{Z}\|_{F}^{2}=\sum_{i=1}^{k}\sum_{j=1}^{n}|z_{ij}|^{2}=\sum_{j=1}^{n}\|\boldsymbol{z}_{j}\|_{2}^{2}.$
(4.6)

Since we aim to preserve the property that $\boldsymbol{B}$ is a
$(2m+1)^{d}$-sparse matrix, we rewrite the norm in (4.5) by (4.6) in terms of
the columns of $\boldsymbol{\tilde{B}}$ considering only the nonzero entries
of each column. To this end, analogously to (2.12) we define the index set

$\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell}):=\left\\{j\in\left\\{1,\dots,N\right\\}:\exists\,\boldsymbol{z}\in\mathbb{Z}^{d}\
\text{with}-m\boldsymbol{1}\leq\boldsymbol{M}_{\boldsymbol{\sigma}}\odot\boldsymbol{x}_{j}-\boldsymbol{\ell}+\boldsymbol{z}\leq
m\boldsymbol{1}\right\\}$ (4.7)

of the nonzero entries of the $\boldsymbol{\ell}$-th column of
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$.
Thus, we have

$\big{\|}\boldsymbol{A}^{*}\boldsymbol{\tilde{B}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\big{\|}_{\mathrm{F}}^{2}=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}\big{\|}\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{\tilde{b}}_{\boldsymbol{\ell}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}}\big{\|}_{2}^{2},$
(4.8)

where
$\boldsymbol{\tilde{b}}_{\boldsymbol{\ell}}\in\mathbb{R}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|}$
denotes the vectors of the nonzeros of each column of
$\boldsymbol{\tilde{B}}$,

$\boldsymbol{H}_{\boldsymbol{\ell}}\coloneqq\left(\mathrm{e}^{-2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}_{j}}\right)_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}},\,j\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|}$
(4.9)

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

$\underset{\boldsymbol{\tilde{b}}_{\boldsymbol{\ell}}\in\mathbb{R}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|}}{\text{Minimize
}}\
\big{\|}\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{\tilde{b}}_{\boldsymbol{\ell}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}}\big{\|}_{2}^{2},\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}.$
(4.10)

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
$\boldsymbol{H}_{\boldsymbol{\ell}}^{\dagger}$ as

$\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}\coloneqq\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\left(\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{H}_{\boldsymbol{\ell}}\right)^{-1}\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}},\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}.$
(4.11)

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

$\boldsymbol{\hat{f}}\approx\boldsymbol{h}_{\mathrm{opt}}\coloneqq\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{f}.$
(4.12)

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
$\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{H}_{\boldsymbol{\ell}}$. By
introducing the $d$-dimensional Dirichlet kernel

$D_{\boldsymbol{M}}(\boldsymbol{x})\coloneqq\sum_{k_{1}=-\frac{M}{2}+1}^{\frac{M}{2}-1}\dots\sum_{k_{d}=-\frac{M}{2}+1}^{\frac{M}{2}-1}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}=\prod_{t=1}^{d}D_{\frac{M}{2}-1}(x_{t})=\prod_{t=1}^{d}\frac{\sin((M-1)\pi
x_{t})}{\sin(\pi x_{t})},$

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

$\displaystyle\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{H}_{\boldsymbol{\ell}}$
$\displaystyle=\Bigg{[}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}(\boldsymbol{x}_{h}-\boldsymbol{x}_{j})}\Bigg{]}_{h,j\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})}$
$\displaystyle=\left[\prod_{t=1}^{d}\left(D_{\frac{M}{2}-1}({x_{h}}_{t}-{x_{j}}_{t})+\mathrm{e}^{-M\pi\mathrm{i}({x_{h}}_{t}-{x_{j}}_{t})}\right)\right]_{h,j\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})},$
(4.13)

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

$\displaystyle\boldsymbol{v}_{\boldsymbol{\ell}}\coloneqq\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}}$
$\displaystyle=\Bigg{(}\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\hat{w}(\boldsymbol{k})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\left(\boldsymbol{x}_{j}-\boldsymbol{M}_{\boldsymbol{\sigma}}^{-1}\odot\,\boldsymbol{\ell}\right)}\Bigg{)}_{j\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})},\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}.$
(4.14)

Thus, since
$\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}=\mathrm{diag}(\hat{w}(\boldsymbol{k}))_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$,
the computation of $\boldsymbol{v}_{\boldsymbol{\ell}}$ involves neither
multiplication with nor division by the (possibly) huge number
$|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|$ and is therefore
numerically stable.

This leads to the following algorithm.

######  Algorithm 4.3 (Optimization of the sparse matrix $\boldsymbol{B}$).

For $d,N\in\mathbb{N}$ let $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N,$ be given points. Further let $\boldsymbol{M}\coloneqq
M\cdot\boldsymbol{1}_{d}$ with $M\in 2\mathbb{N}$, an oversampling factor
$\sigma\geq 1$ with $2\mathbb{N}\ni M_{\sigma}\coloneqq 2\lceil\lceil\sigma
M\rceil/2\rceil$ and $\boldsymbol{M}_{\boldsymbol{\sigma}}\coloneqq
M_{\sigma}\cdot\boldsymbol{1}_{d}$ as well as a truncation parameter $m\ll
M_{\sigma}$ be given.

  1. 1.

For $\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}$:

     * Determine the index set $\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})$, cf. (4.7). $\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|)$

     * Compute the right side $\boldsymbol{v}_{\boldsymbol{\ell}}$ via (4.14). $\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|)$

     * Determine $\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{H}_{\boldsymbol{\ell}}$ via (4.2). $\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|^{2})$

     * Solve $\left(\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{H}_{\boldsymbol{\ell}}\right)\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}=\boldsymbol{v}_{\boldsymbol{\ell}}$, i. e., compute $\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}$, cf. (4.11). $\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|^{3})$

  2. 2.

Compose
$\boldsymbol{B}_{\mathrm{opt}}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
columnwise of the
$\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}\in\mathbb{R}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|}$.
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|)$

Output: optimized matrix $\boldsymbol{B}_{\mathrm{opt}}$ Complexity:
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\cdot|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|^{3})$

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

(a) Equispaced points

(b) Jittered points

(c) Chebyshev points

(d) Random points

Figure 4.1:  Nonzero entries of the matrix
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
for several choices of the points $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N,$ with $d=1$, $M_{\sigma}=M=16$, $N=2M$ and $m=2$.

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

######  Algorithm 4.4 (iNFFT – optimization approach).

For $d,N\in\mathbb{N}$ let $\boldsymbol{x}_{j}\in\mathbb{T}^{d}$,
$j=1,\dots,N,$ be given points as well as $\boldsymbol{f}\in\mathbb{C}^{N}$.
Further let $\boldsymbol{M}\coloneqq M\cdot\boldsymbol{1}_{d}$ with $M\in
2\mathbb{N}$, an oversampling factor $\sigma\geq 1$ with $2\mathbb{N}\ni
M_{\sigma}\coloneqq 2\lceil\lceil\sigma M\rceil/2\rceil$ and
$\boldsymbol{M}_{\boldsymbol{\sigma}}\coloneqq
M_{\sigma}\cdot\boldsymbol{1}_{d}$ as well as a truncation parameter $m\ll
M_{\sigma}$ be given.

  1. 0.

Precompute the optimal sparse matrix $\boldsymbol{B}_{\mathrm{opt}}$ using
Algorithm 4.3.

  2. 1.

Compute
$\boldsymbol{h}_{\mathrm{opt}}\coloneqq\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{f}$,
cf. (4.12), by means of a modified adjoint NFFT.

Output:
$\boldsymbol{h}_{\mathrm{opt}}\approx\boldsymbol{\hat{f}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}$,
cf. (3.4). Complexity:
$\mathcal{O}(|\mathcal{I}_{\boldsymbol{M}}|\log(|\mathcal{I}_{\boldsymbol{M}}|)+N)$

######  Theorem 4.5.

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

$\displaystyle\big{\|}\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}}\big{\|}_{2}^{2}=\varepsilon_{\boldsymbol{\ell}}\geq
0,\quad\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}.$
(4.15)

Then there exists an $\varepsilon\geq 0$ such that

$\displaystyle\big{\|}\boldsymbol{h}_{\mathrm{opt}}-\boldsymbol{\hat{f}}\big{\|}_{2}^{2}$
$\displaystyle\leq\varepsilon\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\frac{1}{\hat{w}(\boldsymbol{k})^{2}}\cdot\big{\|}\boldsymbol{\hat{f}}\big{\|}_{2}^{2}.$
(4.16)

Moreover, the (asymmetric) Dirichlet kernel

$\displaystyle
w_{\mathrm{D}}\coloneqq\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{k}\boldsymbol{x}}=\prod_{t=1}^{d}\left(D_{\frac{M}{2}-1}(x_{t})+\mathrm{e}^{-M\pi\mathrm{i}x_{t}}\right)$
(4.17)

is the optimal window function for the inverse NFFT in Algorithm 4.4.

###### Proof.

As in (4) the approximation error can be estimated by

$\displaystyle\big{\|}\boldsymbol{h}_{\mathrm{opt}}-\boldsymbol{\hat{f}}\big{\|}_{2}^{2}$
$\displaystyle=\big{\|}\big{(}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{)}\boldsymbol{\hat{f}}\big{\|}_{2}^{2}\leq\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}\,\big{\|}\boldsymbol{\hat{f}}\big{\|}_{2}^{2}.$
(4.18)

Using the same arguments as for (4.3) and (4) we proceed with

$\displaystyle\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}$
$\displaystyle=\big{\|}\boldsymbol{A}^{*}\boldsymbol{B}_{\mathrm{opt}}\boldsymbol{F}\boldsymbol{D}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}$
$\displaystyle\leq\big{\|}\boldsymbol{A}^{*}\boldsymbol{B}_{\mathrm{opt}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\big{\|}_{\mathrm{F}}^{2}\,\big{\|}\boldsymbol{F}\boldsymbol{D}\big{\|}_{\mathrm{F}}^{2}.$
(4.19)

To estimate the first Frobenius norm in (4), we rewrite it analogously to
(4.8) columnwise as

$\displaystyle\big{\|}\boldsymbol{A}^{*}\boldsymbol{B}_{\mathrm{opt}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{F}^{*}\big{\|}_{\mathrm{F}}^{2}=\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}\big{\|}\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}}\big{\|}_{2}^{2},$

where
$\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}\in\mathbb{R}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|}$
are the nonzeros of the columns of $\boldsymbol{B}_{\mathrm{opt}}$,
$\boldsymbol{H}_{\boldsymbol{\ell}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})|}$
in (4.9) are the corresponding submatrices of
$\boldsymbol{A}^{*}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times N}$,
cf. (2.7), and
$\boldsymbol{f}_{\boldsymbol{\ell}}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|}$
are the columns of
$\boldsymbol{F}^{*}\in\mathbb{C}^{|\mathcal{I}_{\boldsymbol{M}}|\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$,
cf. (2.14). Since
$\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}\in\mathbb{R}^{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}|}$
as solutions to the least squares problems (4.10) satisfy (4.15), we can find
$\varepsilon\coloneqq\max_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}\varepsilon_{\boldsymbol{\ell}}\geq
0$, such that $\varepsilon_{\boldsymbol{\ell}}\leq\varepsilon$,
$\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}$, and
thereby

$\displaystyle\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}\big{\|}\boldsymbol{H}_{\boldsymbol{\ell}}\boldsymbol{b}_{\boldsymbol{\ell}}^{\mathrm{opt}}-\tfrac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}\boldsymbol{D}^{-1}\boldsymbol{f}_{\boldsymbol{\ell}}\big{\|}_{2}^{2}\leq\sum_{\boldsymbol{\ell}\in\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}}\varepsilon_{\boldsymbol{\ell}}\leq\varepsilon\,|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|.$

Therefore, we may write (4) as

$\displaystyle\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}$
$\displaystyle\leq\varepsilon\,|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\cdot\,\big{\|}\boldsymbol{F}\boldsymbol{D}\big{\|}_{\mathrm{F}}^{2}.$
(4.20)

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

$\displaystyle\big{\|}\boldsymbol{F}\boldsymbol{D}\big{\|}_{\mathrm{F}}^{2}=\mathrm{tr}(\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{F}\boldsymbol{D})=|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\cdot\mathrm{tr}(\boldsymbol{D}^{*}\boldsymbol{D})=|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|\cdot\big{\|}\boldsymbol{D}\big{\|}_{\mathrm{F}}^{2}.$
(4.21)

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

$\displaystyle\big{\|}\boldsymbol{D}\big{\|}_{\mathrm{F}}^{2}=\frac{1}{|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|^{2}}\,\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\frac{1}{\hat{w}(\boldsymbol{k})^{2}}.$
(4.22)

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

$\displaystyle\big{\|}\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{A}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\big{\|}_{\mathrm{F}}^{2}$
$\displaystyle\leq\varepsilon\,\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\frac{1}{\hat{w}(\boldsymbol{k})^{2}},$

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

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

$\displaystyle
1\leq\frac{1}{\hat{w}(\boldsymbol{k})}\leq\frac{1}{\hat{w}(\boldsymbol{k})^{2}}$

and therefore

$\displaystyle\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}\frac{1}{\hat{w}(\boldsymbol{k})^{2}}\geq\sum_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}1=|\mathcal{I}_{\boldsymbol{M}}|.$

Hence, the smallest constant is achieved in (4.16) when
$\hat{w}(\boldsymbol{k})=1$, $\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{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.39)
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}_{\boldsymbol{M}}|\leq 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

$\boldsymbol{A}^{*}\boldsymbol{A}\boldsymbol{y}=\boldsymbol{h},\quad\boldsymbol{f}=\boldsymbol{A}\boldsymbol{y}.$

Incorporating the matrix decomposition of the NFFT, cf. Section 2.1, we
recognize that a modification of the matrix
$\boldsymbol{B}\in\mathbb{R}^{N\times|\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}}}|}$
such that
$\boldsymbol{A}^{*}\boldsymbol{B}\boldsymbol{F}\boldsymbol{D}\approx\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}$
implies $\boldsymbol{y}\approx\boldsymbol{h}$ and hence
$\boldsymbol{f}\approx\boldsymbol{B}\boldsymbol{F}\boldsymbol{D}\boldsymbol{h}$.
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

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$.

  1. _(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.,

$x_{t,j}\coloneqq\left(-\frac{1}{2}+\frac{2t-1}{N_{1}},-\frac{1}{2}+\frac{2j-1}{N_{2}}\right)^{T}+\left(\frac{1}{N_{1}}\,\eta_{1},\frac{1}{N_{2}}\,\eta_{2}\right)^{T},$
(5.1)

$t=1,\dots,N_{1}$, $j=1,\dots,N_{2}$, and $\eta_{1},\eta_{2}\sim 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 1(a). Additionally,
we also consider the random grid
$x_{t,j}\coloneqq\frac{1}{2}\left(\eta_{1},\eta_{2}\right)^{T},$ see Figure
1(b).

(a) Jittered grid

(b) Random grid

Figure 5.1:  Exemplary randomized grids of size $N_{1}=N_{2}=12$.

  2. _(ii)_

Moreover, we examine grids of polar kind, as mentioned in [23]. For $R,T\in
2\mathbb{N}$ the points of the polar grid are given by a signed radius
$r_{j}\coloneqq\frac{j}{R}\in\left[-\frac{1}{2},\frac{1}{2}\right)$ and an
angle $\theta_{t}\coloneqq\frac{\pi
t}{T}\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right)$ as

$x_{t,j}\coloneqq
r_{j}\left(\cos\theta_{t},\sin\theta_{t}\right)^{T},\quad(j,t)^{T}\in\mathcal{I}_{R}\times\mathcal{I}_{T}.$
(5.2)

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}\coloneqq
r_{j}\left(\cos\theta_{t},\sin\theta_{t}\right)^{T},\quad(j,t)^{T}\in\mathcal{I}_{\sqrt{2}R}\times\mathcal{I}_{T},$
(5.3)

i. e., we added more concentric circles and excluded the points outside the
unit square, see Figure 2(a). 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 2(b). There we distinguish two sets of points, i. e.,

$x_{t,j}^{BH}\coloneqq\left(\frac{j}{R},\frac{4t}{T}\frac{j}{R}\right)^{T},\quad
x_{t,j}^{BV}\coloneqq\left(-\frac{4t}{T}\frac{j}{R},\frac{j}{R}\right)^{T},\quad(j,t)^{T}\in\mathcal{I}_{R}\times\mathcal{I}_{\frac{T}{2}}.$
(5.4)

(a) Polar (blue) and modified polar (red) grid

(b) Linogram / pseudo-polar grid

(c) Golden angle linogram grid

Figure 5.2:  Polar grids of size $R=12$ and $T=2R$.

  3. _(iii)_

Another modification of these polar grids was introduced in [33], where the
angles $\theta_{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

$\theta_{t}=\mod\left(\frac{\pi}{2}+t\,\frac{2\pi}{1+\sqrt{5}},\pi\right)-\frac{\pi}{2},\quad
t=0,\dots,T-1.$ (5.5)

The golden angle linogram grid is given by

$x_{t,j}\coloneqq\left\\{\begin{matrix}\left(\frac{2j+1}{2R},\,\frac{2j+1}{2R}\tan\left(\theta_{t}-\frac{\pi}{4}\right)\right)^{T}&\colon\theta_{t}\in\left[0,\frac{\pi}{2}\right)\\\
\left(-\frac{2j+1}{2R}\cot\left(\theta_{t}-\frac{\pi}{4}\right),\,\frac{2j+1}{2R}\right)^{T}&\colon\theta_{t}\in\left[-\frac{\pi}{2},0\right)\\\
\end{matrix}\right\\},\quad j\in\mathcal{I}_{R},$

with $\theta_{t}$ in (5.5), as illustrated in Figure 2(c).

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 coefficients $\hat{f}_{\boldsymbol{k}}\in[1,10]$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$.  In this test we consider
several $d\in\\{1,2,3\\}$. For the corresponding function evaluations of (2.8)
at given points
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{d}$,
$j=1,\dots,N$, we test how well these Fourier coefficients can be
approximated. More precisely, we consider the estimate
$\boldsymbol{\tilde{h}}^{\mathrm{w}}=\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}^{*}\boldsymbol{W}\boldsymbol{f}$,
cf. (3.7), with the matrix
$\boldsymbol{W}=\operatorname*{diag}\left(w_{j}\right)_{j=1}^{N}$ of density
compensation factors computed by means of  Algorithm 3.6, i. e., by (3.17), in
case $|\mathcal{I}_{\boldsymbol{2M}}|\leq N$, or by (3.18), if
$|\mathcal{I}_{\boldsymbol{2M}}|>N$, and compute the relative errors

$e_{p}\coloneqq\frac{\|\boldsymbol{\tilde{h}}^{\mathrm{w}}-\boldsymbol{\hat{f}}\|_{p}}{\|\boldsymbol{\hat{f}}\|_{p}},\quad
p\in\\{2,\infty\\}.$ (5.6)

By (3.39) it is known that

$\displaystyle\frac{\big{\|}\boldsymbol{\hat{f}}-\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{f}\big{\|}_{p}}{\big{\|}\boldsymbol{\hat{f}}\big{\|}_{p}}$
$\displaystyle\leq|\mathcal{I}_{\boldsymbol{M}}|\,\varepsilon,\quad
p\in\\{2,\infty\\},$

with the residual
$\varepsilon=\big{\|}\boldsymbol{A}_{|\mathcal{I}_{\boldsymbol{2M}}|}^{T}\,\boldsymbol{w}-\boldsymbol{e}_{\boldsymbol{0}}\big{\|}_{\infty}\geq
0$, cf. (3.38).

In our experiment we use random points $\boldsymbol{x}_{j}$ with
$N_{t}=2^{9-d}$, $t=1,\dots,d$, cf. Figure 1(b), and, for several problem
sizes  $\boldsymbol{M}=M\cdot\boldsymbol{1}_{d}$, $M=2^{c}$ with
$c=1,\dots,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}11-d},$
we choose random Fourier coefficients $\hat{f}_{\boldsymbol{k}}\in[1,10]$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$. Afterwards, we compute the
evaluations of the trigonometric polynomial (2.8) by means of an NFFT and use
the resulting vector $\boldsymbol{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 $|\mathcal{I}_{\boldsymbol{2M}}|<N$, i. e., as long as
$M<\frac{N_{1}}{2}=2^{8-d}$, the weights computed by means of (3.17) lead to
an exact reconstruction of the given Fourier coefficients. However, as soon as
we are in the setting $|\mathcal{I}_{\boldsymbol{2M}}|>N$ the least squares
approximation via (3.18) does not yield good results anymore.

(a) $d=1$

(b) $d=2$

(c) $d=3$

Figure 5.3:  Relative errors (5.6) of the reconstruction of the Fourier
coefficients of a trigonometric polynomial (2.8) with given
$\hat{f}_{\boldsymbol{k}}\in[1,10]$,
$\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$, computed via the density
compensation method from Algorithm 3.7, for random grids with
$N_{t}=2^{9-d}$, $t=1,\dots,d$, and
$\boldsymbol{M}=M\cdot\boldsymbol{1}_{d}$, $M=2^{c}$ with
$c=1,\dots,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}11-d}$.

######  Example 5.3.

In order to study the quality of our optimization method in Section 4, we
consider the Frobenius norms

$n(w,m,\sigma)\coloneqq\left\|\boldsymbol{A}^{*}\boldsymbol{B}\boldsymbol{F}\boldsymbol{D}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\right\|_{\mathrm{F}},\quad
n^{\textrm{opt}}(w,m,\sigma)\coloneqq\left\|\boldsymbol{A}^{*}\boldsymbol{B}_{\textrm{opt}}\boldsymbol{F}\boldsymbol{D}-\boldsymbol{I}_{|\mathcal{I}_{\boldsymbol{M}}|}\right\|_{\mathrm{F}},$
(5.7)

where $\boldsymbol{B}$ denotes the original matrix from the NFFT in (2.15) and
$\boldsymbol{B}_{\mathrm{opt}}$ the optimized matrix generated by Algorithm
4.3. For the original matrix $\boldsymbol{B}$ we utilize the common B-Spline
window function

$w_{\mathrm{B}}\coloneqq
B_{2m}(\boldsymbol{M}_{\boldsymbol{\sigma}}\boldsymbol{x})$ (5.8)

with the centered B-Spline of order $2m$, cf. [53, p. 388]. The optimized
matrix $\boldsymbol{B}_{\mathrm{opt}}$ 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
$\boldsymbol{M}=(12,12)^{T}$ as well as $N_{1}=N_{2}=R=2^{\mu}$,
$\mu\in\\{2,\dots,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}_{\boldsymbol{M}}|>N$ as well as for the overdetermined setting
$|\mathcal{I}_{\boldsymbol{M}}|\leq N$.

Having a look at the results for the grids in Remark 5.1, it becomes apparent
that they separate into two groups. Figure 4(a) 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 4(b)
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 $|\mathcal{I}_{\boldsymbol{M}}|\leq N$.

(a) Polar grid (5.2)

(b) Modified polar grid (5.3)

Figure 5.4:  Frobenius norms (5.7) of the original matrix $\boldsymbol{B}$
(violet) and the optimized matrix $\boldsymbol{B}_{\textrm{opt}}$ generated by
Algorithm 4.3 using the B-Spline $w_{\mathrm{B}}$ (orange) as well as the
Dirichlet window $w_{\mathrm{D}}$ (cyan) with $R=2^{\mu}$,
$\mu\in\\{2,\dots,7\\}$, and $T=2R$ as well as $\boldsymbol{M}=(12,12)^{T}$,
$m\in\\{2,4\\}$ and $\sigma\in\\{1,2\\}$.

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
$\boldsymbol{H}_{\boldsymbol{\ell}}^{*}\boldsymbol{H}_{\boldsymbol{\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 2(a), there are no grid points in the corners of the
unit square. Therefore, some of the index sets
$\mathcal{I}_{\boldsymbol{M}_{\boldsymbol{\sigma}},m}(\boldsymbol{\ell})$, 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 $\boldsymbol{B}$ from Section 4. To this
end, we show an example concerning trigonometric polynomials (2.8) of degree
$\boldsymbol{M}$ and a second one that deals with bandlimited functions of
bandwidth $\boldsymbol{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(a). Here we treat the phantom data as given Fourier
coefficients
$\boldsymbol{\hat{f}}\coloneqq(\hat{f}_{\boldsymbol{k}})_{\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}}$
of a trigonometric polynomial (2.8). For given points
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{2}$,
$j=1,\dots,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
$\boldsymbol{M}=(M,M)^{T}$, $M=2^{c}$ with
$c=3,\dots,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}10},$
and linogram grids (5.4) of size $R=2M$, $T=2R$, i. e., we consider the
setting $|\mathcal{I}_{\boldsymbol{2M}}|<N$. For using Algorithm 4.4 we choose
the oversampling factor $\sigma=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 $\boldsymbol{W}$ or the
computation of the optimized sparse matrix $\boldsymbol{B}_{\textrm{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 reconstruction
$\boldsymbol{\tilde{h}}\in\\{\boldsymbol{\tilde{h}}^{\mathrm{w}},\boldsymbol{h}_{\mathrm{opt}}\\}$,
cf. (3.7) and (4.12), we consider the relative errors

$e_{2}\coloneqq\frac{\|\boldsymbol{\tilde{h}}-\boldsymbol{\hat{f}}\|_{2}}{\|\boldsymbol{\hat{f}}\|_{2}}.$
(5.9)

The corresponding results can be found in Table 5.1. We remark that since we
are in the setting $|\mathcal{I}_{\boldsymbol{2M}}|<N$, the density
compensation method in Algorithm 3.7 with weights computed by (3.17) 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 $\boldsymbol{x}_{j}$ for
different input values $\boldsymbol{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.

|  Relative error $e_{2}$ | Precomputation time | Reconstruction time
---|---|---|---
$M$ |  Alg. 3.7 |  Alg. 4.4 |  Alg. 3.7 |  Alg. 4.4 |  Alg. 3.7 |  Alg. 4.4
8 | 1.3332e-15 | 6.8606e-14 | 9.8254e-02 | 1.9220e+00 | 6.2500e-04 | 2.2360e-03
16 | 7.2315e-15 | 1.5718e-07 | 1.6157e-01 | 8.3276e+00 | 2.5100e-03 | 3.4760e-03
32 | 2.3383e-14 | 4.5778e-07 | 3.3032e-01 | 4.3169e+01 | 3.1860e-03 | 7.4790e-03
64 | 2.5859e-14 | 4.7505e-07 | 3.4324e+00 | 2.4103e+02 | 5.0420e-03 | 4.9310e-03
128 | 7.9006e-14 | 5.9962e-07 | 9.4725e+00 | 1.2045e+03 | 2.9860e-02 | 5.5123e-02
256 | 2.6386e-13 | 4.0943e-06 | 3.8365e+01 | 5.8347e+03 | 6.6443e-02 | 6.7810e-01
512 | 1.0917e-12 | 2.0184e-06 | 1.4020e+02 | 2.9235e+04 | 2.1674e-01 | 3.2674e+00
1024 | 4.2563e-12 | 1.3491e-05 | 7.2153e+02 | 1.4342e+05 | 7.4896e-01 | 1.6114e+01
Table 5.1:  Relative errors (5.9) of the reconstruction of the Shepp-Logan
phantom of size $M$ as well as the runtime in seconds for the density
compensation method from Algorithm 3.7 compared to Algorithm 4.4 with
$\sigma=1.0$ and $m=4$, using linogram grids (5.4) of size $R=2M$, $T=2R$.

As a second experiment we aim to decrease the amount of overdetermination, i.
e., we want to keep the size $|\mathcal{I}_{\boldsymbol{M}}|$ of the phantom,
but reduce the number $N$ of the points $\boldsymbol{x}_{j}$, $j=1,\dots,N$.
To this end, we now consider linogram grids (5.4) of the smaller size $R=M$,
$T=2R$, i. e., now we have $|\mathcal{I}_{\boldsymbol{2M}}|>N$. The
reconstruction of the phantom of size $1024\times 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.18), 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
$|\mathcal{I}_{\boldsymbol{2M}}|<N$ (cf. Section 3.1.1) is violated, it can be
seen in Figure 5(c) that the density compensation method using weights
computed by means of (3.18) does not yield an exact reconstruction is this
setting. On the contrary, we recognize that our optimization method, see
Figure 5(d), achieves a huge improvement in comparison to the density
compensation techniques in Figure 5(b) and 5(c) since no artifacts are
visible. Presumably, this arises because there are more degrees of freedom in
the optimization of the matrix $\boldsymbol{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(a)
and 5(d) 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.

(a) Original phantom

(b) Voronoi weights
with $e_{2}$=5.3040e-01

(c) Algorithm 3.7
with $e_{2}$=5.0585e-01

(d) Algorithm 4.4
with $e_{2}$=2.2737e-03

Figure 5.5:  Reconstruction of the Shepp-Logan phantom of size $1024\times
1024$ (top) via the density compensation method from Section 3.1 using Voronoi
weights and Algorithm 3.7 with weights computed by (3.18) compared to
Algorithm 4.4 for the linogram grid (5.4) of size $R=M=1024$, $T=2R$; as well
as a detailed view of the 832nd row each (bottom).

######  Example 5.5.

Finally, we examine the reconstruction properties for bandlimited functions
$f\in L_{1}(\mathbb{R}^{d})\cap C_{0}(\mathbb{R}^{d})$ with maximum bandwidth
$\boldsymbol{M}$. To this end, we firstly specify a compactly supported
function $\hat{f}$ and consequently compute its inverse Fourier transform $f$,
such that its samples $f(\boldsymbol{x}_{j})$ for given
$\boldsymbol{x}_{j}\in\left[-\frac{1}{2},\frac{1}{2}\right)^{2}$,
$j=1,\dots,N$, can be used for the reconstruction of the samples
$\hat{f}(\boldsymbol{k})$, $\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$.
Here we consider the tensorized function
$\hat{f}(\boldsymbol{v})=g(v_{1})\cdot g(v_{2})$, where $g(v)$ is the one-
dimensional triangular pulse
$g(v)\coloneqq(1-\left|\tfrac{v}{b}\right|)\cdot\chi_{[-b,b]}(v)$. Then for
all $b\in\mathbb{N}$ with $b\leq\frac{M}{2}$ the associated inverse Fourier
transform

$\displaystyle
f(\boldsymbol{x})=\int_{\mathbb{R}^{2}}\hat{f}(\boldsymbol{v})\,\mathrm{e}^{2\pi\mathrm{i}\boldsymbol{v}\boldsymbol{x}}\,\mathrm{d}\boldsymbol{v}=b^{2}\,\mathrm{sinc}^{2}(b\pi\boldsymbol{x})=b^{2}\,\mathrm{sinc}^{2}(b\pi
x_{1})\,\mathrm{sinc}^{2}(b\pi x_{2}),\quad\boldsymbol{x}\in\mathbb{R}^{2},$

is bandlimited with bandwidth $\boldsymbol{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
$|\mathcal{I}_{\boldsymbol{2M}}|\leq 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.60), the weights in (3.64) and Algorithm 3.7 with weights computed via
(3.17), as well as Algorithm 4.4. For the reconstruction
$\boldsymbol{\tilde{h}}\in\\{\boldsymbol{\tilde{h}}^{\mathrm{w}},\boldsymbol{h}_{\mathrm{opt}}\\}$,
cf. (3.7) and (4.12), we then compute the pointwise absolute errors
$\big{|}\boldsymbol{\tilde{h}}-\boldsymbol{\hat{f}}\big{|}$. The corresponding
results are displayed in Figure 5.6. It can easily be seen that Voronoi
weights, see Figure 6(a), and the weights in (3.64), see Figure 6(d), 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.22) 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.60), see Figure 6(c), as well as the weights
in (3.64), see Figure 6(d), is affordable, cf. Section 3.4.2. In contrast, our
new methods using Algorithm 3.7, see Figure 6(b), or Algorithm 4.4, see Figure
6(e), are much more effective and therefore better suited for the given
problem.

(a) Voronoi weights

(b) Algorithm 3.7

(c) Weights via (3.60)

(d) Weights (3.64)

(e) Algorithm 4.4

Figure 5.6:  Pointwise absolute error
$\big{|}\boldsymbol{\tilde{h}}-\boldsymbol{\hat{f}}\big{|}$ of the
reconstruction of samples $\hat{f}(\boldsymbol{k})$ of the tensorized
triangular pulse function with $M=64$ and $b=24$, using the density
compensation methods considered in Section 3 as well as the optimization
approach from Algorithm 4.4 for the jittered grid (5.1) of size
$N_{1}=N_{2}=144$.

##  6 Conclusion

In the present paper we considered several direct methods for computing an
inverse NFFT, i. e., reconstructing the Fourier coefficients
$\hat{f}_{\boldsymbol{k}}$, $\boldsymbol{k}\in\mathcal{I}_{\boldsymbol{M}}$,
from given nonequispaced data $f(\boldsymbol{x}_{j})$, $j=1,\dots,N$. Being a
direct method here means, that for a fixed set of points $\boldsymbol{x}_{j}$,
$j=1,\dots,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 $\boldsymbol{x}_{j}$, $j=1,\dots,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\geq 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}\in\mathbb{C}$, $j=1,\dots,N,$ are
precomputed, such that the reconstruction can be realized by means of an
adjoint NFFT applied to the scaled data $w_{j}f(\boldsymbol{x}_{j})$. We
started our investigations with trigonometric polynomials in Section 3.1. In
Corollary 3.4 we introduced the main formula (3.16), that yields exact
reconstruction for all trigonometric polynomials of degree $\boldsymbol{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
$\boldsymbol{A}^{*}\boldsymbol{W}\boldsymbol{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
$\boldsymbol{A}\approx\boldsymbol{B}\boldsymbol{F}\boldsymbol{D}$ of the NFFT
is used to modify the sparse matrix $\boldsymbol{B}$, such that a
reconstruction is given by
$\boldsymbol{\hat{f}}\approx\boldsymbol{D}^{*}\boldsymbol{F}^{*}\boldsymbol{B}_{\mathrm{opt}}^{*}\boldsymbol{f}$.
In other words, the inversion is done by a modified adjoint NFFT, while the
optimization of the matrix $\boldsymbol{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
$|\mathcal{I}_{\boldsymbol{M}}|\leq N$ and work for many different sampling
patterns. More specifically, in the highly overdetermined case
$|\mathcal{I}_{\boldsymbol{2M}}|\leq 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 $|\mathcal{I}_{\boldsymbol{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.

## Acknowledgments

Melanie Kircheis gratefully acknowledges the support from the BMBF grant
01$\mid$S20053A (project SA$\ell$E). Daniel Potts acknowledges the funding by
Deutsche Forschungsgemeinschaft (German Research Foundation) – Project–ID
416228727 – SFB 1410.

Moreover, the authors thank the referees and the editor for their very useful
suggestions for improvements.

## References

  * [1] R. Archibald, A. Gelb, and R. B. Platte.  Image reconstruction from undersampled Fourier data using the polynomial annihilation transform.  J. Sci. Comput., 67(2):432–452, 2016.
  * [2] R. Archibald, A. Gelb, and J. Yoon.  Polynomial fitting for edge detection in irregularly sampled signals and images.  SINUM, 43(1):259–279, 2005.
  * [3] A. Averbuch, R. Coifman, D. Donoho, M. Israeli, and Y. Shkolnisky.  A framework for discrete integral transformations I – the pseudopolar Fourier transform.  SIAM J. Sci. Comput., 30:764–784, 2008.
  * [4] A. Averbuch, R. Coifman, D. L. Donoho, M. Elad, and M. Israeli.  Fast and accurate polar Fourier transform.  Appl. Comput. Harmon. Anal., 21:145–167, 2006.
  * [5] A. Averbuch, G. Shabat, and Y. Shkolnisky.  Direct inversion of the three-dimensional pseudo-polar Fourier transform.  SIAM J. Sci. Comput., 38(2):A1100–A1120, 2016.
  * [6] R. F. Bass and K. Gröchenig.  Random sampling of multivariate trigonometric polynomials.  SIAM J. Math. Anal., 36:773–795, 2004.
  * [7] G. Beylkin.  On the fast Fourier transform of functions with singularities.  Appl. Comput. Harmon. Anal., 2:363–381, 1995.
  * [8] Å. Björck.  Numerical Methods for Least Squares Problems.  SIAM, Philadelphia, PA, USA, 1996.
  * [9] A. Böttcher and D. Potts.  Probability against condition number and sampling of multivariate trigonometric random polynomials.  Electron. Trans. Numer. Anal., 26:178–189, 2007.
  * [10] E. J. Candès, M. B. Wakin, and S. P. Boyd.  Enhancing sparsity by reweighted $\ell_{1}$ minimization.  J. Fourier Anal. Appl., 14(5):877–905, 2008.
  * [11] R. Chartrand and W. Yin.  Iteratively reweighted algorithms for compressive sensing.  In IEEE International Conference on Acoustics, Speech and Signal Processing, pages 3869–3872. IEEE, 2008.
  * [12] H. Choi and D. C. Munson.  Analysis and design of minimax-optimal interpolators.  IEEE Trans. Signal Process., 46(6):1571–1579, 1998.
  * [13] V. Churchill, R. Archibald, and A. Gelb.  Edge-adaptive $\ell_{2}$ regularization image reconstruction from non-uniform Fourier data.  Inverse Probl. Imaging, 13(5):931–958, 2019.
  * [14] I. Daubechies, R. DeVore, M. Fornasier, and C. S. Güntürk.  Iteratively reweighted least squares minimization for sparse recovery.  Commun. Pure Appl. Math, 63(1):1–38, 2010.
  * [15] M. Doneva, M. Akcakaya, and C. Prieto, editors.  Magnetic Resonance Image Reconstruction: Theory, Methods and Applications, volume 6.  Academic Press, 2022.
  * [16] A. J. W. Duijndam and M. A. Schonewille.  Nonuniform fast Fourier transform.  Geophysics, 64:539–551, 1999.
  * [17] A. Dutt and V. Rokhlin.  Fast Fourier transforms for nonequispaced data.  SIAM J. Sci. Stat. Comput., 14:1368–1393, 1993.
  * [18] A. Dutt and V. Rokhlin.  Fast Fourier transforms for nonequispaced data II.  Appl. Comput. Harmon. Anal., 2:85–100, 1995.
  * [19] H. Eggers, M. Kircheis, and D. Potts.  Non-Cartesian MRI reconstruction.  In M. Doneva, M. Akcakaya, and C. Prieto, editors, Magnetic Resonance Image Reconstruction: Theory, Methods and Applications, volume 6. Academic Press, 2022.
  * [20] B. Elbel and G. Steidl.  Fast Fourier transform for nonequispaced data.  In C. K. Chui and L. L. Schumaker, editors, Approximation Theory IX, pages 39–46, Nashville, 1998. Vanderbilt University Press.
  * [21] G. E. Fasshauer.  Meshfree approximation methods with MATLAB.  World Scientific Publishers, 2007.
  * [22] H. G. Feichtinger, K. Gröchenig, and T. Strohmer.  Efficient numerical methods in non-uniform sampling theory.  Numer. Math., 69:423–440, 1995.
  * [23] M. Fenn, S. Kunis, and D. Potts.  On the computation of the polar FFT.  Appl. Comput. Harmon. Anal., 22:257–263, 2007.
  * [24] J. A. Fessler and B. P. Sutton.  Nonuniform fast Fourier transforms using min-max interpolation.  IEEE Trans. Signal Process., 51:560–574, 2003.
  * [25] K. Fourmont.  Non equispaced fast Fourier transforms with applications to tomography.  J. Fourier Anal. Appl., 9:431–450, 2003.
  * [26] A. Gelb and G. Song.  A frame theoretic approach to the nonuniform fast Fourier transform.  SIAM J. Numer. Anal., 52(3):1222–1242, 2014.
  * [27] A. Gopal and V. Rokhlin.  A fast procedure for the construction of quadrature formulas for bandlimited functions.  Technical Report, YALEU/DCS/TR-1563, 2022.
  * [28] M. Gräf, S. Kunis, and D. Potts.  On the computation of nonnegative quadrature weights on the sphere.  Appl. Comput. Harmon. Anal., 27:124–132, 2009.
  * [29] L. Greengard and J.-Y. Lee.  Accelerating the nonuniform fast Fourier transform.  SIAM Rev., 46:443–454, 2004.
  * [30] L. Greengard, J.-Y. Lee, and S. Inati.  The fast sinc transform and image reconstruction from nonuniform samples in $k$-space.  Commun. Appl. Math. Comput. Sci., 1:121–131, 2006.
  * [31] K. Gröchenig.  Sampling, Marcinkiewicz–Zygmund inequalities, approximation, and quadrature rules.  J. Approx. Theory, 257:105455, 2020.
  * [32] G. Heinig and K. Rost.  Algebraic methods for Toeplitz-like matrices and operators, volume 19 of Mathematical Research.  Akademie-Verlag, Berlin, 1984.
  * [33] E. S. Helou et al.  The discrete Fourier transform for golden angle linogram sampling.  Inverse Problems, 35(125004), 2019.
  * [34] D. Huybrechs.  Stable high-order quadrature rules with equidistant points.  J. Comput. Appl. Math., 231(2):933–947, 2009.
  * [35] J. Keiner, S. Kunis, and D. Potts.  NFFT 3.5, C subroutine library.  http://www.tu-chemnitz.de/~potts/nfft.  Contributors: F. Bartel, M. Fenn, T. Görner, M. Kircheis, T. Knopp, M. Quellmalz, M. Schmischke, T. Volkmer, A. Vollrath.
  * [36] J. Keiner, S. Kunis, and D. Potts.  Using NFFT3 - a software library for various nonequispaced fast Fourier transforms.  ACM Trans. Math. Software, 36:Article 19, 1–30, 2009.
  * [37] M. Kircheis and D. Potts.  Direct inversion of the nonequispaced fast Fourier transform.  Linear Algebra Appl., 575:106–140, 2019.
  * [38] M. Kircheis and D. Potts.  Efficient multivariate inversion of the nonequispaced fast Fourier transform.  PAMM, 20(1):e202000120, 2021.
  * [39] M. Kircheis, D. Potts, and M. Tasche.  Nonuniform fast Fourier transforms with nonequispaced spatial and frequency data and fast sinc transforms.  Numer. Algor., 2022.
  * [40] M. Kircheis, D. Potts, and M. Tasche.  On regularized Shannon sampling formulas with localized sampling.  Sampl. Theory Signal Process. Data Anal., 20(20), 2022.
  * [41] T. Knopp, S. Kunis, and D. Potts.  A note on the iterative MRI reconstruction from nonuniform k-space data.  Int. J. Biomed. Imag., 2007, 2007.  ID 24727.
  * [42] V. A. Kotelnikov.  On the transmission capacity of the “ether” and wire in electrocommunications.  In Modern Sampling Theory: Mathematics and Application, pages 27–45. Birkhäuser, Boston, 2001.  Translated from Russian.
  * [43] S. Kunis and D. Nagel.  On the condition number of Vandermonde matrices with pairs of nearly-colliding nodes.  Numer. Algor., 87:473–496, 2021.
  * [44] S. Kunis and D. Potts.  Stability results for scattered data interpolation by trigonometric polynomials.  SIAM J. Sci. Comput., 29:1403–1419, 2007.
  * [45] S. Kunis and D. Potts.  Time and memory requirements of the nonequispaced FFT.  Sampl. Theory Signal Image Process., 7:77–100, 2008.
  * [46] J.-Y. Lee and L. Greengard.  The type 3 nonuniform FFT and its applications.  J. Comput. Physics, 206:1–5, 2005.
  * [47] R. Lin and H. Zhang.  Convergence analysis of the Gaussian regularized Shannon sampling formula.  Numer. Funct. Anal. Optim., 38(2):224–247, 2017.
  * [48] Y. Liu, J. Ma, Y. Fan, and Z. Liang.  Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction.  Phys. Med. Biol., 57(23):7923, 2012.
  * [49] J. Lund and K. L. Bowers.  Sinc Methods for Quadrature and Differential Equations.  Society for Industrial and Applied Mathematics, 1992.
  * [50] C. Micchelli, Y. Xu, and H. Zhang.  Optimal learning of bandlimited functions from localized sampling.  J. Complexity, 25(2):85–114, 2009.
  * [51] A. Nieslony and G. Steidl.  Approximate factorizations of Fourier matrices with nonequispaced knots.  Linear Algebra Appl., 266:337–351, 2003.
  * [52] J. G. Pipe and P. Menon.  Sampling density compensation in MRI: rationale and an iterative numerical solution.  Magn. Reson. Med., 41:179–186, 1999.
  * [53] G. Plonka, D. Potts, G. Steidl, and M. Tasche.  Numerical Fourier Analysis.  Applied and Numerical Harmonic Analysis. Birkhäuser, 2018.
  * [54] D. Potts and M. Tasche.  Continuous window functions for NFFT.  Adv. Comput. Math., 47(53):1–34, 2021.
  * [55] D. Potts and M. Tasche.  Uniform error estimates for nonequispaced fast Fourier transforms.  Sampl. Theory Signal Process. Data Anal., 19(17):1–42, 2021.
  * [56] K. P. Pruessmann and F. T. A. W. Wayer.  Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories.  In Proc. Intl. Soc. Mag. Reson. Med. 9, page 767. Glasgow, Scotland, 2001.
  * [57] L. Qian.  On the regularized Whittaker–Kotelnikov–Shannon sampling formula.  Proc. Amer. Math. Soc., 131(4):1169–1176, 2003.
  * [58] V. Rasche, R. Proksa, R. Sinkus, P. Börnert, and H. Eggers.  Resampling of data between arbitrary grids using convolution interpolation.  IEEE Trans. Med. Imag., 18:385–392, 1999.
  * [59] D. Rosenfeld.  An optimal and efficient new gridding algorithm using singular value decomposition.  Magn. Reson. Med., 40(1):14–23, 1998.
  * [60] D. Ruiz-Antolin and A. Townsend.  A nonuniform fast Fourier transform based on low rank approximation.  SIAM J. Sci. Comput., 40(1):A529–A547, 2018.
  * [61] T. Scarnati and A. Gelb.  Accelerated variance based joint sparsity recovery of images from Fourier data.  arXiv, 2019.
  * [62] H. Sedarat and D. G. Nishimura.  On the optimality of the gridding reconstruction algorithm.  IEEE Trans. Med. Imaging, 19(4):306–317, 2000.
  * [63] J. Selva.  Efficient type-4 and type-5 non-uniform FFT methods in the one-dimensional case.  IET Signal Processing, 12(1):74–81, 2018.
  * [64] C. E. Shannon.  Communication in the presence of noise.  Proc. I.R.E., 37:10–21, 1949.
  * [65] W. Stefan, R. A. Renaut, and A. Gelb.  Improved total variation-type regularization using higher order edge detectors.  SIIMS, 3(2):232–251, 2010.
  * [66] G. Steidl.  A note on fast Fourier transforms for nonequispaced grids.  Adv. Comput. Math., 9:337–353, 1998.
  * [67] G. Stewart.  Matrix Algorithms: Basic decompositions.  SIAM, Philadelphia, 1998.
  * [68] T. Strohmer and J. Tanner.  Fast reconstruction methods for bandlimited functions from periodic nonuniform sampling.  SIAM J. Numer. Anal., 44(3):1071–1094, 2006.
  * [69] B. P. Sutton, J. A. Fessler, and D. C. Noll.  A min-max approach to the nonuniform $N$-dimensional FFT for rapid iterative reconstruction of MR images.  In Proc. ISMRM 9th Scientific Meeting, page 763, 2001.
  * [70] V. Tchakaloff.  Formules de cubature mécaniques à coefficients non négatifs.  Bull. Sci. Math., 81:123–134, 1957.
  * [71] G. Wasserman, R. Archibald, and A. Gelb.  Image reconstruction from Fourier data using sparsity of edges.  J. Sci. Comput., 65(2):533–552, 2015.
  * [72] E. T. Whittaker.  On the functions which are represented by the expansions of the interpolation theory.  Proc. R. Soc. Edinb., 35:181–194, 1915.


# Analyzing Political Figures in Real-Time: Leveraging YouTube Metadata for
Sentiment Analysis

1st Danendra Athallariq Harya Putra  Master Program of Informatics
School of Electrical Engineering and Informatics
Institut Teknologi Bandung
Bandung, Indonesia
<EMAIL_ADDRESS>2nd Arief Purnama Muharram  Master Program of
Informatics
School of Electrical Engineering and Informatics
Institut Teknologi Bandung
Bandung, Indonesia
<EMAIL_ADDRESS>
###### Abstract

Sentiment analysis using big data from YouTube videos metadata can be
conducted to analyze public opinions on various political figures who
represent political parties. This is possible because YouTube has become one
of the platforms for people to express themselves, including their opinions on
various political figures. The resulting sentiment analysis can be useful for
political executives to gain an understanding of public sentiment and develop
appropriate and effective political strategies. This study aimed to build a
sentiment analysis system leveraging YouTube videos metadata. The sentiment
analysis system was built using Apache Kafka, Apache PySpark, and Hadoop for
big data handling; TensorFlow for deep learning handling; and FastAPI for
deployment on the server. The YouTube videos metadata used in this study is
the video description. The sentiment analysis model was built using LSTM
algorithm and produces two types of sentiments: positive and negative
sentiments. The sentiment analysis results are then visualized in the form a
simple web-based dashboard.

###### Index Terms:

sentiment analysis, big data, politics, YouTube

##  I Introduction

General Elections (Pemilu) is one of the concrete manifestations of a
democratic system. Through Pemilu, the public has the opportunity to
participate in governance by electing their representatives who will represent
them in the government structure [1]. Among the various types of elections,
the Presidential Election (Pilpres) is always a highly anticipated moment and
dubbed as the largest ”democratic party”. In 2024, Indonesia will hold a
Pilpres to determine the candidate for the presidency who will lead Indonesia
for the next 5 years.

Welcoming the Pilpres 2024, every political party is competing to determine
the best presidential and vice-presidential candidate to be endorsed. For
political parties, Pilpres is not only about the positions of President and
Vice President, but also determines their seats in the future government
structure. Therefore, it is crucial for political parties to devise the best
political campaign strategies to win the hearts of the public. One of the
efforts that political parties can undertake to evaluate the quality of their
political figures is through public sentiment analysis.

Sentiment analysis, also known as opinion mining, is a field that studies the
analysis of opinions, sentiments, evaluations, judgments, attitudes, and
emotions of people towards entities such as products, services, organizations,
individuals, issues, events, and others related [2]. Public sentiment analysis
can be used as a tool for political parties to gain a better understanding of
the opinions and views of the public towards their endorsed political
candidates. With public sentiment analysis, political parties can design
effective and responsive campaign strategies to meet the needs of the public.

In the current digital era, social media has become a platform for the public
to express various things, including their views towards various political
figures. Such expressions can be in the form of support or rejection, and can
be expressed in various media such as text, audio, or video. Such expressions
can be used as indicators of public sentiment for political parties to assess
the quality of their endorsed political figures.

This research aims to build a ’real-time’ sentiment analysis system for
political figures in Pilpres 2024 from various videos on YouTube through its
video description as the metadata. The selection of YouTube as a source of big
data is due to the fact that YouTube has become one of the means of political
expression in the form of videos with various purposes [3, 4, 5]. The system
is designed to be ’real-time’ in order to perform sentiment analysis on
various YouTube videos metadata in ’real-time’. The resulting system is
intended for political executives to help gain an understanding of public
sentiment towards their endorsed political figures so that they can devise
appropriate and effective political strategies.

##  II Methodology

Figure 1: System design

###  II-A System Design

The sentiment analysis system was built using Apache Kafka, Apache PySpark,
and Hadoop for handling big data; TensorFlow for deep learning; and FastAPI
for deployment on the server. In terms of architecture, the system was built
using a module-based approach and consists of modules including producer,
streamer, HDFS, inference, and visualizer (Table I). The system’s workflow
involves data retrieval (crawling) through the YouTube API by the producer
module; internal data streaming by the streamer module and storing it into the
Hadoop Distributed File System (HDFS) by the HDFS module; sentiment inference
by the inference module; and displaying the sentiment inference results in a
simple web dashboard by the visualizer module (Figure 1). The producer module
can be set to perform data crawling on a scheduled and regular basis and then
store the results into HDFS to ensure that the data in the system is always
up-to-date in real-time.

TABLE I: Modules in the System Module |  Function |  Technology Used
---|---|---
Producer Module |  Responsible for scheduled data crawling through the YouTube API. |  Apache Kafka
Streamer Module |  Responsible for capturing data produced by the producer and storing it in HDFS. |  Apache PySpark
HDFS Module |  Responsible for organizing data within HDFS. |  Hadoop
Inference Module |  Responsible for conducting sentiment inference. |  TensorFlow
Visualizer Module |  Responsible for presenting the sentiment inference results in the form of a web dashboard. |  Python, FastApi (Backend); HTML, CSS, Javascript, Bootstrap and Chart.js (Frontend)

Figure 2: Data gathering

Figure 3: Dashboard design

Figure 4: Visualizer implementation design

Figure 5: Dashboard page testing results using production data

####  II-A1 Data Gathering via Youtube API

To retrieve data streams from the YouTube API, a Kafka producer with a topic
is needed that will send the metadata of a YouTube video. This task will then
be performed by the producer module. Metadata is obtained using the search
method provided by the YouTube search API, and the final result is in JSON
format.

The search method requires a keyword that will be used to perform the search.
This kind of keyword is similar to what we do when typing keywords for
searching videos on YouTube. The keywords used will be adjusted based on each
political figure, so that one political figure can use many keywords according
to the wishes of the user. In this study, only the name of the political
figure was used as a keyword (’anies’ for Anies Rasyid Baswedan, ’ganjar’ for
Ganjar Pranowo, ’prabowo’ for Prabowo Subianto, and ’puan’ for Puan Maharani).

When using the search method from the YouTube API, the results obtained are
not all videos related to the given keywords. Instead, these related videos
will be divided into several sections (pages). Each page contains a maximum of
50 related videos, so an iterative method is needed to continue the results
that have been obtained before. This can be done by passing the pageToken
parameter when searching. The pageToken parameter is obtained from the
metadata of the search method, specifically the nextPage section. Therefore, a
keyword will also be iterated as long as nextPage from the previous metadata
is not None.

The metadata properties taken for this project from the response data are
videoId, channelId, title, description, and uri. VideoId was used as a
differentiator between videos, so there is no duplicate data saved. When
saving to HDFS, the metadata will be combined with the political figures and
also the keywords data that were searched. Although there are several items of
metadata information saved, only the description will be further used for
sentiment analysis. Figure 2 illustrates this entire process.

####  II-A2 Storing Stream Data to HDFS

The output of the previous process (by the producer module) will be captured
using Kafka in the form of dstream data. First, this dstream data must be
converted into a Spark DataFrame format. Next, the Spark DataFrame is then
saved in the form of a parquet file in the HDFS folder. This task is performed
by the streamer module.

####  II-A3 Sentiment Analysis Model

To perform sentiment inference on each YouTube video metadata, a sentiment
analysis model is required. The metadata used is only the video description.
The model used in this research is based on the Long Short-Term Memory (LSTM)
algorithm [6] built using the TensorFlow library. LSTM is used because it is
one of the popular methods for solving text-related problems. This is because
LSTM is a type of Recurrent Neural Network (RNN) that is suitable for text-
related problems due to its consideration of the previous input context for
the next input. This property is in line with the natural property of text,
which needs to consider the context of previous words to understand the
meaning of the current word in a sentence.

The model is then trained using the sentiment dataset from the 2017 Jakarta
gubernatorial election (Pilkada) [7]. The dataset was taken from social media
Twitter related to the 2017 Jakarta gubernatorial election and consists of two
types of sentiment, positive and negative, with an equal number of each type.
The dataset used has undergone two preprocessing steps, namely replacing
emojis with special markers and removing stopwords. The dataset will be
divided into two: training data and validation data. The training data is used
to train the sentiment analysis model, while the validation data is used to
test the performance of the model on data that has not been seen before by the
model. We used a training to validation data ratio of 0.8:0.2 and proceeded to
train our model on Google Colab.

Before training or inference on the model, the data used needs to undergo
preprocessing first. The necessary preprocessing includes cleaning the text by
removing URLs, hashtags, mentions, and emojis. The cleaned data will then be
tokenized using text vectorization from TensorFlow. In addition to text
preprocessing, label conversion to one-hot encoding is also required. The
cleaned data will then be fed into the model for training.

The architecture of the sentiment analysis model used is as follows: First,
the data will enter an embedding layer to convert the tokenized text into an
embedding vector. Then, this embedding vector will be passed through an LSTM
layer and two dense layers. The last dense layer will be used for
classification.

####  II-A4 The Visualization of Sentiment Inference Results

The aggregation of sentiment inference results is displayed through a simple
web dashboard by the visualizer module. The visualizer module consists of two
submodules, namely backend and frontend. The backend module is responsible for
preparing and aggregating the required data by communicating directly with
HDFS, while the frontend module is responsible for visualizing the data that
has been prepared by the backend into the dashboard page. The backend module
is developed using FastAPI, while the frontend module is developed using
Bootstrap and Chart.js. On the dashboard page, the sentiment results will be
displayed in the form of a doughnut chart comparing the number of positive and
negative sentiments to facilitate readability (Figure 3). In the production
stage implementation, the frontend can be placed on the app server, while the
backend can be placed on the big data server (Figure 4).

###  II-B Evaluation Strategy

####  II-B1 Sentiment Model Evaluation

The evaluation of the sentiment model is performed using the validation
dataset. The evaluation metrics used are precision (1), recall (2), and
F1-score (3) for each type of sentiment, and accuracy (4) for overall
performance.

$precision=\frac{TP}{FP+TP}$ (1)

$recall=\frac{TP}{FN+TP}$ (2)

$F1score=\frac{2\times precision\times recall}{precision+recall}$ (3)

$accuracy=\frac{TP+TN}{TP+FP+FN+TN}$ (4)

####  II-B2 System Evaluation

The evaluation of the developed system is conducted through a usability
testing approach. In this evaluation, the system is locally deployed and
assessed for its functionality and potential weaknesses.

##  III Result

Using validation data from the 2017 Jakarta Gubernatorial Election dataset,
the F1-Score was obtained as shown in Table 3. The LSTM model algorithm
training with 8 epochs was able to provide overall accuracy performance of
0.7056. In terms of sentiment label evaluation, the resulting model has the
same precision value between the two labels. However, the recall value for the
negative label is better than the positive label. The difference in recall
values has an impact on the F1-score, where the negative label has a better
F1-score than the positive label.

TABLE II: Performance of Sentiment Analysis Model with Validation Dataset Sentiment | Precision | Recall | F1 Score
---|---|---|---
Negative | 0.72 | 0.82 | 0.76
Positive | 0.72 | 0.59 | 0.65
Overall | 0.72 | 0.7 | 0.71

We then tested the system by deploying it locally. The results of sentiment
inference are displayed in the form of aggregated numbers of negative and
positive sentiments for each political figure up to the date of the page
request. This information is then presented in the form of a doughnut chart.
Figure 5 illustrates the system. Positive sentiment is given a green color,
while negative sentiment is given a red color.

##  IV Discussion

The sentiment analysis model used in this study is LSTM, trained with the
Pilkada DKI 2017 sentiment dataset [5]. The testing results with that dataset
were able to produce a model performance of 0.7056. In the analysis of the
sentiment label, although both labels have the same precision value (0.72),
there is a significant difference between the recall values for negative
sentiment (0.82) and positive sentiment (0.59). Based on the recall formula
(2), the higher the recall value, the more the model can classify correctly.
The high recall value for negative sentiment implies a tendency for the model
to classify negative sentiment more than positive sentiment.

However, despite the acceptable performance of the model, further studies are
needed to improve the model’s performance, as there are several other factors
that might affect the system’s performance and become limitations in this
study.

  * •

The video search process is highly influenced by the selected keywords, so it
is necessary to choose the appropriate keywords to increase the expected
number of relevant video searches.

  * •

The video search process is limited by the YouTube API call limitations on the
free tier, which is 2,000 calls per day.

  * •

The model inference only used the video description, assuming there is a
correspondence between the content and the video description (not clickbait).

##  V Conclusion

YouTube, as the largest video platform, has been used as a political
expression medium for the public. This study has successfully developed a
political sentiment analysis system that leverages YouTube as a big data
source. Using the LSTM algorithm, the built-in inference model can provide
accuracy performance up to 0.7056. The keyword selection and the use of other
model algorithms (such as deep learning) can be considered for future research
to improve the resulting inference model performance.

## Acknowledgment

We would like to express our gratitude to the lecturer of IF5270 Big Data
Systems course at the Master of Informatics Study Program, School of
Electrical Engineering and Informatics, Institut Teknologi Bandung, who has
taught us the fundamental concepts and technologies of big data systems. This
research was initiated as a task for the course. We also would like to thank
our colleagues who have supported our learning process during the class.

## References

  * [1] Suryanto, Pengantar Ilmu Politik. Bandung: Pustaka Setia, 2018.
  * [2] B. Liu, Sentiment Analysis and Opinion Mining. Morgan & Claypool Publishers, 2012.
  * [3] E. W. Wirga, ”Analisis Konten pada Media Sosial Video YouTube untuk Mendukung Strategi Kampanye Politik,” in Jurnal Ilmiah Informatika dan Komputer, vol. 21, no. 1, pp. 14-26, 2016. [Online]. Available: https://ejournal.gunadarma.ac.id/index.php/infokom/article/view/1716. Accessed: September 20, 2023.
  * [4] K. Arofah, ”YouTube sebagai Media Klarifikasi dan Pernyataan Tokoh Politik,” in Jurnal Ilmu Komunikasi, vol. 13, pp. 111-123, 2015. [Online]. Available: http://jurnal.upnyk.ac.id/index.php/komunikasi/article/view/1442. Accessed: September 20, 2023.
  * [5] S. Thomas et al., ”Studi Analisis Metode Analisis Sentimen pada YouTube,” in Journal of Information Technology, vol. 1, no. 1, pp. 1-7, 2021, doi:10.46229/jifotech.v1i1.201.
  * [6] S. Hochreiter and J. Schmidhuber, ”Long Short-Term Memory,” in Neural Computation, vol. 9, no. 8, pp. 1735-1780, 1997. doi: 10.1162/neco.1997.9.8.1735
  * [7] A. R. T. Lestari, R. S. Perdana, and M. A. Fauzi, ”Analisis Sentimen Tentang Opini Pilkada DKI 2017 Pada Dokumen Twitter Berbahasa Indonesia Menggunakan Naive Bayes dan Pembobotan Emoji,” in Jurnal Pengembangan Teknologi Informasi Dan Ilmu Komputer, vol. 1, no. 12, pp. 1718-1724, 2017. [Online]. Available: https://j-ptiik.ub.ac.id/index.php/j-ptiik/article/view/627. Accessed: September 20, 2023