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
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