Interpolating Gaussian mixture. How many points is sufficient for complete reconstruction

Question

Consider the following Gaussian mixture of $N$ components: \begin{align} f(x)= \sum_{i=1}^N p (s_i) e^{-\frac{(x-s_i)^2}{2}}\big/ \sqrt{2\pi} \end{align} where we assume that

• $\max_{i} |s_i| \le C$
• $p (s_i)$ is probability distribution (i.e., $\sum_i p(s_i) =1 $ and $1 \ge p(s_i) \ge 0 $
• $p(s_i)$ is symmetric. That is, if $s_i \neq 0$, then there exists $j$ such that $i \neq j$ and $s_i=-s_j$ and $p(s_i)=p(s_j)$.
• we know the value of $N$

Otherwise, we have no knowledge of $s_i$'s or $p(s_i)$

My question is the following: Suppose we are allowed to evaluate/sample $f$ at any $K$ points we desire. How many samples $K$ do we require to completely reconstruct $f$?

Unfortunately, $f$ is not bandlimited in the Fourier domain otherwise we could have used Nyquist criteria. However, because of the Gaussian tail, it is almost bandlimited, plus we have additional information about its structure (i.e., we know $N$). So, therefore, I am curious if $K<\infty$ is enough.

This question must have been asked before but I couldn't find any reference and don't really know for what keywords to search.

Some thoughts

One can set this up as a linear algebra problem. First, choose $K$ samples $x_i$ at which we will evaluate $f$ then and let \begin{align} F=[f(x_1),\ldots, f(x_k) ]^T \end{align} also define a matrix $E$ such that \begin{align} E_{ij}= e^{-\frac{(x_i-s_j)^2}{2}}\big/ \sqrt{2\pi} \end{align} and let \begin{align} P=[p(s_1),...p(s_N)] \end{align} Therefore, our problem can be written as \begin{align} E P=F \end{align} we are interested in recovering $E$ and $P$ based on the knowledge of $F$.

Answer 1

Since $$\sum_{i=1}^K p(s_i)e^{-(x_j-s_i)^2/2}=e^{-x_j^2/2}\sum_{i=1}^K p(s_i)e^{x_js_i-s_i^2/2}=e^{-x_j^2/2}\sum_{i=1}^K \underbrace{p(s_i)e^{-s_i^2/2}}_{q_i}e^{x_js_i}$$ it is equivalent to consider the system of equations $$\require{cancel}\cancel{e^{x_j^2/2}}\cancel{e^{-x_j^2/2}}\sum_{i=1}^K q_i\left[e^{s_i}\right]^{x_j}=\alpha_je^{x_j^2/2}\qquad j=1,\ldots,N$$ which should be solvable with $N=2K$ interpolation points $x_j$ taken at random in $(-C,C)$ (without considering the further constraints).

Answer 2

I doubt whether exact recovery of the function from a finite number of test points is possible in general. And if it is possible, it will be difficult. But there are very practical ways to bound the error and based on this error bound to devise designs of test points to minimise the error. For simplicity I will largely ignore the constraints on the coefficients and the values for $s$ as I do not think they change the basics of the problem.

The vector space spanned by the set of functions $\{f(x)=\exp(-\frac {(x-s)^2} {2})\mid s\in\mathbb{R}\}$ is a reproducing Hilbert space (RKHS) with kernel $K(x,s)=\exp(-\frac {(x-s)^2} {2})$. Since $K$ is a positive definite function, this space is infinite dimensional. Your equation $EP=F$ will always have a unique solution, no matter how many (distinct) interpolation points you choose in $F$. Furthermore, this solution can be characterised as the function with minimum Hilbert-space norm among all interpolating functions and of course by construction it is a N-component function.

But this Hilbert space is large and in general there will be other interpolating functions, some of them also with $N$-components. This can be demonstrated by a simple example with N=2. In the picture below I plot the functions $\frac{1}{2}\exp(-\frac {(x+ 1)^2}{2})+\frac{1}{2}\exp(-\frac {(x- 1)^2}{2})$ and $\frac{1}{2}\exp(-\frac {(x+ 2)^2}{2})+\frac{1}{2}\exp(-\frac {(x- 2)^2}{2})$. Note that the graphs intersect. If you happen to choose these points for interpolation you have two possible 2-component interpolants and you will not be able to distinguish between them from the test points alone.

One may argue that these intersections look like very special points, which may be avoided. But if you allow larger $N$ these functions can get more and more complicated and the intersections harder and harder to describe. The full space (without constraints) is actually universal, which means you can approximate every continuous(!) function on a compact interval in the uniform norm.

That said, if you do not require exact recovery but can live with approximation then RKHS theory provides you with pointwise error estimates. Those are based on the kernel function $K$, the fill distance of your design of test points, which is basically the largest possible distance to any of your test points and the norm of your target function $f$. For details see Chapter 14 of Gregory E. Fasshauer, "Meshfree Approximation Methods with MATLAB" and in particular Example 15.1 for the approximation order of the Gaussian kernel.

Finally, the Gaussian kernel is known to be numerically fickle. See cit.loc. Example 16.1 for a discussion of the condition number of Gaussians. This can be addressed by various methods, such as pre-conditioning and choosing a better basis. See cit.loc. Chapter 34 for a discussion of those issues.