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