For an application I have to approximate a continuous (and hopefully smooth) positive even function that decays at infinity with a sum of sums of gaussians, preferably orthogonal ones.

That is, a function $$f(x)|f(x)=f(-x),\lim_{x\to \infty}f(x)=0$$ I'd like to approximate with $$\Sigma_ng_n(x)|g_n(x)=\Sigma_ka_{nk}e^{-(x/s_{nk})^2}$$ so that $$\int g_ng_m=\delta_{mn}$$ and $$lim_{n\to\infty}\Sigma_ng_n(x)=f(x)$$

In addition, I would prefer for each $g_n$ to have bounded, eh, "conditioning number": $$\forall n \frac{\max_k(s_{nk})}{\min_k(s_{nk})}<C$$ where $C$ does not depend on $n$.

One way to go about it is to set $$y=e^{-x^2}$$ then approximate $$h(y):=f(x)$$ with Chebyshev polynomials over $y$, and then change variables back.

However, I don't like this solution very much because it limits $s_{nk}$ to square roots of integers and therefore unlikely to be optimal. In addition, the above "conditional number" won't be bounded.

Since I'm doing this for an application the speed of convergence matters very much. Could you suggest another method for producing the above approximation or a reference that considers such problems?

EDIT: the above seems to be related to Radial basis function approximations.