Returns a manifold struct. to optimize over complex orthonormal matrices. function M = stiefelcomplexfactory(n, p) function M = stiefelcomplexfactory(n, p, k) The complex Stiefel manifold is the set of complex orthonormal nxp matrices. If k is larger than 1, this is the Cartesian product of the complex Stiefel manifold taken k times. The metric is such that the manifold is a Riemannian submanifold of C^nxp equipped with the usual real-trace inner product, that is, it is the usual metric for the complex plane identified with R^2. Points are represented as matrices X of size n x p x k (or n x p if k=1, which is the default) such that each complex 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. Tangent vectors are represented as matrices the same size as points. By default, k = 1. Please cite the Manopt paper as well as either of these research papers pertaining to this specific geometry: @InProceedings{sato2013complex, Title = {A complex singular value decomposition algorithm based on the {R}iemannian {N}ewton method}, Author = {Sato, H. and Iwai, T.}, Booktitle = {Decision and Control ({CDC}), 2013 {IEEE} 52nd Annual Conference on}, Year = {2013}, Organization = {IEEE}, Pages = {2972--2978} } @InProceedings{sato2014Riemannian, Title = {{R}iemannian conjugate gradient method for complex singular value decomposition problem}, Author = {Sato, H.}, Booktitle = {Decision and Control ({CDC}), 2014 {IEEE} 53rd Annual Conference on}, Year = {2014}, Organization = {IEEE}, Pages = {5849--5854} } See also: stiefelfactory
0001 function M = stiefelcomplexfactory(n, p, k) 0002 % Returns a manifold struct. to optimize over complex orthonormal matrices. 0003 % 0004 % function M = stiefelcomplexfactory(n, p) 0005 % function M = stiefelcomplexfactory(n, p, k) 0006 % 0007 % The complex Stiefel manifold is the set of complex orthonormal nxp 0008 % matrices. If k is larger than 1, this is the Cartesian product of the 0009 % complex Stiefel manifold taken k times. The metric is such that the 0010 % manifold is a Riemannian submanifold of C^nxp equipped with the usual 0011 % real-trace inner product, that is, it is the usual metric for the complex 0012 % plane identified with R^2. 0013 % 0014 % Points are represented as matrices X of size n x p x k (or n x p if k=1, 0015 % which is the default) such that each complex n x p matrix is orthonormal, 0016 % i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for 0017 % i = 1 : k if k > 1. Tangent vectors are represented as matrices the same 0018 % size as points. 0019 % 0020 % By default, k = 1. 0021 % 0022 % 0023 % Please cite the Manopt paper as well as either of these research papers 0024 % pertaining to this specific geometry: 0025 % @InProceedings{sato2013complex, 0026 % Title = {A complex singular value decomposition algorithm based on the {R}iemannian {N}ewton method}, 0027 % Author = {Sato, H. and Iwai, T.}, 0028 % Booktitle = {Decision and Control ({CDC}), 2013 {IEEE} 52nd Annual Conference on}, 0029 % Year = {2013}, 0030 % Organization = {IEEE}, 0031 % Pages = {2972--2978} 0032 % } 0033 % @InProceedings{sato2014Riemannian, 0034 % Title = {{R}iemannian conjugate gradient method for complex singular value decomposition problem}, 0035 % Author = {Sato, H.}, 0036 % Booktitle = {Decision and Control ({CDC}), 2014 {IEEE} 53rd Annual Conference on}, 0037 % Year = {2014}, 0038 % Organization = {IEEE}, 0039 % Pages = {5849--5854} 0040 % } 0041 % 0042 % 0043 % See also: stiefelfactory 0044 0045 % This file is part of Manopt: www.manopt.org. 0046 % Original author: Hiroyuki Sato, April 27, 2015. 0047 % Contributors: 0048 % Change log: 0049 % June 18, 2019 (NB) : Using qr_unique for retr and rand. 0050 0051 if ~exist('k', 'var') || isempty(k) 0052 k = 1; 0053 end 0054 0055 if k == 1 0056 M.name = @() sprintf('Complex Stiefel manifold St(%d, %d)', n, p); 0057 elseif k > 1 0058 M.name = @() sprintf('Product complex Stiefel manifold St(%d, %d)^%d', n, p, k); 0059 else 0060 error('k must be an integer no less than 1.'); 0061 end 0062 0063 M.dim = @() k*(2*n*p - p^2); %! k*(n*p - .5*p*(p+1)) -> k*(2*n*p - p^2) 0064 0065 M.inner = @(x, d1, d2) real(d1(:)'*d2(:)); %! trace -> real-trace 0066 0067 M.norm = @(x, d) norm(d(:)); 0068 0069 M.dist = @(x, y) error('stiefel.dist not implemented yet.'); 0070 0071 M.typicaldist = @() sqrt(p*k); 0072 0073 M.proj = @projection; 0074 function Up = projection(X, U) 0075 0076 XHU = multiprod(multihconj(X), U); %! XtU -> XHU, multitransp -> multihconj 0077 herXHU = multiherm(XHU); %! symXtU -> herXHU, multisym -> multiherm 0078 Up = U - multiprod(X, herXHU); %! symXtU -> herXHU 0079 0080 end 0081 0082 M.tangent = M.proj; 0083 0084 % For Riemannian submanifolds, converting a Euclidean gradient into a 0085 % Riemannian gradient amounts to an orthogonal projection. 0086 M.egrad2rgrad = M.proj; 0087 0088 M.ehess2rhess = @ehess2rhess; 0089 function rhess = ehess2rhess(X, egrad, ehess, H) 0090 XHG = multiprod(multihconj(X), egrad); %! XtG -> XHG, multitransp -> multihconj 0091 herXHG = multiherm(XHG); %! symXtG -> herXHG, multisym(XtG) -> multiherm(XHG) 0092 HherXHG = multiprod(H, herXHG); %! HsymXtG -> HherXHG, symXtG -> herXHG 0093 rhess = projection(X, ehess - HherXHG); %! HsymXtG -> HherXHG 0094 end 0095 0096 M.retr = @retraction; 0097 function Y = retraction(X, U, t) 0098 % It is necessary to call qr_unique rather than simply qr to ensure 0099 % this is a retraction, to avoid spurious column sign flips. 0100 % This is only a first-order retraction. 0101 if nargin < 3 0102 Y = qr_unique(X + U); 0103 else 0104 Y = qr_unique(X + t*U); 0105 end 0106 end 0107 0108 M.exp = @exponential; 0109 function Y = exponential(X, U, t) 0110 if nargin == 2 0111 t = 1; 0112 end 0113 tU = t*U; 0114 Y = zeros(size(X)); 0115 for i = 1 : k 0116 % From a formula by Ross Lippert, Example 5.4.2 in AMS08. 0117 Xi = X(:, :, i); 0118 Ui = tU(:, :, i); 0119 Y(:, :, i) = [Xi Ui] * ... 0120 expm([Xi'*Ui , -Ui'*Ui ; eye(p) , Xi'*Ui]) * ... 0121 [ expm(-Xi'*Ui) ; zeros(p) ]; 0122 end 0123 0124 end 0125 0126 M.hash = @(X) ['z' hashmd5([real(X(:)) ; imag(X(:))])]; %! X(:) -> [real(X(:)) ; imag(X(:))] 0127 0128 M.rand = @() qr_unique(randn(n, p, k) + 1i*randn(n, p, k)); 0129 0130 M.randvec = @randomvec; 0131 function U = randomvec(X) 0132 U = projection(X, randn(n, p, k) + 1i*randn(n, p, k)); %! Complex version 0133 U = U / norm(U(:)); 0134 end 0135 0136 M.lincomb = @matrixlincomb; 0137 0138 M.zerovec = @(x) zeros(n, p, k); 0139 0140 M.transp = @(x1, x2, d) projection(x2, d); 0141 0142 M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))]; 0143 M.mat = @(x, u_vec) reshape(u_vec(1:(n*p*k)) + 1i*u_vec((n*p*k+1):end), [n, p, k]); 0144 M.vecmatareisometries = @() true; 0145 0146 end