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.

 See also: stiefelfactory obliquefactory multiprod multitransp

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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 %
0030 % See also: stiefelfactory obliquefactory multiprod multitransp
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     
0091     M.egrad2rgrad = M.proj;
0092     
0093     M.ehess2rhess = @ehess2rhess;
0094     function rhess = ehess2rhess(Y, egrad, ehess, Ydot)
0095         Y3 = to3D(Y);
0096         Ydot3 = to3D(Ydot);
0097         egrad3 = to3D(egrad);
0098         C = symbdiag(Y3, egrad3);
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