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, gpuflag)

DESCRIPTION ^

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

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

 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.

 Set gpuflag = true to have points, tangent vectors and ambient vectors
 stored on the GPU. If so, computations can be done on the GPU directly.
 
 See also: obliquefactory spherecomplexfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = spherefactory(n, m, gpuflag)
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 % function M = spherefactory(n, m, gpuflag)
0007 %
0008 % Manifold of n-by-m real matrices of unit Frobenius norm.
0009 % By default, m = 1, which corresponds to the unit sphere in R^n. The
0010 % metric is such that the sphere is a Riemannian submanifold of the space
0011 % of nxm matrices with the usual trace inner product, i.e., the usual
0012 % metric.
0013 %
0014 % Set gpuflag = true to have points, tangent vectors and ambient vectors
0015 % stored on the GPU. If so, computations can be done on the GPU directly.
0016 %
0017 % See also: obliquefactory spherecomplexfactory
0018 
0019 % This file is part of Manopt: www.manopt.org.
0020 % Original author: Nicolas Boumal, Dec. 30, 2012.
0021 % Contributors:
0022 % Change log:
0023 %
0024 %   Oct. 8, 2016 (NB)
0025 %       Code for exponential was simplified to only treat the zero vector
0026 %       as a particular case.
0027 %
0028 %   Oct. 22, 2016 (NB)
0029 %       Distance function dist now significantly more accurate for points
0030 %       within 1e-7 and less from each other.
0031 %
0032 %   July 20, 2017 (NB)
0033 %       Following conversations with Bruno Iannazzo and P.-A. Absil,
0034 %       the distance function is now even more accurate.
0035 %
0036 %   Sep. 7, 2017 (NB)
0037 %       New isometric vector transport available in M.isotransp,
0038 %       contributed by Changshuo Liu.
0039 %
0040 %   April 17, 2018 (NB)
0041 %       ehess2rhess: Used to compute projection of ehess, then subtract a
0042 %       multiple of u (which is assumed tangent.) Now, similarly to what
0043 %       happens in stiefelfactory, we first subtract the multiple of u from
0044 %       ehess, then we project. Mathematically, these operations are the
0045 %       same. Numerically, the former version used to be better because tCG
0046 %       in trustregions had some drift near fine convergence. Now that the
0047 %       drift in tCG has been fixed, it is reasonable to apply the
0048 %       projection last, to ensure best tangency of the output.
0049 %
0050 %   July 18, 2018 (NB)
0051 %       Added the inverse retraction (M.invretr) for the sphere.
0052 %
0053 %   Aug. 3, 2018 (NB)
0054 %       Added GPU support: just set gpuflag = true.
0055     
0056     
0057     if ~exist('m', 'var') || isempty(m)
0058         m = 1;
0059     end
0060     if ~exist('gpuflag', 'var') || isempty(gpuflag)
0061         gpuflag = false;
0062     end
0063     
0064     % If gpuflag is active, new arrays (e.g., via rand, randn, zeros, ones)
0065     % are created directly on the GPU; otherwise, they are created in the
0066     % usual way (in double precision).
0067     if gpuflag
0068         array_type = 'gpuArray';
0069     else
0070         array_type = 'double';
0071     end
0072         
0073 
0074     if m == 1
0075         M.name = @() sprintf('Sphere S^%d', n-1);
0076     else
0077         M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m);
0078     end
0079     
0080     M.dim = @() n*m-1;
0081     
0082     M.inner = @(x, d1, d2) d1(:)'*d2(:);
0083     
0084     M.norm = @(x, d) norm(d, 'fro');
0085     
0086     M.dist = @dist;
0087     function d = dist(x, y)
0088         
0089         % The following code is mathematically equivalent to the
0090         % computation d = acos(x(:)'*y(:)) but is much more accurate when
0091         % x and y are close.
0092         
0093         chordal_distance = norm(x - y, 'fro');
0094         d = real(2*asin(.5*chordal_distance));
0095         
0096         % Note: for x and y almost antipodal, the accuracy is good but not
0097         % as good as possible. One way to improve it is by using the
0098         % following branching:
0099         % % if chordal_distance > 1.9
0100         % %     d = pi - dist(x, -y);
0101         % % end
0102         % It is rarely necessary to compute the distance between
0103         % almost-antipodal points with full accuracy in Manopt, hence we
0104         % favor a simpler code.
0105         
0106     end
0107     
0108     M.typicaldist = @() pi;
0109     
0110     M.proj = @(x, d) d - x*(x(:)'*d(:));
0111     
0112     M.tangent = M.proj;
0113     
0114     % For Riemannian submanifolds, converting a Euclidean gradient into a
0115     % Riemannian gradient amounts to an orthogonal projection.
0116     M.egrad2rgrad = M.proj;
0117     
0118     M.ehess2rhess = @ehess2rhess;
0119     function rhess = ehess2rhess(x, egrad, ehess, u)
0120         rhess = M.proj(x, ehess - (x(:)'*egrad(:))*u);
0121     end
0122     
0123     M.exp = @exponential;
0124     
0125     M.retr = @retraction;
0126     M.invretr = @inverse_retraction;
0127 
0128     M.log = @logarithm;
0129     function v = logarithm(x1, x2)
0130         v = M.proj(x1, x2 - x1);
0131         di = M.dist(x1, x2);
0132         % If the two points are "far apart", correct the norm.
0133         if di > 1e-6
0134             nv = norm(v, 'fro');
0135             v = v * (di / nv);
0136         end
0137     end
0138     
0139     M.hash = @(x) ['z' hashmd5(x(:))];
0140     
0141     M.rand = @() random(n, m, array_type);
0142     
0143     M.randvec = @(x) randomvec(n, m, x, array_type);
0144     
0145     M.zerovec = @(x) zeros(n, m, array_type);
0146     
0147     M.lincomb = @matrixlincomb;
0148     
0149     M.transp = @(x1, x2, d) M.proj(x2, d);
0150     
0151     % Isometric vector transport of d from the tangent space at x1 to x2.
0152     % This is actually a parallel vector transport, see 5 in
0153     % http://epubs.siam.org/doi/pdf/10.1137/16M1069298
0154     % "A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization
0155     %  on Manifolds", by Hosseini and Uschmajew, SIOPT 2017
0156     M.isotransp = @(x1, x2, d) isometricTransp(x1, x2, d);
0157     function Td = isometricTransp(x1, x2, d)
0158         v = logarithm(x1, x2);
0159         dist_x1x2 = norm(v, 'fro');
0160         if dist_x1x2 > 0
0161             u = v / dist_x1x2;
0162             utd = u(:)'*d(:);
0163             Td = d + (cos(dist_x1x2)-1)*utd*u ...
0164                     -  sin(dist_x1x2)   *utd*x1;
0165         else
0166             % x1 == x2, so the transport is identity
0167             Td = d;
0168         end
0169     end
0170     
0171     M.pairmean = @pairmean;
0172     function y = pairmean(x1, x2)
0173         y = x1+x2;
0174         y = y / norm(y, 'fro');
0175     end
0176 
0177     M.vec = @(x, u_mat) u_mat(:);
0178     M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
0179     M.vecmatareisometries = @() true;
0180     
0181     
0182     % Automatically convert a number of tools to support GPU.
0183     if gpuflag
0184         M = factorygpuhelper(M);
0185     end
0186     
0187 
0188 end
0189 
0190 % Exponential on the sphere
0191 function y = exponential(x, d, t)
0192 
0193     if nargin == 2
0194         % t = 1
0195         td = d;
0196     else
0197         td = t*d;
0198     end
0199     
0200     nrm_td = norm(td, 'fro');
0201     
0202     % Former versions of Manopt avoided the computation of sin(a)/a for
0203     % small a, but further investigations suggest this computation is
0204     % well-behaved numerically.
0205     if nrm_td > 0
0206         y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
0207     else
0208         y = x;
0209     end
0210 
0211 end
0212 
0213 % Retraction on the sphere
0214 function y = retraction(x, d, t)
0215 
0216     if nargin == 2
0217         % t = 1;
0218         td = d;
0219     else
0220         td = t*d;
0221     end
0222     
0223     y = x + td;
0224     y = y / norm(y, 'fro');
0225 
0226 end
0227 
0228 % Given x and y two points on the manifold, if there exists a tangent
0229 % vector d at x such that Retr_x(d) = y, this function returns d.
0230 function d = inverse_retraction(x, y)
0231 
0232     % Since
0233     %   x + d = y*||x + d||
0234     % and x'd = 0, multiply the above by x' on the left:
0235     %   1 + 0 = x'y * ||x + d||
0236     % Then solve for d:
0237     
0238     d = y/(x(:)'*y(:)) - x;
0239 
0240 end
0241 
0242 % Uniform random sampling on the sphere.
0243 function x = random(n, m, array_type)
0244 
0245     x = randn(n, m, array_type);
0246     x = x / norm(x, 'fro');
0247 
0248 end
0249 
0250 % Random normalized tangent vector at x.
0251 function d = randomvec(n, m, x, array_type)
0252 
0253     d = randn(n, m, array_type);
0254     d = d - x*(x(:)'*d(:));
0255     d = d / norm(d, 'fro');
0256 
0257 end

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005