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:
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multiprod Matrix multiply 2-D slices of N-D arrays
• multisym Returns the symmetric parts of the matrices in a 3D array
• multitransp Transpose the matrix slices of an N-D array (no complex conjugate)
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             % Alternatively, one could also use qr_unique as retraction.
0115         end
0116         Y = to2D(Y3);
0117     end
0118
0119     M.exp = @exponential;
0120     function Y = exponential(Y, U, t)
0121         if nargin == 2
0122             t = 1;
0123         end
0124         tU3 = multitransp(to3D(t*U));
0125         Y3 = multitransp(to3D(Y));
0126         % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
0127         for i = 1 : m
0128             X = Y3(:, :, i);
0129             Z = tU3(:, :, i);
0130             Y3(:, :, i) = [X, Z] * ...
0131                           expm([  X'*Z , -Z'*Z ; eye(d) , X'*Z]) * ...
0132                           [ expm(-X'*Z) ; zeros(d) ];
0133             % We may loose orthonormality here. Just to be sure:
0134             [u, s, v] = svd(Y3(:, :, i), 'econ'); %#ok<ASGLU>
0135             Y3(:, :, i) = u*v';
0136         end
0137         Y = to2D(multitransp(Y3));
0138     end
0139
0140     M.hash = @(Y) ['z' hashmd5(Y(:))];
0141
0142     M.rand = @random;
0143     function Y = random()
0144         Y3 = zeros(d, k, m);
0145         for i = 1 : m
0146             [Q, unused] = qr(randn(k, d), 0); %#ok<ASGLU>
0147             Y3(:, :, i) = Q';
0148         end
0149         Y = to2D(Y3);
0150     end
0151
0152     M.randvec = @randomvec;
0153     function U = randomvec(Y)
0154         U = projection(Y, randn(n, k));
0155         U = U / M.norm(Y, U);
0156     end
0157
0158     M.lincomb = @matrixlincomb;
0159
0160     M.zerovec = @(x) zeros(n, k);
0161
0162     M.transp = @(x1, x2, u) projection(x2, u);
0163
0164     M.vec = @(x, u_mat) u_mat(:);
0165     M.mat = @(x, u_vec) reshape(u_vec, [n, k]);
0166     M.vecmatareisometries = @() true;
0167
0168 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005