Approximation of a gaussian function
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.