Returns a manifold struct. for the sphere on the tangent space to M at x. N = tangentspherefactory(M, x) N defines a manifold that is the unit sphere on the tangent space to M at x. Points are represented as tangent vectors of unit norm. Tangent vectors are represented as tangent vectors orthogonal to the root point, with respect to the Riemannian metric on the tangent space. This is chiefly useful to solve optimization problems involving unit norm tangent vectors to M at x, which notably comes up when looking for extreme eigenvectors of the Hessian of a cost function on M at x, for example. The Riemannian structure on this sphere is that of a Riemannian submanifold of the (Euclidean) tangent space, equipped with the Riemannian metric of M at that point. See also: hessianextreme
0001 function N = tangentspherefactory(M, x) 0002 % Returns a manifold struct. for the sphere on the tangent space to M at x. 0003 % 0004 % N = tangentspherefactory(M, x) 0005 % 0006 % N defines a manifold that is the unit sphere on the tangent space to M 0007 % at x. Points are represented as tangent vectors of unit norm. Tangent 0008 % vectors are represented as tangent vectors orthogonal to the root point, 0009 % with respect to the Riemannian metric on the tangent space. 0010 % 0011 % This is chiefly useful to solve optimization problems involving unit norm 0012 % tangent vectors to M at x, which notably comes up when looking for 0013 % extreme eigenvectors of the Hessian of a cost function on M at x, for 0014 % example. The Riemannian structure on this sphere is that of a Riemannian 0015 % submanifold of the (Euclidean) tangent space, equipped with the 0016 % Riemannian metric of M at that point. 0017 % 0018 % See also: hessianextreme 0019 0020 % This file is part of Manopt: www.manopt.org. 0021 % Original author: Nicolas Boumal, March 16, 2015. 0022 % Contributors: 0023 % Change log: 0024 % 0025 % Nov 27, 2015 (NB): 0026 % Extra projection added in the retraction, to prevent numerical 0027 % drift. 0028 % 0029 % Jun 23, 2021 (QR): 0030 % Extra projection in retraction changed to M.tangent. 0031 % 0032 % Jun 23, 2021 (NB): 0033 % Several fixes to the logic of this factory that should help for 0034 % more sophisticated manifolds where representations of points, 0035 % tangent vectors and ambient vectors are not straightforward. 0036 0037 % N is the manifold we build. 0038 % y is a point on N, thus also a tangent vector to M at x. 0039 % This is a typical Riemannian submanifold of a Euclidean space, 0040 % hence it is easy to describe in terms of the tools available for M. 0041 N = struct(); 0042 N.name = @() sprintf('Sphere in a tangent space of [%s]', M.name()); 0043 0044 % u, u1 and u2 are tangent vectors to N at y. 0045 % The tangent space to N at y is a subspace of the tangent space to M 0046 % at x, thus u, u1 and u2 are also tangent vectors to M at x. 0047 0048 N.dim = @() M.dim() - 1; 0049 N.inner = @(y, u1, u2) M.inner(x, u1, u2); 0050 N.norm = @(y, u) M.norm(x, u); 0051 N.proj = @(y, v) M.lincomb(x, 1, v, -M.inner(x, v, y), y); 0052 N.typicaldist = @() 1; 0053 N.tangent = N.proj; 0054 N.egrad2rgrad = N.proj; 0055 N.retr = @retraction; 0056 N.exp = N.retr; 0057 function yy = retraction(y, u, t) 0058 if nargin == 2 0059 t = 1; 0060 end 0061 y_plus_tu = M.lincomb(x, 1, y, t, u); 0062 % Mathematically, y_plus_tu is exactly in the tangent space to M at 0063 % x. However, numerically, it may 'leave' the tangent space 0064 % slightly. The extra 'projection' on the next line is not required 0065 % mathematically but helps prevent numerical issues sometimes. 0066 % If this proves to be a huge slow down, one could consider adding 0067 % a type of counter that only executes this extra step every so 0068 % often, instead of at every call. 0069 y_plus_tu = M.tangent(x, y_plus_tu); 0070 nrm = M.norm(x, y_plus_tu); 0071 yy = M.lincomb(x, 1/nrm, y_plus_tu); 0072 end 0073 N.rand = @random; 0074 function y = random() 0075 y = M.randvec(x); 0076 nrm = M.norm(x, y); 0077 y = M.lincomb(x, 1/nrm, y); 0078 end 0079 N.randvec = @randvec; 0080 function u = randvec(y) 0081 u = N.proj(y, M.randvec(x)); 0082 nrm = N.norm(y, u); 0083 u = M.lincomb(x, 1/nrm, u); 0084 end 0085 N.zerovec = @(y) M.zerovec(x); 0086 N.transp = @(y1, y2, u) N.proj(y2, u); 0087 N.hash = @(y) ['z' hashmd5(M.vec(x, y))]; 0088 0089 N.lincomb = @Nlincomb; 0090 function v = Nlincomb(y, a1, d1, a2, d2) %#ok<INUSL> 0091 if nargin == 3 0092 v = M.lincomb(x, a1, d1); 0093 elseif nargin == 5 0094 v = M.lincomb(x, a1, d1, a2, d2); 0095 else 0096 error('lincomb takes either 3 or 5 inputs.'); 0097 end 0098 end 0099 0100 end