Interpolation
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}. \]
We focus on the under-determined and consistent linear system \({\boldsymbol{A}} {\boldsymbol{\hat f}} = {\boldsymbol{y}}\), i.e., we expect to interpolate the given data \(y_j\in\mathbb{C}\), \(j=0, \dots,M-1\), exactly.
In particular, we incorporate damping factors \(\hat w_{{\boldsymbol{k}}}> 0\), \({\boldsymbol{k}} \in I_{{\boldsymbol{N}}}^d\), and consider the optimal interpolation problem \[ \|{\boldsymbol{\hat f}}\|_{{\boldsymbol{\hat W^{-1}}}}^2 =\sum_{{\boldsymbol{k}}\in I_{{\boldsymbol{N}}}^d} \frac{|\hat f_{{\boldsymbol{k}}}|^2}{\hat w_{{\boldsymbol{k}}}} \stackrel{{\boldsymbol{\hat f}}}{\rightarrow} \min \quad \text{subject to} \quad {\boldsymbol{A}} {\boldsymbol{\hat f}} = {\boldsymbol{y}}, \tag{2} \] where \({\boldsymbol{\hat W}}:= {\rm diag}(\hat w_{{\boldsymbol{k}}})_{{\boldsymbol{k}}\in I_{{\boldsymbol{N}}}^d}\).
This library of C functions computes approximations of (2) with the CGNE method.
The algorithms are implemented by Stefan Kunis in ./solver. Related paper are