Home > examples > packing_on_the_sphere.m

packing_on_the_sphere

PURPOSE ^

Return a set of points spread out on the sphere.

SYNOPSIS ^

function [X, maxdot] = packing_on_the_sphere(d, n, epsilon, X0)

DESCRIPTION ^

 Return a set of points spread out on the sphere.

 function [X, maxdot] = packing_on_the_sphere(d, n, epsilon, X0)

 Using optimization on the oblique manifold, that is, the product of
 spheres, this function returns a set of n points with unit norm in R^d in
 the form of a matrix X of size nxd, such that the points are spread out
 on the sphere. Ideally, we would minimize the maximum inner product
 between any two points X(i, :) and X(j, :), i~=j, but that is a nonsmooth
 cost function. Instead, we replace the max function by a classical
 log-sum-exp approximation and (attempt to) solve:

 min_{X in OB(d, n)} log( .5*sum_{i~=j} exp( xi'*xj/epsilon ) ),

 with xi = X(:, i) and epsilon is some "diffusion constant". As epsilon
 goes to zero, the cost function is a sharper approximation of the max
 function (under some assumptions), but the cost function becomes stiffer
 and hence harder to optimize.

 The second output, maxdot, is the maximum inner product between any two
 points in the returned X. This number is the one we truly are trying to
 minimize.

 Notice that this cost function is invariant under rotation of X:
 f(X) = f(XQ) for all orthogonal Q in O(d).
 This calls for optimization over the set of symmetric positive
 semidefinite matrices of size n and rank d with unit diagonal, which can
 be thought of as the quotient of the oblique manifold OB(d, n) by O(d):
 See elliptopefactory.

 This is known as the Thomson or, more specifically, the Tammes problem:
 http://en.wikipedia.org/wiki/Tammes_problem
 An interesting page by Neil Sloane collecting best known packings is
 available here http://neilsloane.com/packings/

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [X, maxdot] = packing_on_the_sphere(d, n, epsilon, X0)
0002 % Return a set of points spread out on the sphere.
0003 %
0004 % function [X, maxdot] = packing_on_the_sphere(d, n, epsilon, X0)
0005 %
0006 % Using optimization on the oblique manifold, that is, the product of
0007 % spheres, this function returns a set of n points with unit norm in R^d in
0008 % the form of a matrix X of size nxd, such that the points are spread out
0009 % on the sphere. Ideally, we would minimize the maximum inner product
0010 % between any two points X(i, :) and X(j, :), i~=j, but that is a nonsmooth
0011 % cost function. Instead, we replace the max function by a classical
0012 % log-sum-exp approximation and (attempt to) solve:
0013 %
0014 % min_{X in OB(d, n)} log( .5*sum_{i~=j} exp( xi'*xj/epsilon ) ),
0015 %
0016 % with xi = X(:, i) and epsilon is some "diffusion constant". As epsilon
0017 % goes to zero, the cost function is a sharper approximation of the max
0018 % function (under some assumptions), but the cost function becomes stiffer
0019 % and hence harder to optimize.
0020 %
0021 % The second output, maxdot, is the maximum inner product between any two
0022 % points in the returned X. This number is the one we truly are trying to
0023 % minimize.
0024 %
0025 % Notice that this cost function is invariant under rotation of X:
0026 % f(X) = f(XQ) for all orthogonal Q in O(d).
0027 % This calls for optimization over the set of symmetric positive
0028 % semidefinite matrices of size n and rank d with unit diagonal, which can
0029 % be thought of as the quotient of the oblique manifold OB(d, n) by O(d):
0030 % See elliptopefactory.
0031 %
0032 % This is known as the Thomson or, more specifically, the Tammes problem:
0033 % http://en.wikipedia.org/wiki/Tammes_problem
0034 % An interesting page by Neil Sloane collecting best known packings is
0035 % available here http://neilsloane.com/packings/
0036 
0037 % This file is part of Manopt and is copyrighted. See the license file.
0038 %
0039 % Main author: Nicolas Boumal, July 2, 2013
0040 % Contributors:
0041 %
0042 % Change log:
0043 %   Aug. 14, 2013 (NB) : Code now compatible to experiment with both the
0044 %                        obliquefactory and the elliptopefactory.
0045 %
0046 %   Jan.  7, 2014 (NB) : Added reference to Neil Sloane's page and the
0047 %                        maxdot output.
0048 %
0049 %   June 24, 2014 (NB) : Now shifting exponentials to alleviate numerical
0050 %                        trouble when epsilon is too small.
0051 %
0052     
0053     if ~exist('d', 'var') || isempty(d)
0054         % Dimension of the embedding space: R^d
0055         d = 3;
0056     end
0057     if ~exist('n', 'var') || isempty(n)
0058         % Number n of points to place of the sphere in R^d.
0059         % For example, n=12 yields an icosahedron:
0060         % https://en.wikipedia.org/wiki/Icosahedron
0061         % Notice though that platonic solids are not always optimal.
0062         % Try for example n = 8: you don't get a cube.
0063         n = 24;
0064     end
0065     if ~exist('epsilon', 'var') || isempty(epsilon)
0066         % This value should be as close to 0 as affordable.
0067         % If it is too close to zero, optimization first becomes much
0068         % slower, than simply doesn't work anymore becomes of floating
0069         % point overflow errors (NaN's and Inf's start to appear).
0070         % If it is too large, then log-sum-exp is a poor approximation of
0071         % the max function, and the spread will be less uniform.
0072         % An okay value seems to be 0.01 or 0.001 for example. Note that a
0073         % better strategy than using a small epsilon straightaway is to
0074         % reduce epsilon bit by bit and to warm-start subsequent
0075         % optimization in that way. Trustregions will be more appropriate
0076         % for these fine tunings.
0077         epsilon = 0.0015;
0078     end
0079     
0080     % Pick your manifold (the elliptope factory quotients out the global
0081     % rotation invariance of the problem, which is more natural but
0082     % conceptually a bit more complicated --- for usage with the toolbox it
0083     % is the same though: just uncomment the appropriate line).
0084     manifold = obliquefactory(d, n, true);
0085     % manifold = elliptopefactory(n, d);
0086     
0087     % Generate a random initial guess if none was given.
0088     if ~exist('X0', 'var') || isempty(X0)
0089         X0 = manifold.rand();
0090     end
0091 
0092     % Define the cost function with caching system used: the store
0093     % structure we receive as input is tied to the input point X. Everytime
0094     % this cost function is called at this point X, we will receive the
0095     % same store structure back. We may modify the store structure inside
0096     % the function and return it: the changes will be remembered for next
0097     % time.
0098     function [f, store] = cost(X, store)
0099         if ~isfield(store, 'ready')
0100             XXt = X*X';
0101             % Shift the exponentials by the maximum value to reduce
0102             % numerical trouble due to possible overflows.
0103             s = max(max(triu(XXt, 1)));
0104             expXXt = exp((XXt-s)/epsilon);
0105             % Zero out the diagonal
0106             expXXt(1:(n+1):end) = 0;
0107             u = sum(sum(triu(expXXt, 1)));
0108             store.XXt = XXt;
0109             store.s = s;
0110             store.expXXt = expXXt;
0111             store.u = u;
0112             store.ready = true;
0113         end
0114         u = store.u;
0115         s = store.s;
0116         f = s + epsilon*log(u);
0117     end
0118 
0119     % Define the gradient of the cost. When the gradient is called at a
0120     % point X for which the cost was already called, the store structure we
0121     % receive remember everything that the cost function stored in it, so
0122     % we can reuse previously computed elements.
0123     function [g, store] = grad(X, store)
0124         if ~isfield(store, 'ready')
0125             [~, store] = cost(X, store);
0126         end
0127         % Compute the Euclidean gradient
0128         eg = store.expXXt*X / store.u;
0129         % Convert to the Riemannian gradient (by projection)
0130         g = manifold.egrad2rgrad(X, eg);
0131     end
0132 
0133     % Setup the problem structure with its manifold M and cost+grad
0134     % functions.
0135     problem.M = manifold;
0136     problem.cost = @cost;
0137     problem.grad = @grad;
0138 
0139     % For debugging, it's always nice to check the gradient a few times.
0140     % checkgradient(problem);
0141     % pause;
0142     
0143     % Call a solver on our problem with a few options defined. We did not
0144     % specify the Hessian but it is still okay to call trustregion: Manopt
0145     % will approximate the Hessian with finite differences of the gradient.
0146     opts.tolgradnorm = 1e-8;
0147     opts.maxtime = 1200;
0148     opts.maxiter = 1e5;
0149     % X = trustregions(problem, X0, opts);
0150     X = conjugategradient(problem, X0, opts);
0151     
0152     % Evaluate the maximum inner product between any two points of X.
0153     XXt = X*X';
0154     dots = XXt(find(triu(ones(n), 1))); %#ok<FNDSB>
0155     maxdot = max(dots);
0156     
0157     % Similarly, even though we did not specify the Hessian, we may still
0158     % estimate its spectrum at the solution. It should reflect the
0159     % invariance of the cost function under a global rotatioon of the
0160     % sphere, which is an invariance under the group O(d) of dimension
0161     % d(d-1)/2 : this translates into d(d-1)/2 zero eigenvalues in the
0162     % spectrum of the Hessian.
0163     % The approximate Hessian is not a linear operator, and is it a
0164     % fortiori not symmetric. The result of this computation is thus not
0165     % reliable. It does display the zero eigenvalues as expected though.
0166     if manifold.dim() < 300
0167         evs = real(hessianspectrum(problem, X));
0168         figure;
0169         stem(1:length(evs), sort(evs), '.');
0170         title(['Eigenvalues of the approximate Hessian of the cost ' ...
0171                'function at the solution']);
0172     end
0173     
0174     
0175     % Show how the inner products X(:, i)'*X(:, j) are distributed.
0176     figure;
0177     hist(real(acos(dots)), 20);
0178     title('Histogram of the geodesic distances');
0179     
0180     % This is the quantity we actually want to minimize.
0181     fprintf('Maximum inner product between two points: %g\n', maxdot);
0182     
0183     
0184     % Give some visualization if the dimension allows
0185     if d == 2
0186         % For the circle, the optimal solution consists in spreading the
0187         % points with angles uniformly sampled in (0, 2pi). This
0188         % corresponds to the following value for the max inner product:
0189         fprintf('Optimal value for the max inner product: %g\n', cos(2*pi/n));
0190         figure;
0191         t = linspace(-pi, pi, 201);
0192         plot(cos(t), sin(t), '-', 'LineWidth', 3, 'Color', [152,186,220]/255);
0193         daspect([1 1 1]);
0194         box off;
0195         axis off;
0196         hold on;
0197         plot(X(:, 1), X(:, 2), 'r.', 'MarkerSize', 25);
0198         hold off;
0199     end
0200     if d == 3
0201         figure;
0202         % Plot the sphere
0203         [sphere_x, sphere_y, sphere_z] = sphere(50);
0204         handle = surf(sphere_x, sphere_y, sphere_z);
0205         set(handle, 'FaceColor', [152,186,220]/255);
0206         set(handle, 'FaceAlpha', .5);
0207         set(handle, 'EdgeColor', [152,186,220]/255);
0208         set(handle, 'EdgeAlpha', .5);
0209         daspect([1 1 1]);
0210         box off;
0211         axis off;
0212         hold on;
0213         % Add the chosen points
0214         Y = 1.02*X';
0215         plot3(Y(1, :), Y(2, :), Y(3, :), 'r.', 'MarkerSize', 25);
0216         % And connect the points which are at minimal distance,
0217         % within some tolerance.
0218         min_distance = real(acos(maxdot));
0219         connected = real(acos(XXt)) <= 1.20*min_distance;
0220         [Ic, Jc] = find(triu(connected, 1));
0221         for k = 1 : length(Ic)
0222             i = Ic(k); j = Jc(k);
0223             plot3(Y(1, [i j]), Y(2, [i j]), Y(3, [i j]), 'k-');
0224         end
0225         hold off;
0226     end
0227 
0228 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005