Least squares
Let \({\boldsymbol{N}}=(N_1,\ldots,N_d)^{\top}\in 2\mathbb{N}^d\). For given samples \((\boldsymbol{x}_j,y_j)\in \mathbb{T}^d\times\mathbb{C},\; j=0,\dots,M-1\), the index set \(I_{\boldsymbol{N}}^d:= \big\{-\frac{N_1}{2},\dots,\frac{N_1}{2}-1\big\}\times \dots\times\big\{-\frac{N_d}{2},\dots,\frac{N_d}{2}-1\big\}\), of frequencies, we construct a \( d\)-variate trigonometric polynomial \[ f\left({\boldsymbol{x}}\right):=\sum\limits_{\boldsymbol{k} \in I_{{\boldsymbol{N}}}^d} \hat f_{\boldsymbol{k}} {\rm e}^{ 2\pi\mathrm i \boldsymbol{k} \boldsymbol{x}} \] such that \(f(\boldsymbol{x}_j) \approx y_j,\;j=0,\dots,M-1\). Turning this into matrix vector notation, we aim to solve the system of linear equations \[ \boldsymbol{A} {\boldsymbol{\hat f}} \approx {\boldsymbol{y}} \tag{1} \] for the unknown vector of Fourier coefficients \({\boldsymbol{\hat f}}:=(\hat f_{{\boldsymbol{k}}} )_{{\boldsymbol{k}}\in I_{{\boldsymbol{N}}}^d} \in \mathbb{C}^{N^d}\). We denote the vector of the given sample values by \(\boldsymbol{y}:=(y_{j})_{j=0,\ldots,M-1}\in \mathbb{C}^{M}\) and the nonequispaced Fourier matrix by \[ {\boldsymbol{A}}:=\left({\rm e}^{ 2\pi{\mathrm i} {\boldsymbol{k}} {\boldsymbol{x}}_j}\right)_{j=0,\dots,M-1;{\boldsymbol{k}} \in I_{{\boldsymbol{N}}}^d} \in \mathbb{C}^{M\times N^d}. \] For \(|I_{\boldsymbol{N}}^d|\lt M\), the linear system (1) is over-determined, so that in general the given data \(\boldsymbol{y}\) will be only approximated up to a residual \({\boldsymbol{r}}:={\boldsymbol{y}}-{\boldsymbol{A}} {\boldsymbol{\hat f}}\). In order to compensate for clusters in the sampling set \( {\mathcal X}\), it is useful to incorporate weights \( w_j> 0\) and to consider the weighted approximation problem \[ \|{\boldsymbol{y}} - {\boldsymbol{A}} {\boldsymbol{\hat f}}\|_{{\boldsymbol{W}}}^2 = \sum_{j=0}^{M-1} w_j |y_j-f({\boldsymbol{x}}_j)|^2 \stackrel{{\boldsymbol{\hat f}}}{\rightarrow} \min, \tag{2} \] where \(\boldsymbol{W}:={\rm diag}(w_j)_{j=0,\dots,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