# Investigate maxima of Gaussian integral over sphere.

Let $$\alpha>0$$ be a positive parameter and consider the function

$$f(x) = \int_{\mathbb S^{n-1}} e^{-\alpha \left\lVert x-y \right\rVert^2} dS(y)$$ for $$x \in \mathbb R^n.$$ So, since this was asked, although we integrate over the unit sphere, the function "lives" on $$\mathbb R^n.$$

This function is clearly rotationally symmetric.

I would like to show that the global maxima are attained at one single radius $$r$$, only.

The rotational symmetry implies that we can consider it as a one-dimensional function by choosing $$x=(x_1,0....,0)$$, this way the exponent simplifies to $$e^{-\alpha \left\lVert x-y \right\rVert^2}=e^{-\alpha (x_1-y_1)^2+1-y_1^2}.$$

EDIT: Thanks to some interesting comments below, one can say that the global maximum is always attained at some radius $$r \in [0,1]$$ where for small $$\alpha$$ it seems to be attained close to zero and for large $$\alpha$$ it is attained closer to one.

The question remains however why is there only one radius at which the global maximum is attained? -In fact as George Lowther points out in the comments, for $$\alpha \le n/2$$ the unique maximum is attained at $$r=0$$ which leaves the case $$\alpha >n/2$$ when this does not hold true.

• This function has no global minimum. It's positive everywhere and $\to 0$ at $\infty.$ Perhaps you meant "global maximum"? – zhw. Nov 1 '18 at 19:20
• thank you, perhaps that is why nobody ever answered me :-) – Sascha Nov 1 '18 at 20:07
• Is $\lVert x \rVert = r$ and is the radius of the $(n - 1)$-sphere equal to $1$? – Maxim Nov 3 '18 at 15:23
• it is the unit sphere, so radius of the sphere is equal to $1$, but $x \in \mathbb R^n$. – Sascha Nov 3 '18 at 15:39
• Then what is $r$? – Maxim Nov 3 '18 at 16:16

Let us write $$\tilde{f}(x) = \frac{1}{|\mathbf{S}^{n-1}|}\int_{\mathbf{S}^{n-1}}e^{-\alpha|x-\omega|^2}dS(\omega)$$. Since $$\tilde{f}$$ is radially symetric, we may assume $$x = re$$ where $$r = |x|$$ and $$e = (1,0,\ldots,0)$$. Then we may write

$$\tilde{f}(re) = \frac{1}{|\mathbf{S}^{n-1}|}\int_{\mathbf{S}^{n-1}}e^{-\alpha(r^2-2r\langle e,\omega\rangle+1)}dS(\omega) = \frac{e^{-\alpha(r^2+1)} }{|\mathbf{S}^{n-1}|}\int_{\mathbf{S}^{n-1}}e^{2\alpha r\langle e,\omega\rangle}dS(\omega).$$
Disregard the factor $$e^{-\alpha}$$ and define $$f(r) = \frac{e^{-\alpha r^2} }{|\mathbf{S}^{n-1}|}\int_{\mathbf{S}^{n-1}}e^{2\alpha r\langle e,\omega\rangle}dS(\omega)$$. We can observe that $$X: \omega \mapsto \langle e,\omega\rangle\in [-1,1]$$ is a random variable with probability density function $$p_n(x) = c_n (1-x^2)^{\frac{n-3}{2}}$$ on $$[-1,1].$$($$c_n$$ is a positive constant such that $$\int_{[-1,1]}p_n(x)dx =1$$.)

Thus we have $$f(r) = e^{-\alpha r^2}\int_{-1}^{1} e^{2\alpha r x}p_n(x)dx\;(*)$$ and note that $$M_X(t) = \int_{-1}^{1} e^{t x}p_n(x)dx$$ is a moment-generating function(mgf) of $$X$$ , and $$C_X(t) = \log M_X(t)$$ is cumulant-generating function of $$X$$.

By taking logarithm on both sides of $$(*)$$, we maximize the following functional $$-\alpha r^2 + C_X(2\alpha r).$$

The first-order condition for this is:

$$r = C'_X(2\alpha r) = \frac{M'_X(2\alpha r)}{M_X(2\alpha r)}\quad \cdots (**).$$ Since $$r=0$$ is a trivial solution of $$(**)$$ we assume $$r>0$$.

What follows next is a probabilistic interpretation of this equation. Firstly, $$C'_X(t) = \frac{M'_X(t)}{M_X(t)} = E\left[X \frac{e^{tX}}{M_X(t)}\right] = \int x\cdot \frac{e^{tx}p(x)}{\int e^{ty}p(y)dy}dx$$ is an expected value of another random variable $$X_t$$ with the new density $$\frac{e^{tx}p(x)}{\int e^{ty}p(y)dy}$$. Thus $$r = C'_X(2\alpha r)$$ should be less than 1.

Secondly, $$C''_X(t) = \frac{M''_X(t)}{M_X(t)}- \left(\frac{M'_X(t)}{M_X(t)}\right)^2$$ in the similar way can be seen as variance of $$X_t$$ and thus is greater than $$0$$ (not zero since the distribution is not degenerate.) So $$C'_X(t)$$ is a strictly increasing function. Moreover, it should be that $$C'_X(t) \uparrow 1$$ as $$t\to \infty$$ since the distribution of $$X_t$$ converges to that of degenerate random variable $$X\equiv 1$$.

From the above argument, we see that there is unique $$\alpha= A(r) := \frac{(C'_X)^{-1}(r)}{2r}$$ for each $$0 satisfying $$(**)$$. The uniqueness of $$r$$'s corresponding to each $$\alpha$$ follows once if we show that $$A(r)$$ is strictly increasing on $$(0,1)$$. Since $$C'_X$$ is an increasing bijection from $$(0,\infty)$$ to $$(0,1)$$, we can substitute $$r$$ for $$r = C'_X(u)$$, and showing $$1/A(r) = \frac{2C'_X(u)}{u}$$ is decreasing in $$u\in(0,\infty)$$ is okay.

Our (almost) final claim is that $$r \mapsto C'_X(r)/r$$ is a decreasing function. Let us write by change of variable $$u = (x+1)/2$$, $$M_X(r) = e^{-r}\int_0^1 e^{2ru} \frac{u^{\frac{n-3}{2}}(1-u)^{\frac{n-3}{2}}}{B((n-1)/2,(n-1)/2)}du =: e^{-r}\Phi(2r).$$(Here, $$B(x,y)$$ is a beta function and $$\Phi$$ is a mgf of the beta distribution.)
We can calculate $$\Phi(r)$$ explicitly and omitting the procedure , we have

$$\Phi(r) =\sum_{k=0}^{\infty} \frac{a_k}{k!}r^k,\quad a_k := \prod_{j=0}^{k-1} \frac{\frac{n-1}{2}+j}{n-1+j}.$$ (https://en.wikipedia.org/wiki/Beta_distribution) From this, we have $$C'_X(r/2) = 2\left(\frac{\Phi'(r)}{\Phi(r)}-1/2\right).$$ Some more labor yields:

$$\Phi'(r) - \frac{1}{2}\Phi(r) = \sum_{k\geq 0} (a_{k+1}-a_k/2)\frac{r^k}{k!} = r \sum_{k=0}^{\infty}(a_{k+2}-a_{k+1}/2)\frac{r^k}{(k+1)!},$$ $$\frac{1}{4}\frac{C'_X(r/2)}{r/2} = \frac{\sum_{k=0}^{\infty}(a_{k+2}-a_{k+1}/2)\frac{r^k}{(k+1)!}}{\sum_{k=0}^{\infty} a_k\frac{r^k}{k!}}.$$

Using $$\frac{a_{k+1}}{a_k} = \frac{1}{2}(1 + \frac{k}{n+k-1})$$, we have $$\sum_{k=0}^{\infty}(a_{k+2}-a_{k+1}/2)\frac{r^k}{(k+1)!} = \sum_{k=0}^{\infty}\frac{a_{k+1}}{n+k}\frac{r^k}{k!}.$$ Finally, $$\frac{d}{dr}\left(\frac{1}{4}\frac{C'_X(r/2)}{r/2}\right) = \frac{\sum_{k=0}^{\infty}\frac{a_{k+2}}{n+k+1}\frac{r^k}{k!}\sum_{k=0}^{\infty} a_k\frac{r^k}{k!} -\sum_{k=0}^{\infty}\frac{a_{k+1}}{n+k}\frac{r^k}{k!}\sum_{k=0}^{\infty} a_{k+1}\frac{r^k}{k!} }{\left(\sum_{k=0}^{\infty} a_k\frac{r^k}{k!}\right)^2}.$$ Now the denominator is itself a power series whose $$s$$-th coefficient $$c_s$$is given by $$c_s = \sum_{j=0}^s \left(\frac{a_{j+2} a_{s-j}}{n+j+1}-\frac{a_{j+1} a_{s-j+1}}{n+j}\right)\frac{1}{j!(s-j)!}.$$ Using the recursion formula of $$a_k$$, $$\begin{multline} c_s= -\sum_{j=0}^s \frac{(s-2j)(n-1)+2(s-j)}{(n+j)(n+j+1)(n+2s-2j-1)}\frac{a_{j+1} a_{s-j+1}}{j!(s-j)!}\\ = -\frac{1}{2}\sum_{j=0}^s \left [\frac{(s-2j)(n-1)+2(s-j)}{(n+j)(n+j+1)(n+2s-2j-1)} +\frac{(2j-s)(n-1)+2j}{(n+2j-1)(n+s-j)(n+s-j+1)} \right]\frac{a_{j+1} a_{s-j+1}}{j!(s-j)!}. \end{multline}$$ And, $$\begin{multline}\frac{(s-2j)(n-1)+2(s-j)}{(n+j)(n+j+1)(n+2s-2j-1)} +\frac{(2j-s)(n-1)+2j}{(n+2j-1)(n+s-j)(n+s-j+1)} \\ =(n-1)(s-2j)\left[ \frac{1}{(n+j)(n+j+1)(n+2s-2j-1)} - \frac{1}{(n+2j-1)(n+s-j)(n+s-j+1)}\right] + \text{positive terms}\\ >\frac{(n-1)(s-2j)^2\{2j(s-j) +(n-1)(s-3) -4\}}{(n+j)(n+j+1)(n+2j-1)(n+s-j)(n+s-j+1)(n+2s-2j-1)} \geq 0, \end{multline}$$ if $$s\geq 3, 1\leq j\leq s-1$$.

In case $$j=0$$ or $$j=s$$, it is quite easy to check that $$\frac{s}{n(n+2s-1)} - \frac{s}{(n+s)(n+s+1)}\geq 0.$$ So, negativity of $$c_s, s\geq 3$$ is established. Some more labors yield: $$c_0 = 0, c_1 = c_2 = -\frac{1}{4n^2(n+2)}.$$

Since all $$c_s \leq 0$$, the denominator is less than 0, and the claim $$\frac{d}{dr}\left(\frac{C'_X(r)}{r}\right)<0$$ is proved.
Finally, observe that $$\lim_{r\to 0^+}A(r) = \lim_{u\to 0^+}\frac{u}{2C'_X(u)} = \frac{1}{2C''_X(0)} = \frac{n}{2}.$$ This is because $$C''_X(0)$$ is a variance of $$X=\langle e,\omega\rangle=\omega_1$$, and this value is equal to $$\frac{1}{|\mathbf{S}^{n-1}|}\int_{\mathbf{S}^{n-1}}\omega_1^2 dS(\omega)=\frac{1}{|\mathbf{S}^{n-1}|}\int_{\mathbf{S}^{n-1}}\frac{\omega_1^2 +\cdots \omega_n^2}{n} dS(\omega) = \frac{1}{n}.$$ Hence, $$A((0,1)) = (\frac{n}{2}, \infty)$$ and we conclude that:
(1) if $$0<\alpha \leq \frac{n}{2}$$, then the unique global maximizer is $$r=0$$.
(2) if $$\alpha >\frac{n}{2}$$, then the unique global maximizer is $$r = A^{-1}(\alpha)>0$$. Here, $$0$$ is not a maximizer even though it satisfies the FOC, since the second-order condition is not satisfied: $$(-\alpha r^2 + C_X(2\alpha r))''\Big|_{r=0}=2\alpha(\frac{2\alpha}{n}-1) >0.$$

• If the dimension of $S^{n - 1}$ is $n - 1$, the exponent in $p_n(x)$ is off by $1/2$ ($p_2$ cannot be a constant). – Maxim Nov 6 '18 at 15:52
• Thank you for your valuable comment. You're right. I should have been more careful. – Song Nov 7 '18 at 1:52

If $$\mathcal S$$ is the unit $$(n - 1)$$-sphere around the origin and $$\lVert x \rVert = r$$, then choosing hyperspherical coordinates s.t. $$x_1 = r, \,y_1 = \cos \theta$$ gives $$f(r) = \int_{\mathcal S} e^{-\alpha \lVert x - y \rVert^2} dS(y) = C_n \int_0^\pi e^{-\alpha (1 + r^2 - 2 r \cos \theta)} \sin^{n - 2} \theta \,d\theta = \\ C_n \sqrt \pi \,\Gamma {\left( \frac {n - 1} 2 \right)} (\alpha r)^{-n/2 + 1} e^{-\alpha (1 + r^2)} I_{n/2 - 1}(2 \alpha r), \\ f'(r) = 2 C_n \sqrt \pi \,\Gamma {\left( \frac {n - 1} 2 \right)} \alpha (\alpha r)^{-n/2 + 1} e^{-\alpha (1 + r^2)} (I_{n/2}(2 \alpha r) - r I_{n/2 - 1}(2 \alpha r)).$$ The goal is to prove that the ratio $$\omega(r) = r I_{n/2 - 1}(2 \alpha r)/I_{n/2}(2 \alpha r)$$ is monotonous in $$r$$. Since $$\omega'(r) = 2 \alpha r \left( 1 - \frac {I_{n/2 - 1}(2 \alpha r) I_{n/2 + 1}(2 \alpha r)} {I_{n/2}^2(2 \alpha r)} \right),$$ we need to prove $$I_\nu^2(z) > I_{\nu - 1}(z) I_{\nu + 1}(z), \quad\nu > 0, \,z > 0.$$ Multiplying the series expansions yields $$I_\nu^2(z) = \sum_{k \geq 0} \frac {\Gamma {\left( k + \nu + \frac 1 2 \right)}} {\sqrt \pi \,\Gamma(k + 1) \Gamma(k + \nu + 1) \Gamma(k + 2 \nu + 1)} z^{2 (k + \nu)}, \\ I_{\nu - 1}(z) I_{\nu + 1}(z) = \sum_{k \geq 0} \frac {\Gamma \!\left( k + \nu + \frac 1 2 \right)} {\sqrt \pi \,\Gamma(k + 1) \Gamma(k + \nu + 1) \Gamma(k + 2 \nu + 1)} \frac {k + \nu} {k + \nu + 1} z^{2 (k + \nu)}, \\ \frac {[z^{2 (k + \nu)}] I_\nu^2(z)} {[z^{2 (k + \nu)}] I_{\nu - 1}(z) I_{\nu + 1}(z)} = 1 + \frac 1 {k + \nu},$$ proving the inequality. Since $$\omega(\infty) = \infty$$, $$\omega$$ increases to infinity monotonously.

Finally, $$\omega(0) = n/(2 \alpha)$$. If $$\omega(0) < 1$$, there is a unique $$r_0 > 0$$ where $$\omega(r_0) = 1$$ and $$f'(r_0) = 0$$, therefore $$f(r)$$ increases on $$[0, r_0]$$ and decreases on $$[r_0, \infty)$$. If $$\omega(0) \geq 1$$, $$f(r)$$ decreases on $$[0, \infty)$$.

• sorry, what precisely is this constant $C_n$? – Sascha Nov 5 '18 at 17:51
• It's the integral of the area element over the remaining $n - 2$ hyperspherical angles, thus the surface area of a unit $(n - 2)$-sphere. – Maxim Nov 5 '18 at 18:12
• thank you very much. – Sascha Nov 5 '18 at 18:13

This is certainly not a full solution.

First we would like to get an expression for $$f(x)$$. (Below I write $$f(x)$$ or $$f_n(x)$$ when I should probably should have written $$f((x_1, 0,0,\ldots, 0)$$.)

I am going to use the function $$s(n,r) := \frac{2 \pi ^{\frac{n+1}{2}}}{\Gamma \left(\frac{n+1}{2}\right)} r^n,$$ the surface area of an n-dimensional sphere.

We can use $$s$$ to integrate $$f$$,

$$f_n(x) = \int_{\theta=0}^\pi s (n-1 ,\sin(\theta)) \exp\left(-a \left((x-\cos (\theta))^2+\sin ^2(\theta)\right)\right)d\theta.$$

According to Mathematica, $$\;f_2(x) = \frac{\pi e^{-a (x+1)^2} \left(e^{4 a x}-1\right)}{a x}$$.

And $$f_4(x) = \frac{\pi ^2 e^{-a (x+1)^2} \left(2 a x+e^{4 a x} (2 a x-1)+1\right)}{2 a^3 x^3}$$.

For odd values of $$n$$, Mathematica gives Hypergeometric functions. $$f_3(x) = 2 \pi ^2 e^{-a \left(x^2+1\right)} \, _0\tilde{F}_1\left(;2;a^2 x^2\right)$$

In the Mathematica language, $$f_3(x)$$ is

"(2*Pi^2*Hypergeometric0F1Regularized[2, a^2*x^2])/E^(a*(1 + x^2))".

If you try to find the maximum of $$f$$, it seems that the maximum occurs at $$x=0$$ when $$a<2$$. One way to think about this is releasing a million drunkards randomly on the $$S^{n-1}$$ sphere and having them walk randomly. Initially, the most drunkards per unit volume occurs on the $$S^{n-1}$$ sphere, but as time goes on the radius of maximum drunkards moves toward the center.

Another way to look at it is as a convolution of the sphere with a Gaussian with $$\sigma = \sqrt{1/(2 a)}$$. When sigma is large, $$f_n(x)$$ is approximately Gaussian.

I was not able to find closed form solutions for $$\mathrm{argmax}_x\, f_n(x)$$. If we take the derivative when $$n$$ is even, set it equal to zero, and simplify, we get expressions of the form $$p_1(x,a) + \cosh(2 a x)p_2(x,a) +\sinh(2 a x)p_2(x,a) =0$$ where the $$p_i$$ are polynomials of degree of low degree with "nice" coefficients.

For example, for $$n=10$$, Mathematica says that the derivative of $$f_{10}$$ is zero if $$32 a^5 x^6 \sinh (2 a x)-32 a^5 x^5 \cosh (2 a x)-160 a^4 x^5 \cosh (2 a x)+240 a^4 x^4 \sinh (2 a x)+360 a^3 x^4 \sinh (2 a x)-840 a^3 x^3 \cosh (2 a x)-420 a^2 x^3 \cosh (2 a x)+1680 a^2 x^2 \sinh (2 a x)+210 a x^2 \sinh (2 a x)+945 \sinh (2 a x)-1890 a x \cosh (2 a x)=0.$$

Hope that helps.

• yes, I feel that this does not tell us much about the location and number of global maxima so far. The hypergeometric function is still a somewhat opaque object. – Sascha Nov 3 '18 at 18:24