Home > manopt > manifolds > stiefel > stiefelstackedfactory.m

# stiefelstackedfactory

## PURPOSE

Stiefel(k, d)^m, represented as matrices of size m*d-by-k.

## SYNOPSIS

function M = stiefelstackedfactory(m, d, k)

## DESCRIPTION

``` Stiefel(k, d)^m, represented as matrices of size m*d-by-k.

function M = stiefelstackedfactory(m, d, k)

Points on this manifold are matrices Y of size n x k, with n = m*d.
Y is thought of as m matrices of size d x k each, stacked on top of each
other. Call them Y1, ..., Ym. Each Yi is an orthonormal matrix, that is,
its d rows are unit norm and are orthogonal to each other. Thus, this
geometry is a product of Stiefel manifolds.

To easily transform matrices Y to 3D arrays Y3 of size d x k x m such
that each slice Y3(:, :, i) corresponds to one of the matrices Yi, use
the functions

Y3 = M.to3D(Y)   and   Y = M.to2D(Y3).

The ambient space R^(nxk) is endowed with the usual inner product
<A, B> = trace(A'*B). This inner product is restricted to the tangent
spaces of the present manifold, thus making it a Riemannian submanifold
of the Euclidean space R^(nxk). Tangent vectors are represented as
matrices of the same size as Y, and can likewise be converted to 3D
arrays and back using to3D() and to2D().

In dealing with this geometry, especially when dealing with the 3D array
representations of points and tangent vectors, the tools multiprod,
multitransp, multitrace, multiscale etc. available in Manopt are often
useful.

## 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.
• multisym Returns the symmetric parts of the matrices in the 3D matrix X
• multitransp Transposing arrays of matrices.
This function is called by:

## SOURCE CODE

```0001 function M = stiefelstackedfactory(m, d, k)
0002 % Stiefel(k, d)^m, represented as matrices of size m*d-by-k.
0003 %
0004 % function M = stiefelstackedfactory(m, d, k)
0005 %
0006 % Points on this manifold are matrices Y of size n x k, with n = m*d.
0007 % Y is thought of as m matrices of size d x k each, stacked on top of each
0008 % other. Call them Y1, ..., Ym. Each Yi is an orthonormal matrix, that is,
0009 % its d rows are unit norm and are orthogonal to each other. Thus, this
0010 % geometry is a product of Stiefel manifolds.
0011 %
0012 % To easily transform matrices Y to 3D arrays Y3 of size d x k x m such
0013 % that each slice Y3(:, :, i) corresponds to one of the matrices Yi, use
0014 % the functions
0015 %
0016 %    Y3 = M.to3D(Y)   and   Y = M.to2D(Y3).
0017 %
0018 % The ambient space R^(nxk) is endowed with the usual inner product
0019 % <A, B> = trace(A'*B). This inner product is restricted to the tangent
0020 % spaces of the present manifold, thus making it a Riemannian submanifold
0021 % of the Euclidean space R^(nxk). Tangent vectors are represented as
0022 % matrices of the same size as Y, and can likewise be converted to 3D
0023 % arrays and back using to3D() and to2D().
0024 %
0025 % In dealing with this geometry, especially when dealing with the 3D array
0026 % representations of points and tangent vectors, the tools multiprod,
0027 % multitransp, multitrace, multiscale etc. available in Manopt are often
0028 % useful.
0029 %
0031
0032 % This file is part of Manopt: www.manopt.org.
0033 % Original author: Nicolas Boumal, May 4, 2015.
0034 % Contributors:
0035 % Change log:
0036
0037     assert(k >= d, 'k must be at least as large as d.');
0038
0039     n = m*d;
0040
0041     M.name = @() sprintf('Manifold of %d orthonormal matrices of size %dx%d, stacked', m, d, k);
0042
0043     M.dim = @() m*(k*d - .5*d*(d+1));
0044
0045     M.size = @() [m, d, k];
0046
0047     M.inner = @(x, d1, d2) d1(:).'*d2(:);
0048
0049     M.norm = @(x, d) norm(d(:));
0050
0051     M.dist = @(x, y) error('stiefelstackedfactory.dist not implemented yet.');
0052
0053     M.typicaldist = @() sqrt(M.dim());
0054
0055     % Convert a dxkxm matrix to an nxk matrix
0056     M.to2D = @to2D;
0057     function A2 = to2D(A3)
0058         A2 = reshape(multitransp(A3), [k, m*d])';
0059     end
0060
0061     % Convert an nxk matrix to a dxkxm matrix
0062     M.to3D = @to3D;
0063     function A3 = to3D(A2)
0064         A3 = multitransp(reshape(A2', [k, d, m]));
0065     end
0066
0067     % Given 2 3D matrices A and B of size dxkxm, returns a 3D matrix C of
0068     % size dxdxm such that each slice C(:, :, i) is the symmetric part of
0069     % the product A(:, :, i) * B(:, :, i)'. The name is short for
0070     % "symmetric-block-diagonal", because if A and B were transformed to
0071     % their 2D equivalents via to2D, then the output would contain the
0072     % symmetric parts of the diagonal blocks of A*B'.
0073     M.symbdiag = @symbdiag;
0074     function C = symbdiag(A, B)
0075         C = multisym(multiprod(A, multitransp(B)));
0076     end
0077
0078     % Orthogonal projection from the ambient space R^(nxk) to the tangent
0079     % space at X.
0080     M.proj = @projection;
0081     function Zt = projection(Y, Z)
0082         Y3 = to3D(Y);
0083         Z3 = to3D(Z);
0084         Lambda = symbdiag(Y3, Z3);
0085         Zt3 = Z3 - multiprod(Lambda, Y3);
0086         Zt = to2D(Zt3);
0087     end
0088
0089     M.tangent = M.proj;
0090
0092
0093     M.ehess2rhess = @ehess2rhess;
0094     function rhess = ehess2rhess(Y, egrad, ehess, Ydot)
0095         Y3 = to3D(Y);
0096         Ydot3 = to3D(Ydot);
0099         CYdot = to2D(multiprod(C, Ydot3));
0100         rhess = projection(Y, ehess - CYdot);
0101     end
0102
0103     M.retr = @retraction;
0104     function Y = retraction(Y, U, t)
0105         if nargin < 3
0106             t = 1.0;
0107         end
0108         Y = Y + t*U;
0109         Y3 = to3D(Y);
0110         for i = 1 : m
0111             % Orthonormalize the rows of Y3(:, :, i):
0112             [u, s, v] = svd(Y3(:, :, i), 'econ'); %#ok<ASGLU>
0113             Y3(:, :, i) = u*v';
0114             % Alternative code if one desires to use QR instead of SVD.
0115             % The instruction with the signs of R assures we are not
0116             % flipping signs of some columns.
0117             % [Q, R] = qr(Y3(:, :, i)', 0);
0118             % Y3(:, :, i) = (Q * diag(sign(sign(diag(R))+.5)))';
0119         end
0120         Y = to2D(Y3);
0121     end
0122
0123     M.exp = @exponential;
0124     function Y = exponential(Y, U, t)
0125         if nargin == 2
0126             t = 1;
0127         end
0128         tU3 = multitransp(to3D(t*U));
0129         Y3 = multitransp(to3D(Y));
0130         % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
0131         for i = 1 : m
0132             X = Y3(:, :, i);
0133             Z = tU3(:, :, i);
0134             Y3(:, :, i) = [X, Z] * ...
0135                           expm([  X'*Z , -Z'*Z ; eye(d) , X'*Z]) * ...
0136                           [ expm(-X'*Z) ; zeros(d) ];
0137             % We may loose orthonormality here. Just to be sure:
0138             [u, s, v] = svd(Y3(:, :, i), 'econ'); %#ok<ASGLU>
0139             Y3(:, :, i) = u*v';
0140         end
0141         Y = to2D(multitransp(Y3));
0142     end
0143
0144     M.hash = @(Y) ['z' hashmd5(Y(:))];
0145
0146     M.rand = @random;
0147     function Y = random()
0148         Y3 = zeros(d, k, m);
0149         for i = 1 : m
0150             [Q, unused] = qr(randn(k, d), 0); %#ok<NASGU>
0151             Y3(:, :, i) = Q';
0152         end
0153         Y = to2D(Y3);
0154     end
0155
0156     M.randvec = @randomvec;
0157     function U = randomvec(Y)
0158         U = projection(Y, randn(n, k));
0159         U = U / M.norm(Y, U);
0160     end
0161
0162     M.lincomb = @matrixlincomb;
0163
0164     M.zerovec = @(x) zeros(n, k);
0165
0166     M.transp = @(x1, x2, u) projection(x2, u);
0167
0168     M.vec = @(x, u_mat) u_mat(:);
0169     M.mat = @(x, u_vec) reshape(u_vec, [n, k]);
0170     M.vecmatareisometries = @() true;
0171
0172 end```

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