Home > manopt > manifolds > grassmann > grassmannfactory.m

# grassmannfactory

## PURPOSE

Returns a manifold struct to optimize over the space of vector subspaces.

## SYNOPSIS

function M = grassmannfactory(n, p, k)

## DESCRIPTION

``` Returns a manifold struct to optimize over the space of vector subspaces.

function M = grassmannfactory(n, p)
function M = grassmannfactory(n, p, k)

Grassmann manifold: each point on this manifold is a collection of k
vector subspaces of dimension p embedded in R^n.

The metric is obtained by making the Grassmannian a Riemannian quotient
manifold of the Stiefel manifold, i.e., the manifold of orthonormal
matrices, itself endowed with a metric by making it a Riemannian
submanifold of the Euclidean space, endowed with the usual inner product.
In short: it is the usual metric used in most cases.

This structure deals with matrices X of size n x p x k (or n x p if
k = 1, which is the default) such that each n x p matrix is orthonormal,
i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
i = 1 : k if k > 1. Each n x p matrix is a numerical representation of
the vector subspace its columns span.

By default, k = 1.

## CROSS-REFERENCE INFORMATION

This function calls:
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multiprod Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
• multitransp Transposing arrays of matrices.
This function is called by:

## SOURCE CODE

```0001 function M = grassmannfactory(n, p, k)
0002 % Returns a manifold struct to optimize over the space of vector subspaces.
0003 %
0004 % function M = grassmannfactory(n, p)
0005 % function M = grassmannfactory(n, p, k)
0006 %
0007 % Grassmann manifold: each point on this manifold is a collection of k
0008 % vector subspaces of dimension p embedded in R^n.
0009 %
0010 % The metric is obtained by making the Grassmannian a Riemannian quotient
0011 % manifold of the Stiefel manifold, i.e., the manifold of orthonormal
0012 % matrices, itself endowed with a metric by making it a Riemannian
0013 % submanifold of the Euclidean space, endowed with the usual inner product.
0014 % In short: it is the usual metric used in most cases.
0015 %
0016 % This structure deals with matrices X of size n x p x k (or n x p if
0017 % k = 1, which is the default) such that each n x p matrix is orthonormal,
0018 % i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
0019 % i = 1 : k if k > 1. Each n x p matrix is a numerical representation of
0020 % the vector subspace its columns span.
0021 %
0022 % By default, k = 1.
0023 %
0025
0026 % This file is part of Manopt: www.manopt.org.
0027 % Original author: Nicolas Boumal, Dec. 30, 2012.
0028 % Contributors:
0029 % Change log:
0030 %   March 22, 2013 (NB) :
0031 %       Implemented geodesic distance.
0032 %
0033 %   April 17, 2013 (NB) :
0034 %       Retraction changed to the polar decomposition, so that the vector
0035 %       transport is now correct, in the sense that it is compatible with
0036 %       the retraction, i.e., transporting a tangent vector G from U to V
0037 %       where V = Retr(U, H) will give Z, and transporting GQ from UQ to VQ
0038 %       will give ZQ: there is no dependence on the representation, which
0039 %       is as it should be. Notice that the polar factorization requires an
0040 %       SVD whereas the qfactor retraction requires a QR decomposition,
0041 %       which is cheaper. Hence, if the retraction happens to be a
0042 %       bottleneck in your application and you are not using vector
0043 %       transports, you may want to replace the retraction with a qfactor.
0044 %
0045 %   July  4, 2013 (NB) :
0046 %       Added support for the logarithmic map 'log'.
0047 %
0048 %   July  5, 2013 (NB) :
0049 %       Added support for ehess2rhess.
0050 %
0051 %   June 24, 2014 (NB) :
0052 %       Small bug fix in the retraction, and added final
0053 %       re-orthonormalization at the end of the exponential map. This
0054 %       follows discussions on the forum where it appeared there is a
0055 %       significant loss in orthonormality without that extra step. Also
0056 %       changed the randvec function so that it now returns a globally
0057 %       normalized vector, not a vector where each component is normalized
0058 %       (this only matters if k>1).
0059
0060     assert(n >= p, ...
0061            ['The dimension n of the ambient space must be larger ' ...
0062             'than the dimension p of the subspaces.']);
0063
0064     if ~exist('k', 'var') || isempty(k)
0065         k = 1;
0066     end
0067
0068     if k == 1
0069         M.name = @() sprintf('Grassmann manifold Gr(%d, %d)', n, p);
0070     elseif k > 1
0071         M.name = @() sprintf('Multi Grassmann manifold Gr(%d, %d)^%d', ...
0072                              n, p, k);
0073     else
0074         error('k must be an integer no less than 1.');
0075     end
0076
0077     M.dim = @() k*p*(n-p);
0078
0079     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0080
0081     M.norm = @(x, d) norm(d(:));
0082
0083     M.dist = @distance;
0084     function d = distance(x, y)
0085         square_d = 0;
0086         XtY = multiprod(multitransp(x), y);
0087         for i = 1 : k
0088             cos_princ_angle = svd(XtY(:, :, i));
0089             square_d = square_d + sum(real(acos(cos_princ_angle)).^2);
0090         end
0091         d = sqrt(square_d);
0092     end
0093
0094     M.typicaldist = @() sqrt(p*k);
0095
0096     % Orthogonal projection of an ambient vector U to the horizontal space
0097     % at X.
0098     M.proj = @projection;
0099     function Up = projection(X, U)
0100
0101         XtU = multiprod(multitransp(X), U);
0102         Up = U - multiprod(X, XtU);
0103
0104     end
0105
0106     M.tangent = M.proj;
0107
0109
0110     M.ehess2rhess = @ehess2rhess;
0111     function rhess = ehess2rhess(X, egrad, ehess, H)
0112         PXehess = projection(X, ehess);
0114         HXtG = multiprod(H, XtG);
0115         rhess = PXehess - HXtG;
0116     end
0117
0118     M.retr = @retraction;
0119     function Y = retraction(X, U, t)
0120         if nargin < 3
0121             t = 1.0;
0122         end
0123         Y = X + t*U;
0124         for i = 1 : k
0125
0126             % Compute the polar factorization of Y = X+tU
0127             [u, s, v] = svd(Y(:, :, i), 'econ'); %#ok
0128             Y(:, :, i) = u*v';
0129
0130             % Another popular retraction uses QR instead of SVD.
0131             % As compared with the Stiefel factory, we do not need to
0132             % worry about flipping signs of columns here, since only
0133             % the column space is important, not the actual columns.
0134             % [Q, unused] = qr(Y(:, :, i), 0); %#ok
0135             % Y(:, :, i) = Q;
0136
0137         end
0138     end
0139
0140     M.exp = @exponential;
0141     function Y = exponential(X, U, t)
0142         if nargin == 3
0143             tU = t*U;
0144         else
0145             tU = U;
0146         end
0147         Y = zeros(size(X));
0148         for i = 1 : k
0149             [u, s, v] = svd(tU(:, :, i), 0);
0150             cos_s = diag(cos(diag(s)));
0151             sin_s = diag(sin(diag(s)));
0152             Y(:, :, i) = X(:, :, i)*v*cos_s*v' + u*sin_s*v';
0153             % From numerical experiments, it seems necessary to
0154             % re-orthonormalize. This is overall quite expensive.
0155             [q, unused] = qr(Y(:, :, i), 0); %#ok
0156             Y(:, :, i) = q;
0157         end
0158     end
0159
0160     % Test code for the logarithm:
0161     % Gr = grassmannfactory(5, 2, 3);
0162     % x = Gr.rand()
0163     % y = Gr.rand()
0164     % u = Gr.log(x, y)
0165     % Gr.dist(x, y) % These two numbers should
0166     % Gr.norm(x, u) % be the same.
0167     % z = Gr.exp(x, u) % z needs not be the same matrix as y, but it should
0168     % v = Gr.log(x, z) % be the same point as y on Grassmann: dist almost 0.
0169     M.log = @logarithm;
0170     function U = logarithm(X, Y)
0171         U = zeros(n, p, k);
0172         for i = 1 : k
0173             x = X(:, :, i);
0174             y = Y(:, :, i);
0175             ytx = y.'*x;
0176             At = y.'-ytx*x.';
0177             Bt = ytx\At;
0178             [u, s, v] = svd(Bt.', 'econ');
0179
0180             u = u(:, 1:p);
0181             s = diag(s);
0182             s = s(1:p);
0183             v = v(:, 1:p);
0184
0185             U(:, :, i) = u*diag(atan(s))*v.';
0186         end
0187     end
0188
0189     M.hash = @(X) ['z' hashmd5(X(:))];
0190
0191     M.rand = @random;
0192     function X = random()
0193         X = zeros(n, p, k);
0194         for i = 1 : k
0195             [Q, unused] = qr(randn(n, p), 0); %#ok<NASGU>
0196             X(:, :, i) = Q;
0197         end
0198     end
0199
0200     M.randvec = @randomvec;
0201     function U = randomvec(X)
0202         U = projection(X, randn(n, p, k));
0203         U = U / norm(U(:));
0204     end
0205
0206     M.lincomb = @matrixlincomb;
0207
0208     M.zerovec = @(x) zeros(n, p, k);
0209
0210     % This transport is compatible with the polar retraction.
0211     M.transp = @(x1, x2, d) projection(x2, d);
0212
0213     M.vec = @(x, u_mat) u_mat(:);
0214     M.mat = @(x, u_vec) reshape(u_vec, [n, p, k]);
0215     M.vecmatareisometries = @() true;
0216
0217 end```

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