Home > manopt > manifolds > sphere > spheresymmetricfactory.m

spheresymmetricfactory

PURPOSE ^

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

SYNOPSIS ^

function M = spheresymmetricfactory(n)

DESCRIPTION ^

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

 function M = spheresymmetricfactory(n)

 Manifold of n-by-n real symmetric matrices of unit Frobenius norm.
 The metric is such that the sphere is a Riemannian submanifold of the
 space of nxn symmetric matrices with the usual trace inner product, i.e.,
 the usual metric <A, B> = trace(A'*B).
 
 See also: spherefactory obliquefactory spherecomplexfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = spheresymmetricfactory(n)
0002 % Returns a manifold struct to optimize over unit-norm symmetric matrices.
0003 %
0004 % function M = spheresymmetricfactory(n)
0005 %
0006 % Manifold of n-by-n real symmetric matrices of unit Frobenius norm.
0007 % The metric is such that the sphere is a Riemannian submanifold of the
0008 % space of nxn symmetric matrices with the usual trace inner product, i.e.,
0009 % the usual metric <A, B> = trace(A'*B).
0010 %
0011 % See also: spherefactory obliquefactory spherecomplexfactory
0012 
0013 % This file is part of Manopt: www.manopt.org.
0014 % Original author: Nicolas Boumal, April 17, 2015.
0015 % Contributors:
0016 % Change log:
0017 %
0018 %   Oct. 8, 2016 (NB)
0019 %       Code for exponential was simplified to only treat the zero vector
0020 %       as a particular case.
0021 %
0022 %   Oct. 22, 2016 (NB)
0023 %       Distance function dist now significantly more accurate for points
0024 %       within 1e-7 and less from each other.
0025 
0026 
0027     M.name = @() sprintf('Sphere of symmetric matrices of size %d', n);
0028     
0029     M.dim = @() n*(n+1)/2 - 1;
0030     
0031     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0032     
0033     M.norm = @(x, d) norm(d, 'fro');
0034     
0035     M.dist = @dist;
0036     function d = dist(x, y)
0037         d = real(acos(x(:).'*y(:)));
0038         % The above formula is numerically inaccurate if x and y are too
0039         % close together. In that case, norm is a much better proxy.
0040         if d < 1e-6
0041             d = norm(x-y, 'fro');
0042         end
0043     end
0044     
0045     M.typicaldist = @() pi;
0046     
0047     M.proj = @proj;
0048     function xdot = proj(x, d)
0049         d = (d+d.')/2;
0050         xdot = d - x*(x(:).'*d(:));
0051     end
0052     
0053     M.tangent = @proj;
0054     
0055     % For Riemannian submanifolds, converting a Euclidean gradient into a
0056     % Riemannian gradient amounts to an orthogonal projection.
0057     M.egrad2rgrad = @proj;
0058     
0059     M.ehess2rhess = @ehess2rhess;
0060     function rhess = ehess2rhess(x, egrad, ehess, u)
0061         % these are not explicitly required, given the use.
0062         % egrad = (egrad + egrad.')/2;
0063         % ehess = (ehess + ehess.')/2;
0064         rhess = proj(x, ehess) - (x(:)'*egrad(:))*u;
0065     end
0066     
0067     M.exp = @exponential;
0068     
0069     M.retr = @retraction;
0070 
0071     M.log = @logarithm;
0072     function v = logarithm(x1, x2)
0073         v = proj(x1, x2 - x1);
0074         di = M.dist(x1, x2);
0075         % If the two points are "far apart", correct the norm.
0076         if di > 1e-6
0077             nv = norm(v, 'fro');
0078             v = v * (di / nv);
0079         end
0080     end
0081     
0082     M.hash = @(x) ['z' hashmd5(x(:))];
0083     
0084     M.rand = @() random(n);
0085     
0086     M.randvec = @(x) randomvec(n, x);
0087     
0088     M.lincomb = @matrixlincomb;
0089     
0090     M.zerovec = @(x) zeros(n);
0091     
0092     M.transp = @(x1, x2, d) proj(x2, d);
0093     
0094     M.pairmean = @pairmean;
0095     function y = pairmean(x1, x2)
0096         y = x1+x2;
0097         y = y / norm(y, 'fro');
0098     end
0099 
0100     % TODO : check isometry and fix.
0101     M.vec = @(x, u_mat) u_mat(:);
0102     M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
0103     M.vecmatareisometries = @() false;
0104 
0105 end
0106 
0107 % Exponential on the sphere
0108 function y = exponential(x, d, t)
0109 
0110     if nargin == 2
0111         % t = 1;
0112         td = d;
0113     else
0114         td = t*d;
0115     end
0116     
0117     nrm_td = norm(td, 'fro');
0118     
0119     if nrm_td > 0
0120         y = x*cos(nrm_td) + td*(sin(nrm_td)/nrm_td);
0121     else
0122         y = x;
0123     end
0124 
0125 end
0126 
0127 % Retraction on the sphere
0128 function y = retraction(x, d, t)
0129 
0130     if nargin == 2
0131         t = 1;
0132     end
0133     
0134     y = x + t*d;
0135     y = y / norm(y, 'fro');
0136 
0137 end
0138 
0139 % Uniform random sampling on the sphere.
0140 function x = random(n)
0141 
0142     x = randn(n);
0143     x = (x + x.')/2;
0144     x = x/norm(x, 'fro');
0145 
0146 end
0147 
0148 % Random normalized tangent vector at x.
0149 function d = randomvec(n, x)
0150 
0151     d = randn(n);
0152     d = (d + d.')/2;
0153     d = d - x*(x(:).'*d(:));
0154     d = d / norm(d, 'fro');
0155 
0156 end

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