Home > examples > thomson_problem.m

# thomson_problem

## PURPOSE

Simple attempt at computing n well distributed points on a sphere in R^d.

## SYNOPSIS

function X = thomson_problem(n, d)

## DESCRIPTION

``` Simple attempt at computing n well distributed points on a sphere in R^d.

This is an example of how Manopt can approximate the gradient and even
the Hessian of a cost function based on finite differences, even if only
the cost function is specified without its derivatives.

This functionality is provided only as a help for prototyping, and should
not be used to compare algorithms in terms of computation time or
accuracy, because the underlying gradient approximation scheme is slow.

## CROSS-REFERENCE INFORMATION

This function calls:
• obliquefactory Returns a manifold struct to optimize over matrices w/ unit-norm columns.
• rlbfgs Riemannian limited memory BFGS solver for smooth objective functions.
This function is called by:

## SOURCE CODE

```0001 function X = thomson_problem(n, d)
0002 % Simple attempt at computing n well distributed points on a sphere in R^d.
0003 %
0004 % This is an example of how Manopt can approximate the gradient and even
0005 % the Hessian of a cost function based on finite differences, even if only
0006 % the cost function is specified without its derivatives.
0007 %
0008 % This functionality is provided only as a help for prototyping, and should
0009 % not be used to compare algorithms in terms of computation time or
0010 % accuracy, because the underlying gradient approximation scheme is slow.
0011 %
0014
0015 % This file is part of Manopt: www.manopt.org.
0016 % Original author: Nicolas Boumal, Nov. 1, 2016
0017 % Contributors:
0018 % Change log:
0019 %
0020 %   Xiaowen Jiang Aug. 20, 2021
0022
0023 if ~exist('n', 'var') || isempty(n)
0024     n = 50;
0025 end
0026 if ~exist('d', 'var') || isempty(d)
0027     d = 3;
0028 end
0029
0030 % Define the Thomson problem with 1/r^2 potential. That is: find n points
0031 % x_i on a sphere in R^d such that the sum over all pairs (i, j) of the
0032 % potentials 1/||x_i - x_j||^2 is minimized. Since the points are on a
0033 % sphere, each potential is .5/(1-x_i'*x_j).
0034 problem.M = obliquefactory(d, n);
0035 problem.cost = @(X) sum(sum(triu(1./(1-X'*X), 1))) / n^2;
0036
0037 % From Matlab 2021a, computating the egrad and the ehess via automatic
0038 % differentiation is available. Notice that the function triu is not
0039 % supported for AD so far.Replace it with ctriu described in the file
0040 % manoptADhelp.m. Also, in this particular case, 1./(1-X'*X) may contain
0041 % NaN on the diagonal which can cause numerical issues when computing the
0042 % egrad via AD although the cost is not a function of the diagonal
0043 % elements. To avoid this problem, first take the upper triangular part by
0044 % calling ctriu before dot division.
0045 % problem.cost = @(X) sum(sum(ctriu(1./ctriu((1-X'*X), 1),1))) / n^2;
0047
0048 % Attempt to minimize the cost. Since the gradient is not provided, Manopt
0049 % approximates it with finite differences. This is /slow/, since for each
0050 % gradient approximation, problem.M.dim()+1 calls to the cost function are
0051 % necessary, on top of generating an orthonormal basis of the tangent space
0052 % at each iterate.
0053 %
0054 % Note that it is difficult to reach high accuracy critical points with an
0055 % approximate gradient, hence it may be required to set a less ambitious
0056 % value for the gradient norm tolerance.
0058
0059 % Pick a solver. Both work fairly well on this problem.
0060 % X = conjugategradient(problem, [], opts);
0061 X = rlbfgs(problem, [], opts);
0062
0063 % Plot the points on a translucid sphere
0064 if nargout == 0 && d == 3
0065     [x, y, z] = sphere(50);
0066     surf(x, y, z, 'FaceAlpha', .5);
0067     hold all;
0068     plot3(X(1, :), X(2, :), X(3, :), '.', 'MarkerSize', 20);
0069     axis equal;
0070     box off;
0071     axis off;
0072 end
0073
0074 % For much better performance, after an early prototyping phase, the
0075 % gradient of the cost function should be specified, typically in