Find maximum Gaussian function that fits under a curve

Question

I have a bounded, positive function $f(x)$, and I would like to fit a single Gaussian function $g(x)$ "underneath" it such that

$$g(x) \leq f(x) \qquad \forall x$$

while maximising the integral under the Gaussian. (Either in closed-form or by iterative optimisation.) i.e. maximise

$$ \int_\mathbb{R}g(x)dx $$

subject to the above constraint. We could express $g$ as:

$$g(x) = a e^\frac{-(x-b)^2}{2c^2}$$

and in that form, we're seeking to maximise $\int_\mathbb{R}g(x)dx = \sqrt{2 \pi} a |c|$.

Note that $g$ is a Gaussian function but not necessarily a Gaussian distribution, i.e. not constrained to integrate to 1.

In my case $f(x)$ is a Gaussian mixture model, though perhaps there's an approach that works for general $f(x)$. I'd also hope there's an approach that works for the multivariate equivalent of this problem too.

---

P.S. here's a plot showing some arbitrary examples, in this case using "step 1" from @mark-l-stone 's approach.

Answer 1

The complication for this optimization problem is that this is an infinite dimensional problem due to constraint $g(x) \le f(x)$ needing to hold for all x.

Because infinite dimensional optimization solvers basically don't exist, the only way i can think of to handle this problem is imposing constraints at some finite number of values of x in lieu of the infinite dimensional constraint, then performing

Step 1: Numerically maximize, with respect to $\{a,b,c\}$, the integral of g(x), which is proportional to $ac$, and which can be taken as the objective function, subject to the constraint that $c$ be nonnegative, and to the additional constraints that $g(x) \le f(x)$, at each of a finite number of specified values of x.

Step 2: Numerically globally maximize $g^*(x) - f(x)$, with respect to x over x from $-\infty$ to $\infty$, where $g^*(x)$ is g(x) using the optimal parameters found in Step 1.

If the optimal objective value of the step 2 optimization is nonpositive, then the step 1 optimization solved the original problem.

Else Add constraints for additional values of x, and re-run Step 1.Repeat as necessary.

This paradigm could be used regardless of whether f(x) is a Gaussian mixture density.

In principle, this same approach could be applied to multivariate g(x) and f(x). however, the computational challenge is likely to increase substantially. The curse of dimensionality will be going against you in instantiating an adequate, but tractable, set of points at which to impose constraints for step 1.

Answer 2

Seeing as I can't think of anything better, I'll start with a somewhat naive approach. First, notation. Let $N[\lambda,A](x)=A\exp(-\lambda x^2)$ where $\lambda, A>0$. The first step is to find the maximum value that $f$ attains, and the $x$ value this maximum occurs at. Call this point $(x_*,f_\max)$. This assumes two things. 1: $f$ has a maximum at all 2: $f$ has a unique maximum. If $f$ doesn't have a maximum, I don't think there is any limit to how large we can make $\int_{\Bbb{R}} g(x)\mathrm{d}x$. But, I don't yet have a proof of this. But, if $f$ does not have a unique maximum, i.e, it attains the value of $f_\max$ at multiple points, this problem will be much more difficult. So let's assume $f$ does have a unique maximum. So, given that this maximum occurs at $(x_*,f_\max)$, we need to choose the parameters $A, \lambda$ to maximize the integral. It is fairly elementary to see that a choice of $A=f_\max$ is optimal. But what about $\lambda$? In lieu of a better method, we can simply use a midpoint search algorithm. Pick some suitably small value for $\lambda$, say $0.01$, and some suitably large value, say $100$. Now find the midpoint between these two values ($\approx 50$) and check if the equation $f(x)=N[\lambda,f_\max](x)$ has any roots in $\Bbb{R}$. If it does not, restrict the search area to $50\leq \lambda \leq 100$ and repeat the process. If it does, restrict the search area to $0.01 \leq \lambda \leq 50$ and repeat the process. Do this several times, and the process to converge to an optimal $\lambda_\max$. This is arduous with pen and paper, but should be easy enough with a Python program. Once you've found this value, the integral maximizing gaussian should be $$g(x)=N[\lambda_\max, f_\max](x-x_*).$$ If anyone can offer any refinements on this algorithm, do let me know in the comments!