Logarithmic gradient of an exponential squared-distance sum (source code)

= Logarithmic gradient of an exponential squared-distance sum
{title2=$\nabla\log\sum_i e^{\|x-a_i\|^2}=2(x-\sum_iw_i a_i)$}

For fixed <vectors> $a_1,\ldots,a_m$ in a finite-dimensional <Euclidean space>, define
$$
F(x)=\sum_{i=1}^m e^{\|x-a_i\|^2},\quad
w_i(x)=\frac{e^{\|x-a_i\|^2}}{F(x)},\quad
\bar a(x)=\sum_iw_i(x)a_i.
$$
The <gradient> and <Hessian matrix> of $\Phi=\log F$ are
$$
\nabla\Phi(x)=2(x-\bar a),\qquad
\nabla^2\Phi(x)=2I+4\sum_iw_i(a_i-\bar a)(a_i-\bar a)^T.
$$
To derive the first identity, differentiate each <exponential function> and divide the sum of its <gradients> by $F$. For the second, differentiating a normalized weight gives $\nabla w_i=2w_i(\bar a-a_i)$. Thus $J_{\bar a}=-2\operatorname{Cov}_w(a_i)$, where the weighted <covariance matrix> is the sum in the displayed <Hessian matrix>. Differentiating $2(x-\bar a)$ completes the calculation.

The weighted <covariance matrix> is positive semidefinite, so $\nabla^2\Phi\succeq2I$. Also $\Phi(x)\geq\|x-a_1\|^2\to\infty$ as $\|x\|\to\infty$. Therefore $F$ and $\Phi$ share a unique global minimizer, characterized by $x=\bar a(x)$. The weights place that minimizer in the <convex hull> of the centres.

For stable numerical evaluation, subtract the largest exponent before computing the normalized weights. This is the <log-sum-exp> principle and avoids <floating-point overflow> without changing the <gradient>. A naive iteration $x_{n+1}=\bar a(x_n)$ need not converge: in one dimension, centres $-1,1$ give $\bar a(x)=-\tanh(2x)$ with <derivative> $-2$ at its unique fixed point. A <Newton root-finding iteration> or a suitable <Quasi-Newton method> instead uses the positive curvature of $\Phi$.