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/

- obliquefactory Returns a manifold struct to optimize over matrices w/ unit-norm columns.
- conjugategradient Conjugate gradient minimization algorithm for Manopt.
- hessianspectrum Returns the eigenvalues of the (preconditioned) Hessian at x.

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