 Solution

A practical choice is a BFGS method, using exact gradients and a line search. At dimension $d=100$, its $d$-by-$d$ inverse-Hessian matrix approximation is small, whereas forming a million rank-one contributions to the exact Hessian matrix is expensive. If $F(x)={\sum }_i{e}^{\Vert x-a_i{\Vert }^2}$, the explicit derivatives areThe positive-definite Hessian matrix and coercivity show that the objective has one global minimizer.

For numerical stability one can minimize $\Phi =\log F$ instead, which has exactly the same minimizer. Put $r_i=\Vert x-a_i{\Vert }^2$, $w_i=e^{r_i}/{\sum }_j{e}^{r_j}$, and $\bar{a}={\sum }_i{w}_i{a}_i$. Then the logarithmic gradient of an exponential squared-distance sum givesEvaluate the weights as $e^{r_i-r_{\max }}/{\sum }_j{e}^{r_j-r_{\max }}$ to avoid floating-point overflow. A Quasi-Newton method step iswhere $H_n$ is the positive-definite inverse-Hessian matrix approximation updated by BFGS, and ${\alpha }_n$ is chosen by a line search, for example one satisfying the Wolfe conditions. Under the standard smoothness and nondegeneracy assumptions at this strongly convex minimizer, full-memory BFGS has local superlinear convergence. An exact Newton's method in optimization would have local quadratic convergence, but pays for those much more expensive exact second derivatives.^[ref]