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. See also the derivative free solvers for an alternative: pso and neldermead.
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 % 0012 % See also the derivative free solvers for an alternative: 0013 % pso and neldermead. 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 0021 % Added AD to compute the egrad and the ehess 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; 0046 % problem = manoptAD(problem); 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. 0057 opts.tolgradnorm = 1e-4; 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 0076 % problem.grad or problem.egrad. See the online document at 0077 % 0078 % http://www.manopt.org 0079 % 0080 % for more information. 0081 0082 end