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.

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).

## CROSS-REFERENCE INFORMATION This function calls:
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
This function is called by:

## 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 %
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 %
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
0069
0070         % First, scale egrad according the to the scaled metric in the
0071         % Euclidean space.
0073
0074         % Second, project onto the tangent space.
0076         %
0078
0080     end
0081
0082
0083
0084     M.ehess2rhess = @ehess2rhess;
0085     function rhess = ehess2rhess(X, egrad, ehess, H)
0087         Xdot = H;
0088
0089         % Directional derivative of the Riemannian gradient.
0094
0095         % Project onto the tangent space.
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