Global maximum of a sum of Gaussians

Question

Suppose I have the weighted sum of $n$ Gaussian functions with varying means and standard deviations, $$f(x) = \sum_{i=1}^n a_i\exp\left(-\frac{(x-\mu_i)^2}{2\sigma_i^2}\right).$$ All the weights $a_i$ are positive. I want to find the position of the global maximum of $f$.

Clearly $f$ is non-concave and may have multiple local maxima. However, it feels like each local maximum should be close to the peak $\mu_i$ of one of the Gaussians.

Suppose I perform gradient ascent $n$ times starting from $x = \mu_i$ for $i = 1, 2,\ldots, n$. Am I guaranteed to find all the local maxima, and thus the global maximum? Is there a better strategy to find the global maximum?

Answer 1

Here's a partial result making the intuition in my question rigorous.

Let $f(x) = \sum\limits_{i=1}^n g_i(x)$, where $$g_i(x) = a_i\exp\left(-\frac{(x-\mu_i)^2}{2\sigma_i^2}\right)$$ are the $n$ Gaussians being summed together.

Suppose $x$ is a critical point of $f$, that is, $f'(x)=0$. For it to be a local maximum, $f''(x)$ must be negative. But $f''(x)=\sum\limits_{i=1}^ng_i''(x)$, and $$g_i''(x)=\frac{(x-\mu_i)^2-\sigma_i^2}{\sigma_i^4}g(x)$$ is negative only when $x\in[\mu_i-\sigma_i,\mu_i+\sigma_i]$. So any local maximum can only occur within $\sigma_i$ of a $\mu_i$.

This allows us to only search within the potentially smaller union of intervals $\bigcap\limits_{i=1}^n [\mu_i-\sigma_i,\mu_i+\sigma_i]$ rather than all of the approximately $[\min\mu_i,\max\mu_i]$ as I interpret Yuandong's answer. This is helpful if the $\mu_i$ are widely separated. Nevertheless, $f$ is not guaranteed to be concave in any of these intervals, so we cannot be assured to find a global maximum simply by gradient ascent. One can get within $\epsilon$ of the global maximum by a grid search with step size $\epsilon/M$, where $M$ is as defined in Yuandong's answer, but it would be nice to know if any faster convergence is possible.

Answer 2

Since $x$ is one dimensional, you can do a simple grid search with fine enough grids, which guarantees to find global optimum. It is not necessary the case that each local maximum is very close to $\mu_i$. In fact, new local maxima can be created between two distant $\mu_i$s if their $\sigma_i$s are large enough.

EDIT

To see how the grid search will give guarantees on global optimal solution. We first compute the bound of $f'(x)$:

$$|f'(x)| \le \sum_i a_i \left|\frac{x-\mu_i}{\sigma_i^2}\right|\exp\left(-\frac{(x-\mu_i)^2}{2\sigma_i^2}\right)$$

since (by taking derivative)

$$\frac{|x|}{\sigma^2}\exp\left(-\frac{x^2}{2\sigma^2}\right) \le \frac{1}{\sigma}e^{-1/2}$$

we have $$|f'(x)| \le e^{-1/2}\sum_i \frac{a_i}{\sigma_i} \equiv M$$

On the other hand, there exists a sufficiently large $K$ so that $|f(x)| \le \max_i a_i$ for $x \notin [-K,K]$ and cannot be the global optimum solution. We thus divide $[-K, K]$ into $\lceil 2KM/\epsilon \rceil$ bins $[l_j, u_j]$ so that $u_j - l_j \le \epsilon/M$. Thus since $|f'(x)|\le M$, within each bin $j$ there exists $c_j$ so that

$$c_j - \frac{\epsilon}{2}\le f(x) \le c_j + \frac{\epsilon}{2}$$

Let $j^* = \arg\max_j c_j$. Then any $\hat x^* \in [l_{j^*}, u_{j^*}]$ gives $f(\hat x^*)$ that is at most $\epsilon$ worse than the true global optimum.

I admit that from this analysis, it is not possible to provide a finite time algorithm that find the exact global optimum.