In this paper, we propose a method to approximate the Gaussian function on \(\mathbb{R}\) by a short cosine sum. We generalize and extend the differential approximation method proposed in [6, 8] to approximate \(\mathrm{e}^{-t^{2}/2\sigma}\) in the weighted space \(L^{2}(\mathbb{R}, \mathrm{e}^{-t^{2}/2\rho})\) where \(\sigma, \, \rho > 0\) . We prove that the optimal frequency parameters \(\lambda_1, \ldots , \lambda_{N}\) for this method in the approximation problem \( \min\limits_{\lambda_{1},\ldots, \lambda_{N}, \gamma_{1}, \ldots, \gamma_{N}}\|\mathrm{e}^{-\cdot^{2}/2\sigma} - \sum_{j=1}^{N} \gamma_{j} \, \mathrm{e}^{\lambda_{j} \cdot}\|_{L^{2}(\mathbb{R}, \mathrm{e}^{-t^{2}/2\rho})}\) , are zeros of a scaled Hermite polynomial. This observation leads us to a numerically stable approximation method with low computational cost of \(\mathcal{O}(N^{3})\) operations. We derive a direct algorithm to solve this approximation problem based on a matrix pencil method for a special structured matrix. The entries of this matrix are determined by hypergeometric functions. For the weighted \(L^{2}\) -norm, we prove that the approximation error decays exponentially with respect to the length \(N\) of the sum. An exponentially decaying error in the (unweighted) \(L^{2}\) -norm is achieved using a truncated cosine sum. Our new convergence result for approximation of Gaussian functions by exponential sums of length \(N\) shows that exponential error decay rates \(e^{-cN}\) are not only achievable for complete monotone functions.
