Home > manopt > manifolds > sphere > spherefactory.m

spherefactory

PURPOSE ^

Returns a manifold struct to optimize over unit-norm vectors or matrices.

SYNOPSIS ^

function M = spherefactory(n, m)

DESCRIPTION ^

 Returns a manifold struct to optimize over unit-norm vectors or matrices.

 function M = spherefactory(n)
 function M = spherefactory(n, m)

 Manifold of n-by-m real matrices of unit Frobenius norm.
 By default, m = 1, which corresponds to the unit sphere in R^n. The
 metric is such that the sphere is a Riemannian submanifold of the space
 of nxm matrices with the usual trace inner product, i.e., the usual
 metric.
 
 See also: obliquefactory spherecomplexfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = spherefactory(n, m)
0002 % Returns a manifold struct to optimize over unit-norm vectors or matrices.
0003 %
0004 % function M = spherefactory(n)
0005 % function M = spherefactory(n, m)
0006 %
0007 % Manifold of n-by-m real matrices of unit Frobenius norm.
0008 % By default, m = 1, which corresponds to the unit sphere in R^n. The
0009 % metric is such that the sphere is a Riemannian submanifold of the space
0010 % of nxm matrices with the usual trace inner product, i.e., the usual
0011 % metric.
0012 %
0013 % See also: obliquefactory spherecomplexfactory
0014 
0015 % This file is part of Manopt: www.manopt.org.
0016 % Original author: Nicolas Boumal, Dec. 30, 2012.
0017 % Contributors:
0018 % Change log:
0019 %
0020 %   Oct. 8, 2016 (NB)
0021 %       Code for exponential was simplified to only treat the zero vector
0022 %       as a particular case.
0023 %
0024 %   Oct. 22, 2016 (NB)
0025 %       Distance function dist now significantly more accurate for points
0026 %       within 1e-7 and less from each other.
0027 %
0028 %   July 20, 2017 (NB)
0029 %       Following conversations with Bruno Iannazzo and P.-A. Absil,
0030 %       the distance function is now even more accurate.
0031 %
0032 %   Sep. 7, 2017 (NB)
0033 %       New isometric vector transport available in M.isotransp,
0034 %       contributed by Changshuo Liu.
0035 
0036     
0037     if ~exist('m', 'var')
0038         m = 1;
0039     end
0040 
0041     if m == 1
0042         M.name = @() sprintf('Sphere S^%d', n-1);
0043     else
0044         M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m);
0045     end
0046     
0047     M.dim = @() n*m-1;
0048     
0049     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0050     
0051     M.norm = @(x, d) norm(d, 'fro');
0052     
0053     M.dist = @dist;
0054     function d = dist(x, y)
0055         
0056         % The following code is mathematically equivalent to the
0057         % computation d = acos(x(:)'*y(:)) but is much more accurate when
0058         % x and y are close.
0059         
0060         chordal_distance = norm(x - y, 'fro');
0061         d = real(2*asin(.5*chordal_distance));
0062         
0063         % Note: for x and y almost antipodal, the accuracy is good but not
0064         % as good as possible. One way to improve it is by using the
0065         % following branching:
0066         % % if chordal_distance > 1.9
0067         % %     d = pi - dist(x, -y);
0068         % % end
0069         % It is rarely necessary to compute distance between
0070         % almost-antipodal points with full accuracy in Manopt, hence we
0071         % favor a simpler code.
0072         
0073     end
0074     
0075     M.typicaldist = @() pi;
0076     
0077     M.proj = @(x, d) d - x*(x(:).'*d(:));
0078     
0079     M.tangent = M.proj;
0080     
0081     % For Riemannian submanifolds, converting a Euclidean gradient into a
0082     % Riemannian gradient amounts to an orthogonal projection.
0083     M.egrad2rgrad = M.proj;
0084     
0085     M.ehess2rhess = @ehess2rhess;
0086     function rhess = ehess2rhess(x, egrad, ehess, u)
0087         rhess = M.proj(x, ehess) - (x(:)'*egrad(:))*u;
0088     end
0089     
0090     M.exp = @exponential;
0091     
0092     M.retr = @retraction;
0093 
0094     M.log = @logarithm;
0095     function v = logarithm(x1, x2)
0096         v = M.proj(x1, x2 - x1);
0097         di = M.dist(x1, x2);
0098         % If the two points are "far apart", correct the norm.
0099         if di > 1e-6
0100             nv = norm(v, 'fro');
0101             v = v * (di / nv);
0102         end
0103     end
0104     
0105     M.hash = @(x) ['z' hashmd5(x(:))];
0106     
0107     M.rand = @() random(n, m);
0108     
0109     M.randvec = @(x) randomvec(n, m, x);
0110     
0111     M.lincomb = @matrixlincomb;
0112     
0113     M.zerovec = @(x) zeros(n, m);
0114     
0115     M.transp = @(x1, x2, d) M.proj(x2, d);
0116     
0117     % Isometric vector transport of d from the tangent space at x1 to x2.
0118     % This is actually a parallel vector transport, see 5 in
0119     % http://epubs.siam.org/doi/pdf/10.1137/16M1069298
0120     % "A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization
0121     %  on Manifolds", by Hosseini and Uschmajew, SIOPT 2017
0122     M.isotransp = @(x1, x2, d) isometricTransp(x1, x2, d);
0123     function Td = isometricTransp(x1, x2, d)
0124         v = logarithm(x1, x2);
0125         dist_x1x2 = norm(v, 'fro');
0126         if dist_x1x2 > 0
0127             u = v / dist_x1x2;
0128             utd = u(:)'*d(:);
0129             Td = d + (cos(dist_x1x2)-1)*utd*u ...
0130                     -  sin(dist_x1x2)   *utd*x1;
0131         else
0132             % x1 == x2, so the transport is identity
0133             Td = d;
0134         end
0135     end
0136     
0137     M.pairmean = @pairmean;
0138     function y = pairmean(x1, x2)
0139         y = x1+x2;
0140         y = y / norm(y, 'fro');
0141     end
0142 
0143     M.vec = @(x, u_mat) u_mat(:);
0144     M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
0145     M.vecmatareisometries = @() true;
0146 
0147 end
0148 
0149 % Exponential on the sphere
0150 function y = exponential(x, d, t)
0151 
0152     if nargin == 2
0153         % t = 1
0154         td = d;
0155     else
0156         td = t*d;
0157     end
0158     
0159     nrm_td = norm(td, 'fro');
0160     
0161     % Former versions of Manopt avoided the computation of sin(a)/a for
0162     % small a, but further investigations suggest this computation is
0163     % well-behaved numerically.
0164     if nrm_td > 0
0165         y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
0166     else
0167         y = x;
0168     end
0169 
0170 end
0171 
0172 % Retraction on the sphere
0173 function y = retraction(x, d, t)
0174 
0175     if nargin == 2
0176         % t = 1;
0177         td = d;
0178     else
0179         td = t*d;
0180     end
0181     
0182     y = x + td;
0183     y = y / norm(y, 'fro');
0184 
0185 end
0186 
0187 % Uniform random sampling on the sphere.
0188 function x = random(n, m)
0189 
0190     x = randn(n, m);
0191     x = x / norm(x, 'fro');
0192 
0193 end
0194 
0195 % Random normalized tangent vector at x.
0196 function d = randomvec(n, m, x)
0197 
0198     d = randn(n, m);
0199     d = d - x*(x(:).'*d(:));
0200     d = d / norm(d, 'fro');
0201 
0202 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005