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