Home > manopt > manifolds > grassmann > grassmanncomplexfactory.m

grassmanncomplexfactory

PURPOSE ^

Returns a manifold struct to optimize over the set of subspaces in C^n.

SYNOPSIS ^

function M = grassmanncomplexfactory(n, p, k)

DESCRIPTION ^

 Returns a manifold struct to optimize over the set of subspaces in C^n.

 function M = grassmanncomplexfactory(n, p)
 function M = grassmanncomplexfactory(n, p, k)

 Complex Grassmann manifold: each point on this manifold is a collection
 of k vector subspaces of dimension p embedded in C^n.

 The metric is obtained by making the Grassmannian a Riemannian quotient
 manifold of the complex Stiefel manifold, i.e., the manifold of
 orthonormal matrices, itself endowed with a metric by making it a
 Riemannian submanifold of the Euclidean space, endowed with the usual
 real-trace inner product, that is, it is the usual metric for the complex
 plane identified with R^2.
 
 This structure deals with complex matrices X of size n x p x k
 (or n x p if k = 1, which is the default) such that each n x p matrix is
 orthonormal, i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) =
 eye(p) for i = 1 : k if k > 1. Each n x p matrix is a numerical
 representation of the vector subspace its columns span.

 By default, k = 1.

 See also: grassmannfactory, stiefelcomplexfactory, grassmanngeneralizedfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = grassmanncomplexfactory(n, p, k)
0002 % Returns a manifold struct to optimize over the set of subspaces in C^n.
0003 %
0004 % function M = grassmanncomplexfactory(n, p)
0005 % function M = grassmanncomplexfactory(n, p, k)
0006 %
0007 % Complex Grassmann manifold: each point on this manifold is a collection
0008 % of k vector subspaces of dimension p embedded in C^n.
0009 %
0010 % The metric is obtained by making the Grassmannian a Riemannian quotient
0011 % manifold of the complex Stiefel manifold, i.e., the manifold of
0012 % orthonormal matrices, itself endowed with a metric by making it a
0013 % Riemannian submanifold of the Euclidean space, endowed with the usual
0014 % real-trace inner product, that is, it is the usual metric for the complex
0015 % plane identified with R^2.
0016 %
0017 % This structure deals with complex matrices X of size n x p x k
0018 % (or n x p if k = 1, which is the default) such that each n x p matrix is
0019 % orthonormal, i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) =
0020 % eye(p) for i = 1 : k if k > 1. Each n x p matrix is a numerical
0021 % representation of the vector subspace its columns span.
0022 %
0023 % By default, k = 1.
0024 %
0025 % See also: grassmannfactory, stiefelcomplexfactory, grassmanngeneralizedfactory
0026 
0027 % This file is part of Manopt: www.manopt.org.
0028 % Original author: Hiroyuki Sato, May 21, 2015.
0029 % Contributors:
0030 % Change log:
0031 
0032     assert(n >= p, ...
0033            ['The dimension n of the ambient space must be larger ' ...
0034             'than the dimension p of the subspaces.']);
0035     
0036     if ~exist('k', 'var') || isempty(k)
0037         k = 1;
0038     end
0039     
0040     if k == 1
0041         M.name = @() sprintf('Complex Grassmann manifold Gr(%d, %d)', n, p);
0042     elseif k > 1
0043         M.name = @() sprintf(['Multi complex Grassmann manifold ' ...
0044             'Gr(%d, %d)^%d'], n, p, k);
0045     else
0046         error('k must be an integer no less than 1.');
0047     end
0048     
0049     M.dim = @() 2*k*p*(n-p); %! k*p*(n-p) -> 2*k*p*(n-p)
0050     
0051     M.inner = @(x, d1, d2) real(d1(:)'*d2(:)); %! trace -> real-trace
0052     
0053     M.norm = @(x, d) norm(d(:));
0054     
0055     M.dist = @distance;
0056     function d = distance(x, y)
0057         principal_angles = zeros(p, k);
0058         XHY = multiprod(multihconj(x), y); %! XtY -> XHY, multitransp -> multihconj
0059         for i = 1 : k
0060             cos_princ_angle = svd(XHY(:, :, i));
0061             principal_angles(:, i) = acos(cos_princ_angle);
0062         end
0063         d = norm(real(principal_angles), 'fro');
0064     end
0065     
0066     M.typicaldist = @() sqrt(p*k);
0067     
0068     % Orthogonal projection of an ambient vector U to the horizontal space
0069     % at X.
0070     M.proj = @projection;
0071     function Up = projection(X, U)
0072         
0073         XHU = multiprod(multihconj(X), U); %! XtU -> XHU, multitransp -> multihconj
0074         Up = U - multiprod(X, XHU); %! XtU -> XHU
0075 
0076     end
0077     
0078     M.tangent = M.proj;
0079     
0080     M.egrad2rgrad = M.proj;
0081     
0082     M.ehess2rhess = @ehess2rhess;
0083     function rhess = ehess2rhess(X, egrad, ehess, H)
0084         PXehess = projection(X, ehess);
0085         XHG = multiprod(multihconj(X), egrad); %! XtG -> XHG, multitransp -> multihconj
0086         HXHG = multiprod(H, XHG); %! HXtG -> HXHG, XtG -> XHG
0087         rhess = PXehess - HXHG; %! HXtG -> HXHG
0088     end
0089     
0090     M.retr = @retraction;
0091     function Y = retraction(X, U, t)
0092         if nargin < 3
0093             t = 1.0;
0094         end
0095         Y = X + t*U;
0096         for i = 1 : k 
0097         
0098             % Compute the polar factorization of Y = X+tU
0099             [u, s, v] = svd(Y(:, :, i), 'econ'); %#ok
0100             Y(:, :, i) = u*v';
0101             
0102             % Another popular retraction uses QR instead of SVD.
0103             % As compared with the Stiefel factory, we do not need to
0104             % worry about flipping signs of columns here, since only
0105             % the column space is important, not the actual columns.
0106             % [Q, unused] = qr(Y(:, :, i), 0); %#ok
0107             % Y(:, :, i) = Q;
0108             
0109         end
0110     end
0111     
0112     M.exp = @exponential;
0113     function Y = exponential(X, U, t)
0114         if nargin == 3
0115             tU = t*U;
0116         else
0117             tU = U;
0118         end
0119         Y = zeros(size(X));
0120         for i = 1 : k
0121             [u, s, v] = svd(tU(:, :, i), 0);
0122             cos_s = diag(cos(diag(s)));
0123             sin_s = diag(sin(diag(s)));
0124             Y(:, :, i) = X(:, :, i)*v*cos_s*v' + u*sin_s*v';
0125             % From numerical experiments, it seems necessary to
0126             % re-orthonormalize. This is overall quite expensive.
0127             [q, unused] = qr(Y(:, :, i), 0); %#ok
0128             Y(:, :, i) = q;
0129         end
0130     end
0131 
0132     % Test code for the logarithm:
0133     % Gr = grassmanncomplexfactory(5, 2, 3);
0134     % x = Gr.rand()
0135     % y = Gr.rand()
0136     % u = Gr.log(x, y)
0137     % Gr.dist(x, y) % These two numbers should
0138     % Gr.norm(x, u) % be the same.
0139     % z = Gr.exp(x, u) % z needs not be the same matrix as y, but it should
0140     % v = Gr.log(x, z) % be the same point as y on Grassmann: dist almost 0.
0141     M.log = @logarithm;
0142     function U = logarithm(X, Y)
0143         U = zeros(n, p, k);
0144         for i = 1 : k
0145             x = X(:, :, i);
0146             y = Y(:, :, i);
0147             yHx = y'*x; %! ytx -> yHx, y.' -> y'
0148             AH = y'-yHx*x'; %! At -> AH, x.' -> x', y.' -> y'
0149             BH = yHx\AH; %! Bt -> BH, ytx -> yHx, At -> AH
0150             [u, s, v] = svd(BH', 'econ'); %! Bt.' -> BH'
0151 
0152             u = u(:, 1:p);
0153             s = diag(s);
0154             s = s(1:p);
0155             v = v(:, 1:p);
0156 
0157             U(:, :, i) = u*diag(atan(s))*v'; %! v.' -> v'
0158         end
0159     end
0160 
0161     M.hash = @(X) ['z' hashmd5([real(X(:)); imag(X(:))])]; %! X(:) -> [real(X(:)); imag(X(:))]
0162     
0163     M.rand = @random;
0164     function X = random()
0165         X = zeros(n, p, k);
0166         for j = 1 : k
0167             [Q, unused] = qr(randn(n, p) + 1i*randn(n, p), 0); %#ok<NASGU> %! Complex version
0168             X(:, :, j) = Q;
0169         end
0170     end
0171     
0172     M.randvec = @randomvec;
0173     function U = randomvec(X)
0174         U = projection(X, randn(n, p, k) + 1i*randn(n, p, k)); %! Complex version
0175         U = U / norm(U(:));
0176     end
0177     
0178     M.lincomb = @matrixlincomb;
0179     
0180     M.zerovec = @(x) zeros(n, p, k);
0181     
0182     % This transport is compatible with the polar retraction.
0183     M.transp = @(x1, x2, d) projection(x2, d);
0184     
0185     M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))];
0186     M.mat = @(x, u_vec) reshape(u_vec(1:(n*p*k)) + 1i*u_vec((n*p*k+1):end), [n, p, k]);
0187     M.vecmatareisometries = @() true;
0188 
0189 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005