The algorithm

In summary, the following Algorithm 1 is obtained for the fast computation of (2.3) with ${\cal O} \left(\left\vert I_{\ensuremath{\boldsymbol{n}}}\right\vert \log \left\vert I_{\ensuremath{\boldsymbol{n}}}\right\vert + m M\right)$ arithmetic operations.

Algorithm 1 reads in matrix vector notation as

$$\ensuremath{\boldsymbol{A}} \ensuremath{\boldsymbol{\hat f}} \app... ...h{\boldsymbol{F}} \ensuremath{\boldsymbol{D}} \ensuremath{\boldsymbol{\hat f}},$$

where $\ensuremath{\boldsymbol{B}}$ denotes the real $M \times \left\vert I_{\ensuremath{\boldsymbol{n}}}\right\vert$ sparse matrix

$$\ensuremath{\boldsymbol{B}} := \left( \tilde\psi\left(\ensuremath... ...,\hdots,M-1;\, \ensuremath{\boldsymbol{l}}\in I_{\ensuremath{\boldsymbol{n}}}},$$

where $\ensuremath{\boldsymbol{F}}$ is the Fourier matrix of size $\left\vert I_{\ensuremath{\boldsymbol{n}}}\right\vert \times \left\vert I_{\ensuremath{\boldsymbol{n}}}\right\vert$ , and where $\ensuremath{\boldsymbol{D}}$ is the real $\left\vert I_{\ensuremath{\boldsymbol{n}}}\right\vert \times \vert I_{\ensuremath{\boldsymbol{N}}}\vert$ 'diagonal' matrix

$$\ensuremath{\boldsymbol{D}} := \bigotimes_{t=0}^{d-1} \left(\ensu... ...t)_{k_t\in I_{N_t}} \, \vert \, \ensuremath{\boldsymbol{O}}_t \right)^{\rm {T}}$$

with zero matrices $\ensuremath{\boldsymbol{O}}_t$ of size $N_t \times \frac{n_t-N_t}{2}$ .

The corresponding computation of the adjoint matrix vector product reads as

$$\ensuremath{\boldsymbol{A}}^{{\vdash \hspace*{-1.72mm} \dashv}} \... ...ashv}} \ensuremath{\boldsymbol{B}}^{\top} \ensuremath{\boldsymbol{\hat f}} \, .$$

With the help of the transposed index set

$$I_{\ensuremath{\boldsymbol{n}},m}^{\top}\left(\ensuremath{\boldsy... ...ol{x}}_j \le \ensuremath{\boldsymbol{l}}+m \ensuremath{\boldsymbol{1}}\right\},$$

one obtains Algorithm 2 for the adjoint NFFT.

Due to the characterisation of the nonzero elements of the matrix $\ensuremath{\boldsymbol{B}}$ , i.e.,

$$\bigcup\limits_{j=0}^{M-1} j \times I_{\ensuremath{\boldsymbol{n}... ...op}\left(\ensuremath{\boldsymbol{l}}\right) \times \ensuremath{\boldsymbol{l}}.$$

the multiplication with the sparse matrix $\ensuremath{\boldsymbol{B}}^{\top}$ is implemented in a 'transposed' way in the library, summation as outer loop and only using the multi-index sets $I_{\ensuremath{\boldsymbol{n}},m} \left({\ensuremath{\boldsymbol{x}}}_j\right)$ .