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
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
0067
0068         % First, scale egrad according the to the scaled metric in the
0069         % Euclidean space.
0071
0072         % Second, project onto the tangent space.
0074         %
0076
0078     end
0079
0080
0081
0082     M.ehess2rhess = @ehess2rhess;
0083     function rhess = ehess2rhess(X, egrad, ehess, H)
0085         Xdot = H;
0086
0087         % Directional derivative of the Riemannian gradient.
0092
0093         % Project onto the tangent space.
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