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
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 % Jan. 8, 2021 (NB) 0057 % Added tangent2ambient/tangent2ambient_is_identity pair. 0058 0059 0060 if ~exist('m', 'var') || isempty(m) 0061 m = 1; 0062 end 0063 if ~exist('gpuflag', 'var') || isempty(gpuflag) 0064 gpuflag = false; 0065 end 0066 0067 % If gpuflag is active, new arrays (e.g., via rand, randn, zeros, ones) 0068 % are created directly on the GPU; otherwise, they are created in the 0069 % usual way (in double precision). 0070 if gpuflag 0071 array_type = 'gpuArray'; 0072 else 0073 array_type = 'double'; 0074 end 0075 0076 0077 if m == 1 0078 M.name = @() sprintf('Sphere S^%d', n-1); 0079 else 0080 M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m); 0081 end 0082 0083 M.dim = @() n*m-1; 0084 0085 M.inner = @(x, d1, d2) d1(:)'*d2(:); 0086 0087 M.norm = @(x, d) norm(d, 'fro'); 0088 0089 M.dist = @dist; 0090 function d = dist(x, y) 0091 0092 % The following code is mathematically equivalent to the 0093 % computation d = acos(x(:)'*y(:)) but is much more accurate when 0094 % x and y are close. 0095 0096 chordal_distance = norm(x - y, 'fro'); 0097 d = real(2*asin(.5*chordal_distance)); 0098 0099 % Note: for x and y almost antipodal, the accuracy is good but not 0100 % as good as possible. One way to improve it is by using the 0101 % following branching: 0102 % % if chordal_distance > 1.9 0103 % % d = pi - dist(x, -y); 0104 % % end 0105 % It is rarely necessary to compute the distance between 0106 % almost-antipodal points with full accuracy in Manopt, hence we 0107 % favor a simpler code. 0108 0109 end 0110 0111 M.typicaldist = @() pi; 0112 0113 M.proj = @(x, d) d - x*(x(:)'*d(:)); 0114 0115 M.tangent = M.proj; 0116 0117 M.tangent2ambient_is_identity = true; 0118 M.tangent2ambient = @(X, U) U; 0119 0120 % For Riemannian submanifolds, converting a Euclidean gradient into a 0121 % Riemannian gradient amounts to an orthogonal projection. 0122 M.egrad2rgrad = M.proj; 0123 0124 M.ehess2rhess = @ehess2rhess; 0125 function rhess = ehess2rhess(x, egrad, ehess, u) 0126 rhess = M.proj(x, ehess - (x(:)'*egrad(:))*u); 0127 end 0128 0129 M.exp = @exponential; 0130 0131 M.retr = @retraction; 0132 M.invretr = @inverse_retraction; 0133 0134 M.log = @logarithm; 0135 function v = logarithm(x1, x2) 0136 v = M.proj(x1, x2 - x1); 0137 di = M.dist(x1, x2); 0138 % If the two points are "far apart", correct the norm. 0139 if di > 1e-6 0140 nv = norm(v, 'fro'); 0141 v = v * (di / nv); 0142 end 0143 end 0144 0145 M.hash = @(x) ['z' hashmd5(x(:))]; 0146 0147 M.rand = @() random(n, m, array_type); 0148 0149 M.randvec = @(x) randomvec(n, m, x, array_type); 0150 0151 M.zerovec = @(x) zeros(n, m, array_type); 0152 0153 M.lincomb = @matrixlincomb; 0154 0155 M.transp = @(x1, x2, d) M.proj(x2, d); 0156 0157 % Isometric vector transport of d from the tangent space at x1 to x2. 0158 % This is actually a parallel vector transport, see Ch. 5 in 0159 % http://epubs.siam.org/doi/pdf/10.1137/16M1069298 0160 % "A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization 0161 % on Manifolds", by Hosseini and Uschmajew, SIOPT 2017 0162 M.isotransp = @(x1, x2, d) isometricTransp(x1, x2, d); 0163 function Td = isometricTransp(x1, x2, d) 0164 v = logarithm(x1, x2); 0165 dist_x1x2 = norm(v, 'fro'); 0166 if dist_x1x2 > 0 0167 u = v / dist_x1x2; 0168 utd = u(:)'*d(:); 0169 Td = d + (cos(dist_x1x2)-1)*utd*u ... 0170 - sin(dist_x1x2) *utd*x1; 0171 else 0172 % x1 == x2, so the transport is identity 0173 Td = d; 0174 end 0175 end 0176 0177 M.pairmean = @pairmean; 0178 function y = pairmean(x1, x2) 0179 y = x1+x2; 0180 y = y / norm(y, 'fro'); 0181 end 0182 0183 M.vec = @(x, u_mat) u_mat(:); 0184 M.mat = @(x, u_vec) reshape(u_vec, [n, m]); 0185 M.vecmatareisometries = @() true; 0186 0187 0188 % Automatically convert a number of tools to support GPU. 0189 if gpuflag 0190 M = factorygpuhelper(M); 0191 end 0192 0193 0194 end 0195 0196 % Exponential on the sphere 0197 function y = exponential(x, d, t) 0198 0199 if nargin == 2 0200 % t = 1 0201 td = d; 0202 else 0203 td = t*d; 0204 end 0205 0206 nrm_td = norm(td, 'fro'); 0207 y = x*cos(nrm_td) + td*sinxoverx(nrm_td); 0208 0209 end 0210 0211 % Retraction on the sphere 0212 function y = retraction(x, d, t) 0213 0214 if nargin == 2 0215 % t = 1; 0216 td = d; 0217 else 0218 td = t*d; 0219 end 0220 0221 y = x + td; 0222 y = y / norm(y, 'fro'); 0223 0224 end 0225 0226 % Given x and y two points on the manifold, if there exists a tangent 0227 % vector d at x such that Retr_x(d) = y, this function returns d. 0228 function d = inverse_retraction(x, y) 0229 0230 % Since 0231 % x + d = y*||x + d|| 0232 % and x'd = 0, multiply the above by x' on the left: 0233 % 1 + 0 = x'y * ||x + d|| 0234 % Then solve for d: 0235 0236 d = y/(x(:)'*y(:)) - x; 0237 0238 end 0239 0240 % Uniform random sampling on the sphere. 0241 function x = random(n, m, array_type) 0242 0243 x = randn(n, m, array_type); 0244 x = x / norm(x, 'fro'); 0245 0246 end 0247 0248 % Random normalized tangent vector at x. 0249 function d = randomvec(n, m, x, array_type) 0250 0251 d = randn(n, m, array_type); 0252 d = d - x*(x(:)'*d(:)); 0253 d = d / norm(d, 'fro'); 0254 0255 end