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:
• obliquefactory Returns a manifold struct to optimize over matrices w/ unit-norm columns.
• hessianspectrum Returns the eigenvalues of the (preconditioned) Hessian at x.
This function is called by:

## 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)
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;
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)
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)
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;
0138
0139     % For debugging, it's always nice to check the gradient a few times.
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.
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 Mon 10-Sep-2018 11:48:06 by m2html © 2005