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

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$$.

• (1) Could you explain what "completely reconstruct $f$" might mean, given that you cannot hope to estimate its parameters exactly with any finite sample (and, without a positive lower bound on the $p_i,$ you can never be sure you have identified all the components unless you know $N$)? (2) By analyzing the absolute values of the data you don't need the symmetry assumption and $N$ can be reduced to $\lceil N/2 \rceil.$ Moreover, when $N$ is known and odd you already know one of the $s_i=0.$
– whuber
Jul 11, 2022 at 14:43
• @whuber Regarding (1) we want $\sup_x |\hat{f} (x)-f(x)|=0$ where $\hat{f}$ is the reconstruction of $f$. I guess I wasn't specific about knowledge of $N$. We assume that we know $N$. I will add this.
– Boby
Jul 11, 2022 at 14:51
• @Xi'an yes, we know that values of $f(x_1), ..., f(x_k)$ exactly.
– Boby
Jul 11, 2022 at 15:00
• You should then remove mentions of sample and * sampling* as this is confusing, favouring the impression that this is a statistical problem. (I am also unsure the assumption that $p(s)=p(-s)$ is helping. I would first consider the problem without this condition. Jul 11, 2022 at 15:09
• Ok. I will, but I felt this was standard terminology. We have something called the Nyquist sampling theorem. This is along the same lines. What do you think I should use instead of sampling?
– Boby
Jul 11, 2022 at 15:11

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).

• The idea is good, but what happened to the $\exp(-s_i^2/2)$ terms? As a practical matter, if $2C\gg N,$ it won't be possible to evaluate $f$ with sufficient accuracy unless many more than $N$ evaluations are done, making this solution a mathematical curiosity rather than a practical approach. As a practical matter, using double precision floats one would have to sample a bare minimum of $C/8$ equally spaced locations to be sure of detecting a contribution from every component; and to obtain reasonable (single float) precision in the reconstruction, the count needs to be at least $C/6.$
– whuber
Jul 11, 2022 at 16:26
• A strategic choice of the $x_i$ requires the $s_j$'s to be known to some extent. Otherwise, as stated by @whuber, sampling every fourth integer between $-C$ and $C$ would be needed for each component to be detectable (since outside $(s_j-2,s_j+2)$ the component density is numerically [almost] zero). Jul 11, 2022 at 17:57
• Thank you for pointing out what I should have noticed! I guess I view the dependency on $s_i$ through $\exp(-s_i^2/2)$ to be so critical to analyzing the problem that I needed to see it explicitly. But for your argument it's unnecessary. Indeed, all that is needed is to point out there's a system of $K$ equations of the form $f(p_i,s_i)=0$ to be solved for $N$ values $(p_i,s_i),$ so that generically you need $N$ such equations (by symmetry). But if some are redundant you need more; and in very special cases (I believe not possible here) a smaller set of equations can yield a unique solution.
– whuber
Jul 11, 2022 at 18:39
• BTW, some experiments show that as soon as $N$ is as large as $3,$ the condition number of derivative is essentially infinite (in double precision). Thus, as a practical matter, the problem appears hopeless except in the simplest situations.
– whuber
Jul 12, 2022 at 15:24
• In that case, note that $2K-1$ points will generically work due to the sum-to-unity constraint on the $p_i.$
– whuber
Jul 13, 2022 at 2:12

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.

• The function can be exactly found, (purely in mathematical theory--not as a practical matter!) by sampling it at practically any set of $2N-1$ distinct points. As a practical matter, the problem is nearly impossible with any $N\gt 3$ unless you are working in very high precision arithmetic.
– whuber
Jul 25, 2022 at 17:53