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 
0036     
0037     if ~exist('B', 'var') || isempty(B)
0038         B = speye(n); % Standard Stiefel manifold.
0039     end
0040     
0041     M.name = @() sprintf('Generalized Stiefel manifold St(%d, %d)', n, p);
0042     
0043     M.dim = @() (n*p - .5*p*(p+1));
0044     
0045     M.inner = @(X, eta, zeta) trace(eta'*(B*zeta)); % Scaled metric.
0046     
0047     M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
0048     
0049     M.dist = @(X, Y) error('stiefelgeneralizedfactory.dist not implemented yet.');
0050     
0051     M.typicaldist = @() sqrt(p);
0052     
0053     % Orthogonal projection of an ambient vector U to the tangent space
0054     % at X.
0055     M.proj = @projection;
0056     function Up = projection(X, U)
0057         BX = B*X;
0058         
0059         % Projection onto the tangent space
0060         Up = U - X*symm(BX'*U);  
0061     end
0062     
0063     M.tangent = M.proj;
0064     
0065     M.egrad2rgrad = @egrad2rgrad;
0066     function rgrad = egrad2rgrad(X, egrad)
0067         
0068         % First, scale egrad according the to the scaled metric in the
0069         % Euclidean space.
0070         egrad_scaled = B\egrad;
0071         
0072         % Second, project onto the tangent space.
0073         % rgrad = egrad_scaled - X*symm((B*X)'*egrad_scaled);
0074         %
0075         % Verify that symm(BX'*egrad_scaled) = symm(X'*egrad).
0076         
0077         rgrad = egrad_scaled - X*symm(X'*egrad);
0078     end
0079     
0080     
0081     
0082     M.ehess2rhess = @ehess2rhess;
0083     function rhess = ehess2rhess(X, egrad, ehess, H)
0084         egraddot = ehess;
0085         Xdot = H;
0086         
0087         % Directional derivative of the Riemannian gradient.
0088         egrad_scaleddot = B\egraddot;
0089         rgraddot = egrad_scaleddot - Xdot*symm(X'*egrad)...
0090             - X*symm(Xdot'*egrad)...
0091             - X*symm(X'*egraddot);
0092         
0093         % Project onto the tangent space.
0094         rhess = M.proj(X, rgraddot);
0095     end
0096     
0097     
0098     M.retr = @retraction;
0099     function Y = retraction(X, U, t)
0100         if nargin < 3
0101             t = 1.0;
0102         end
0103         Y = guf(X + t*U); % Ensure that Y'*B*Y is identity.
0104     end
0105     
0106     
0107     M.exp = @exponential;
0108     function Y = exponential(X, Z, t)
0109         if nargin < 3
0110             t = 1.0;
0111         end
0112         Y = retraction(X, Z, t);
0113         warning('manopt:stiefelgeneralizedfactory:exp', ...
0114                ['Exponential for generalized Stiefel manifold ' ...
0115                 'manifold not implemented yet. Used retraction instead.']);
0116     end
0117 
0118 
0119     M.hash = @(X) ['z' hashmd5(X(:))];
0120     
0121     M.rand = @random;
0122     function X = random()
0123         X = guf(randn(n, p)); % Ensure that X'*B*X is identity;
0124     end
0125     
0126     M.randvec = @randomvec;
0127     function U = randomvec(X)
0128         U = projection(X, randn(n, p));
0129         U = U / norm(U(:));
0130     end
0131     
0132     M.lincomb = @matrixlincomb;
0133     
0134     M.zerovec = @(X) zeros(n, p);
0135     
0136     % This transport is compatible with the generalized polar retraction.
0137     M.transp = @(X1, X2, d) projection(X2, d);
0138     
0139     M.vec = @(X, u_mat) u_mat(:);
0140     M.mat = @(X, u_vec) reshape(u_vec, [n, p]);
0141     M.vecmatareisometries = @() false;
0142     
0143     % Some auxiliary functions
0144     symm = @(D) (D + D')/2;
0145     
0146     function X = guf(Y)
0147         % Generalized polar decomposition of an n-by-p matrix Y.
0148         % X'*B*X is identity.
0149         
0150         % Method 1
0151         [u, ~, v] = svd(Y, 0);
0152   
0153         % Instead of the following three steps, an equivalent, but an
0154         % expensive way is to do X = u*(sqrtm(u'*(B*u))\(v')).
0155         [q, ssquare] = eig(u'*(B*u));
0156         qsinv = q/sparse(diag(sqrt(diag(ssquare))));
0157         X = u*((qsinv*q')*v'); % X'*B*X is identity.
0158         
0159         
0160         % Another computation using restricted_svd
0161         % [u, ~, v] = restricted_svd(Y);
0162         % X = u*v'; % X'*B*X is identity.
0163         
0164     end
0165     
0166     function [u, s, v] = restricted_svd(Y)
0167         % We compute a thin svd-like decomposition of an n-by-p matrix Y
0168         % into matrices u, s, and v such that u is an n-by-p matrix
0169         % with u'*B*u being identity, s is a p-by-p diagonal matrix
0170         % with positive entries, and v is a p-by-p orthogonal matrix.
0171         % Y = u*s*v'.
0172         [v, ssquare] = eig(symm(Y'*(B*Y))); % Y*B*Y is positive definite
0173         ssquarevec = diag(ssquare);
0174         
0175         s = sparse(diag(abs(sqrt(ssquarevec))));
0176         u = Y*(v/s); % u'*B*u is identity.
0177     end
0178 
0179 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005