# Approximation of a gaussian function

I want to approximate a Gaussian function as shown below $$e^{\frac{-\|x-c\|^2}{2\sigma^2}} \approx \sum_{i=1}^{N}\alpha(c,c_i)e^{\frac{-\|x-c_i\|^2}{2\sigma^2}} \forall x$$ Here c, $\sigma$ and $c_i$ are known. Are there functions $\alpha(c,c_i)$ which gives a good approximation? The limitation I am facing is that I can not optimize as I have to calculate the functions $\alpha$ for various c(around 1024*1024) thus leading to too many optimizations. I used linear solver as this optimisation problem can be framed as $A\alpha=b$. Here $x$ in the above equation is sampled according to each $c$. Hence for every $c$, both $A$ and $b$ in $A\alpha = b$ changes.

I am also fine if functions in RHS are not Gaussians i.e existence of pair of $\alpha$ and $f_\sigma$ such as

$$e^{\frac{-\|x-c\|^2}{2\sigma^2}} \approx \sum_{i=1}^{N}\alpha(c,c_i) {f_\sigma}(x,c_i) \forall x$$

Can someone help me to point out how to go about it? Here N can be treated as a variable. So as N increases number of gaussians used to interpolate increases and hence get a better approximation.

Now to explain what I meant by a good approximation : The function in LHS is used as a filtering kernel to filter my data. Replacing it with RHS (one part of my algorithm) will help in speeding up my filtering algorithm by a huge margin. Now I check my performance by comparing actual filtering algorithm and this approximation by using parameters such as Mean Square Error.

• This seems like a classical problem in interpolation. First you have to explain what you mean by "good approximation", and then open a textbook on numerical analysis or on approximation theory and find the answer. – Yuval Filmus May 22 '17 at 13:06
• 1. I suspect your notation is misleading. I suspect the best value for $\alpha(c,c_i)$ probably depends not just on $c$ and $c_i$ but all of $c_1,\dots,c_n$, so the notation should probably something more like $\alpha(c,c_1,\dots,c_n,i)$. 2. How do you plan to measure the quality of the approximation? – D.W. May 22 '17 at 16:10
• @Yuval. I use this approximation in another algorithm of data reconstruction which will help me to rate this approximation with proper metric. – pravin May 23 '17 at 10:23
• @D.W Here N do come into picture. But I don't find a reason why every interpolation function (total N of them ) should be dependent on other gaussian means as well. – pravin May 23 '17 at 10:28
• You still need to tell us what the proper metric is. We probably can't answer the question without knowing how you want to measure the quality of the approximation. (We can guess and make up a metric of our own, but if that isn't the metric that matters to you, that might be wasting our time.) Also can you edit the question to tell us more about what you tried, when you say optimizing is too slow? What optimization approach did you try? This might give us an idea whether better optimization methods might be faster. How large is $N$? – D.W. May 25 '17 at 16:35

One approach could be to use least-squares fitting.

Your problem is: given functions $f_1,\dots,f_N$ and $f$, we want to find constants $\alpha_1,\dots,\alpha_N$ such that

$$f(x) \approx \sum_i \alpha_i f_i(x).$$

Here

$$f(x) = e^{\frac{-\|x-c\|^2}{2\sigma^2}}$$

and

$$f_i(x) = e^{\frac{-\|x-c_i\|^2}{2\sigma^2}}.$$

To make this problem well-defined, we need to make precise the notion of approximation. One plausible notion is the $L_2$ error of the approximation, i.e., the error is

$$E = \int \big( f(x) - \sum_i \alpha_i f_i(x) \big)^2 \; dx.$$

This is tricky to compute, but we can approximate $E$ by picking values $x_1,\dots,x_m$ and then calculating

$$\hat{E} = \sum_j \big( f(x_j) - \sum_i \alpha_i f_i(x_j) \big)^2.$$

If the $x_j$'s are chosen appropriately (e.g., evenly spaced over the range you care about; or randomly chosen according to the distribution of $x$-values you're likely to evaluate the approximation at), and if $m$ is large enough, then $\hat{E}$ should be a good approximation to $E$.

So, the problem becomes: given functions $f_1,\dots,f_N$ and $f$ and values $x_1,\dots,x_m$, we want to find constants $\alpha_1,\dots,\alpha_N$ that minimize the error $\hat{E}$.

This turns out to be exactly a multiple linear regression problem. We try to find a fit

$$z = \alpha_1 y_1 + \dots + \alpha_N y_N$$

where $y_i = f_i(x)$ and $z = f(x)$. We can evaluate the $f_i$'s and $f$ at any input $x$ of our choice. When we evaluate at $x_j$, we obtain a single input $(y_1,\dots,y_N)$ and a single output $z$; in this way, we obtain $m$ matching input-output pairs. You can now use a standard library to solve this multiple linear regression problem and find $\alpha_1,\dots,\alpha_N$.

All that remains is to choose $m$ and $x_1,\dots,x_m$. I would suggest that you try to figure out what values you'll evaluate the approximation on, and choose the $x_j$'s randomly from the distribution (it might be a uniform distribution over some range, or something else), or choose the $x_j$'s to be evenly spaced over some range. The larger $m$ is, the better the approximation will be, but the longer the algorithm will take, so experiment with some different values of $m$ and see what the tradeoff is. My guess is you'll want $m$ to be a small constant multiple of $N$ (say $m=2N$), but you can see what you get.

Since you said $N$ is small (8--16), this should be very fast.

An optimization: Since you're going to repeat this many times with the same values for $c_1,\dots,c_N$ (i.e., $f_1,\dots,f_N$) but a different value of $c$ each time, you can also speed things up. If the matrix is $X$ and the vector is $y$, the solution is $(X^\top X)^{-1} X^\top y$, where $X$ depends only on the $c_i$'s but not on $c$. Thus, you can precompute the matrix $(X^\top X)^{-1} X^\top$ (once), then for each instance of the problem you want to solve, you just have to do a matrix multiplication -- no need to solve a new linear system.

Lastly, one optional refinement you could also try: use multiple linear regression to choose $\alpha_i$'s, then follow that up with gradient descent to fine-tune the $\alpha_i$'s and improve the quality of the approximation. Often if you choose a good starting point that is close to optimal, gradient descent will converge faster. Thus you can use the solution from multiple linear regression as the starting point for the gradient descent optimization routine.

• Thanks for this valuable answer. Sorry for the delay in replying as I wanted to answer after proper analysis. I am using the trick of optimization you mentioned. Here my ci's are non uniform samples of my data. So I used them as my samples. Hence m=N. This helps but I believe I should take more samples for better approximation. I have a doubt. If I take my whole data as my samples, then my approximation would be the best,right? – pravin Jun 1 '17 at 11:30
• I have no prior knowledge in approximation theory. Can you please help me with some references(preferably papers) which can help me solve this problem? – pravin Jun 1 '17 at 11:32
• @pravin, I apologize, but I didn't understand your first question. (I don't see any reason to expect that $x_i=c_i$ is optimal, nor that $m=N$ is necessarily best, if that's what you are asking, but you can experiment and see what happens.) I don't know of any references I can suggest. – D.W. Jun 1 '17 at 15:16
• Sorry for posing the question in unclear way. I am not expecting that m=N is best. I just experimented using m=N as xi=ci as my data as ci are non uniform samples of my data and it is clear that more samples are required for better approximation. – pravin Jun 1 '17 at 16:10
• My questions: 1) If I take my entire data set as my samples, will the approximation be the best i.e will L2 error be minimum? 2) I was asking for references in approximation theory on how to sample the data. – pravin Jun 1 '17 at 16:14