Jump to main content
Nonequispaced fast Fourier transform
Least squares

Least squares

Let $ \ensuremath{\boldsymbol{N}}=(N_1,\ldots,N_d)^{\top}\in 2\mathbb{N}^d$. For given samples $ (\ensuremath{\boldsymbol{x}}_j,y_j)\in \mathbb{T}^d\times\mathbb{C},\; j=0,\hdots,M-1$, the index set $ I_{\ensuremath{\boldsymbol{N}}}^d:=
\dots\times\big\{-\frac{N_d}{2},\hdots,\frac{N_d}{2}-1\big\}$, of frequencies, we construct a $ d$-variate trigonometric polynomial
$\displaystyle f\left(\ensuremath{\boldsymbol{x}}\right):=\sum\limits_{\ensurema...
...{\scriptsize {\rm i}}} \ensuremath{\boldsymbol{k}} \ensuremath{\boldsymbol{x}}}$    

such that $ f(\ensuremath{\boldsymbol{x}}_j) \approx y_j,\;j=0,\hdots,M-1$. Turning this into matrix vector notation, we aim to solve the system of linear equations
$\displaystyle \ensuremath{\boldsymbol{A}} \ensuremath{\boldsymbol{\hat f}} \approx \ensuremath{\boldsymbol{y}}$ (1)

for the unknown vector of Fourier coefficients $ \ensuremath{\boldsymbol{\hat f}}:=(\hat f_{\ensuremath{\boldsymbol{k}}}
...math{\boldsymbol{k}}\in I_{\ensuremath{\boldsymbol{N}}}^d} \in \mathbb{C}^{N^d}$. We denote the vector of the given sample values by $ \ensuremath{\boldsymbol{y}}:=(y_{j})_{j=0,\ldots,M-1}\in \mathbb{C}^{M}$ and the nonequispaced Fourier matrix by
$\displaystyle \ensuremath{\boldsymbol{A}}:=\left({\rm e}^{ 2\pi{\mbox{\scriptsi...
...symbol{k}} \in I_{\ensuremath{\boldsymbol{N}}}^d} \in \mathbb{C}^{M\times N^d}.$    

For $ \vert I_{\ensuremath{\boldsymbol{N}}}^d\vert<M$, the linear system (1) is over-determined, so that in general the given data $ \ensuremath{\boldsymbol{y}}$ will be only approximated up to a residual $ \ensuremath{\boldsymbol{r}}:=\ensuremath{\boldsymbol{y}}-\ensuremath{\boldsymbol{A}} \ensuremath{\boldsymbol{\hat f}}$. In order to compensate for clusters in the sampling set $ {\cal X}$, it is useful to incorporate weights $ w_j> 0$ and to consider the weighted approximation problem

$\displaystyle \Vert\ensuremath{\boldsymbol{y}} - \ensuremath{\boldsymbol{A}} \e...
...ol{x}}_j)\vert^2 \stackrel{\ensuremath{\boldsymbol{\hat f}}}{\rightarrow} \min,$ (2)

where $ \ensuremath{\boldsymbol{W}}:={\rm diag}(w_j)_{j=0,\ldots,M-1}$. This library of C functions computes approximations of (2) with the CGNR method.
The algorithms are implemented by Stefan Kunis in ./solver. Related paper are

oKunis, S. and Potts, D.
Stability Results for Scattered Data Interpolation by Trigonometric Polynomials.
SIAM J. Sci. Comput. 29, 1403 - 1419,   (full paper ps, pdf)   2007