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     
0029     if ~exist('m', 'var')
0030         m = 1;
0031     end
0032 
0033     if m == 1
0034         M.name = @() sprintf('Sphere S^%d', n-1);
0035     else
0036         M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m);
0037     end
0038     
0039     M.dim = @() n*m-1;
0040     
0041     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0042     
0043     M.norm = @(x, d) norm(d, 'fro');
0044     
0045     M.dist = @dist;
0046     function d = dist(x, y)
0047         
0048         % This computation cannot be accurate below an output of 2e-8.
0049         % The reason is: if two unit-norm vectors x and y are very close to
0050         % one another, their inner product is about 1. The machine
0051         % precision at 1 is eps(1) = 2e-16. The correct value for
0052         % acos(1-eps(1)) is about 2e-8. This can be checked with the
0053         % syms toolbox: syms x; f = acos(1-x); vpa(subs(f, x, eps(1)), 32)
0054         % Thus, if x and y are actually closer to each other than 2e-8,
0055         % their inner product will be even closer to 1, but that cannot be
0056         % represented in IEEE arithmetic. Thus, their inner product will be
0057         % rounded to either 1 (giving 0 distance) or to 1-eps(1), which
0058         % gives a distance of 2e-8, or to something even further from 1. No
0059         % distance between 0 and 2e-8 can thus be computed this way.
0060         d = real(acos(x(:).'*y(:)));
0061         
0062         % Hence, if the distance proves dangerously small so that it is
0063         % possible that we suffered from round-off, we compute the distance
0064         % in the embedding space instead. At this scale, this is quite
0065         % accurate.
0066         if d < 1e-7
0067             d = norm(x-y, 'fro');
0068         end
0069     end
0070     
0071     M.typicaldist = @() pi;
0072     
0073     M.proj = @(x, d) d - x*(x(:).'*d(:));
0074     
0075     M.tangent = M.proj;
0076     
0077     % For Riemannian submanifolds, converting a Euclidean gradient into a
0078     % Riemannian gradient amounts to an orthogonal projection.
0079     M.egrad2rgrad = M.proj;
0080     
0081     M.ehess2rhess = @ehess2rhess;
0082     function rhess = ehess2rhess(x, egrad, ehess, u)
0083         rhess = M.proj(x, ehess) - (x(:)'*egrad(:))*u;
0084     end
0085     
0086     M.exp = @exponential;
0087     
0088     M.retr = @retraction;
0089 
0090     M.log = @logarithm;
0091     function v = logarithm(x1, x2)
0092         v = M.proj(x1, x2 - x1);
0093         di = M.dist(x1, x2);
0094         % If the two points are "far apart", correct the norm.
0095         if di > 1e-6
0096             nv = norm(v, 'fro');
0097             v = v * (di / nv);
0098         end
0099     end
0100     
0101     M.hash = @(x) ['z' hashmd5(x(:))];
0102     
0103     M.rand = @() random(n, m);
0104     
0105     M.randvec = @(x) randomvec(n, m, x);
0106     
0107     M.lincomb = @matrixlincomb;
0108     
0109     M.zerovec = @(x) zeros(n, m);
0110     
0111     M.transp = @(x1, x2, d) M.proj(x2, d);
0112     
0113     M.pairmean = @pairmean;
0114     function y = pairmean(x1, x2)
0115         y = x1+x2;
0116         y = y / norm(y, 'fro');
0117     end
0118 
0119     M.vec = @(x, u_mat) u_mat(:);
0120     M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
0121     M.vecmatareisometries = @() true;
0122 
0123 end
0124 
0125 % Exponential on the sphere
0126 function y = exponential(x, d, t)
0127 
0128     if nargin == 2
0129         % t = 1
0130         td = d;
0131     else
0132         td = t*d;
0133     end
0134     
0135     nrm_td = norm(td, 'fro');
0136     
0137     % Former versions of Manopt avoided the computation of sin(a)/a for
0138     % small a, but further investigations suggest this computation is
0139     % well-behaved numerically.
0140     if nrm_td > 0
0141         y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
0142     else
0143         y = x;
0144     end
0145 
0146 end
0147 
0148 % Retraction on the sphere
0149 function y = retraction(x, d, t)
0150 
0151     if nargin == 2
0152         % t = 1;
0153         td = d;
0154     else
0155         td = t*d;
0156     end
0157     
0158     y = x + td;
0159     y = y / norm(y, 'fro');
0160 
0161 end
0162 
0163 % Uniform random sampling on the sphere.
0164 function x = random(n, m)
0165 
0166     x = randn(n, m);
0167     x = x / norm(x, 'fro');
0168 
0169 end
0170 
0171 % Random normalized tangent vector at x.
0172 function d = randomvec(n, m, x)
0173 
0174     d = randn(n, m);
0175     d = d - x*(x(:).'*d(:));
0176     d = d / norm(d, 'fro');
0177 
0178 end

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