Home > manopt > manifolds > stiefel > stiefelgeneralizedfactory.m

stiefelgeneralizedfactory

PURPOSE ^

Returns a manifold structure of "scaled" orthonormal matrices.

SYNOPSIS ^

function M = stiefelgeneralizedfactory(n, p, B)

DESCRIPTION ^

 Returns a manifold structure of "scaled" orthonormal matrices.

 function M = stiefelgeneralizedfactory(n, p)
 function M = stiefelgeneralizedfactory(n, p, B)

 The generalized Stiefel manifold is the set of "scaled" orthonormal 
 nxp matrices X such that X'*B*X is identity. B must be positive definite.
 If B is identity, then this is the standard Stiefel manifold.

 The generalized Stiefel manifold is endowed with a scaled metric
 by making it a Riemannian submanifold of the Euclidean space,
 again endowed with the scaled inner product.

 Some notions (not all) are from Section 4.5 of the paper
 "The geometry of algorithms with orthogonality constraints",
 A. Edelman, T. A. Arias, S. T. Smith, SIMAX, 1998.

 Paper link: http://arxiv.org/abs/physics/9806030.

 Note: egrad2rgrad and ehess2rhess involve solving linear systems in B. If
 this is a bottleneck for a specific application, then a way forward is to
 create a modified version of this file which preprocesses B to speed this
 up (typically, by computing a Cholesky factorization of it, then calling
 an appropriate solver).

 See also: stiefelfactory  grassmannfactory  grassmanngeneralizedfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = stiefelgeneralizedfactory(n, p, B)
0002 % Returns a manifold structure of "scaled" orthonormal matrices.
0003 %
0004 % function M = stiefelgeneralizedfactory(n, p)
0005 % function M = stiefelgeneralizedfactory(n, p, B)
0006 %
0007 % The generalized Stiefel manifold is the set of "scaled" orthonormal
0008 % nxp matrices X such that X'*B*X is identity. B must be positive definite.
0009 % If B is identity, then this is the standard Stiefel manifold.
0010 %
0011 % The generalized Stiefel manifold is endowed with a scaled metric
0012 % by making it a Riemannian submanifold of the Euclidean space,
0013 % again endowed with the scaled inner product.
0014 %
0015 % Some notions (not all) are from Section 4.5 of the paper
0016 % "The geometry of algorithms with orthogonality constraints",
0017 % A. Edelman, T. A. Arias, S. T. Smith, SIMAX, 1998.
0018 %
0019 % Paper link: http://arxiv.org/abs/physics/9806030.
0020 %
0021 % Note: egrad2rgrad and ehess2rhess involve solving linear systems in B. If
0022 % this is a bottleneck for a specific application, then a way forward is to
0023 % create a modified version of this file which preprocesses B to speed this
0024 % up (typically, by computing a Cholesky factorization of it, then calling
0025 % an appropriate solver).
0026 %
0027 % See also: stiefelfactory  grassmannfactory  grassmanngeneralizedfactory
0028 
0029 % This file is part of Manopt: www.manopt.org.
0030 % Original author: Bamdev Mishra, June 30, 2015.
0031 % Contributors:
0032 %
0033 % Change log:
0034 %
0035 %    Sep. 6, 2018 (NB):
0036 %        Removed M.exp() as it was not implemented.
0037 
0038     
0039     if ~exist('B', 'var') || isempty(B)
0040         B = speye(n); % Standard Stiefel manifold.
0041     end
0042     
0043     M.name = @() sprintf('Generalized Stiefel manifold St(%d, %d)', n, p);
0044     
0045     M.dim = @() (n*p - .5*p*(p+1));
0046     
0047     M.inner = @(X, eta, zeta) trace(eta'*(B*zeta)); % Scaled metric.
0048     
0049     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0050     
0051     M.dist = @(X, Y) error('stiefelgeneralizedfactory.dist not implemented yet.');
0052     
0053     M.typicaldist = @() sqrt(p);
0054     
0055     % Orthogonal projection of an ambient vector U to the tangent space
0056     % at X.
0057     M.proj = @projection;
0058     function Up = projection(X, U)
0059         BX = B*X;
0060         
0061         % Projection onto the tangent space
0062         Up = U - X*symm(BX'*U);  
0063     end
0064     
0065     M.tangent = M.proj;
0066     
0067     M.egrad2rgrad = @egrad2rgrad;
0068     function rgrad = egrad2rgrad(X, egrad)
0069         
0070         % First, scale egrad according the to the scaled metric in the
0071         % Euclidean space.
0072         egrad_scaled = B\egrad;
0073         
0074         % Second, project onto the tangent space.
0075         % rgrad = egrad_scaled - X*symm((B*X)'*egrad_scaled);
0076         %
0077         % Verify that symm(BX'*egrad_scaled) = symm(X'*egrad).
0078         
0079         rgrad = egrad_scaled - X*symm(X'*egrad);
0080     end
0081     
0082     
0083     
0084     M.ehess2rhess = @ehess2rhess;
0085     function rhess = ehess2rhess(X, egrad, ehess, H)
0086         egraddot = ehess;
0087         Xdot = H;
0088         
0089         % Directional derivative of the Riemannian gradient.
0090         egrad_scaleddot = B\egraddot;
0091         rgraddot = egrad_scaleddot - Xdot*symm(X'*egrad)...
0092             - X*symm(Xdot'*egrad)...
0093             - X*symm(X'*egraddot);
0094         
0095         % Project onto the tangent space.
0096         rhess = M.proj(X, rgraddot);
0097     end
0098     
0099     
0100     M.retr = @retraction;
0101     function Y = retraction(X, U, t)
0102         if nargin < 3
0103             t = 1.0;
0104         end
0105         Y = guf(X + t*U); % Ensure that Y'*B*Y is identity.
0106     end
0107 
0108 
0109     M.hash = @(X) ['z' hashmd5(X(:))];
0110     
0111     M.rand = @random;
0112     function X = random()
0113         X = guf(randn(n, p)); % Ensure that X'*B*X is identity;
0114     end
0115     
0116     M.randvec = @randomvec;
0117     function U = randomvec(X)
0118         U = projection(X, randn(n, p));
0119         U = U / norm(U(:));
0120     end
0121     
0122     M.lincomb = @matrixlincomb;
0123     
0124     M.zerovec = @(X) zeros(n, p);
0125     
0126     % This transport is compatible with the generalized polar retraction.
0127     M.transp = @(X1, X2, d) projection(X2, d);
0128     
0129     M.vec = @(X, u_mat) u_mat(:);
0130     M.mat = @(X, u_vec) reshape(u_vec, [n, p]);
0131     M.vecmatareisometries = @() false;
0132     
0133     % Some auxiliary functions
0134     symm = @(D) (D + D')/2;
0135     
0136     function X = guf(Y)
0137         % Generalized polar decomposition of an n-by-p matrix Y.
0138         % X'*B*X is identity.
0139         
0140         % Method 1
0141         [u, ~, v] = svd(Y, 0);
0142   
0143         % Instead of the following three steps, an equivalent, but an
0144         % expensive way is to do X = u*(sqrtm(u'*(B*u))\(v')).
0145         [q, ssquare] = eig(u'*(B*u));
0146         qsinv = q/sparse(diag(sqrt(diag(ssquare))));
0147         X = u*((qsinv*q')*v'); % X'*B*X is identity.
0148         
0149         
0150         % Another computation using restricted_svd
0151         % [u, ~, v] = restricted_svd(Y);
0152         % X = u*v'; % X'*B*X is identity.
0153         
0154     end
0155     
0156     function [u, s, v] = restricted_svd(Y)
0157         % We compute a thin svd-like decomposition of an n-by-p matrix Y
0158         % into matrices u, s, and v such that u is an n-by-p matrix
0159         % with u'*B*u being identity, s is a p-by-p diagonal matrix
0160         % with positive entries, and v is a p-by-p orthogonal matrix.
0161         % Y = u*s*v'.
0162         [v, ssquare] = eig(symm(Y'*(B*Y))); % Y*B*Y is positive definite
0163         ssquarevec = diag(ssquare);
0164         
0165         s = sparse(diag(abs(sqrt(ssquarevec))));
0166         u = Y*(v/s); % u'*B*u is identity.
0167     end
0168 
0169 end

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