Returns a manifold structure to optimize over unitary matrices. function M = unitaryfactory(n) function M = unitaryfactory(n, k) Unitary group: deals with arrays U of size n x n x k (or n x n if k = 1, which is the default) such that each n x n matrix is unitary, that is, X'*X = eye(n) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(n) for i = 1 : k if k > 1. This is a description of U(n)^k with the induced metric from the embedding space (C^nxn)^k, i.e., this manifold is a Riemannian submanifold of (C^nxn)^k endowed with the usual real inner product on C^nxn, namely, <A, B> = real(trace(A'*B)). This is important: Tangent vectors are represented in the Lie algebra, i.e., as skew-Hermitian matrices. Use the function M.tangent2ambient(X, H) to switch from the Lie algebra representation to the embedding space representation. This is often necessary to define problem.ehess(X, H), as the input H will then be a skew-Hermitian matrix (but the output must not be, as the output is the Hessian in the embedding Euclidean space.) By default, the retraction is only a first-order approximation of the exponential. To force the use of a second-order approximation, call M.retr = M.retr_polar after creating M. This switches from a QR-based computation to an SVD-based computation. By default, k = 1. See also: stiefelcomplexfactory rotationsgroup stiefelfactory
0001 function M = unitaryfactory(n, k) 0002 % Returns a manifold structure to optimize over unitary matrices. 0003 % 0004 % function M = unitaryfactory(n) 0005 % function M = unitaryfactory(n, k) 0006 % 0007 % Unitary group: deals with arrays U of size n x n x k (or n x n if k = 1, 0008 % which is the default) such that each n x n matrix is unitary, that is, 0009 % X'*X = eye(n) if k = 1, or 0010 % X(:, :, i)' * X(:, :, i) = eye(n) for i = 1 : k if k > 1. 0011 % 0012 % This is a description of U(n)^k with the induced metric from the 0013 % embedding space (C^nxn)^k, i.e., this manifold is a Riemannian 0014 % submanifold of (C^nxn)^k endowed with the usual real inner product on 0015 % C^nxn, namely, <A, B> = real(trace(A'*B)). 0016 % 0017 % This is important: 0018 % Tangent vectors are represented in the Lie algebra, i.e., as 0019 % skew-Hermitian matrices. Use the function M.tangent2ambient(X, H) to 0020 % switch from the Lie algebra representation to the embedding space 0021 % representation. This is often necessary to define problem.ehess(X, H), 0022 % as the input H will then be a skew-Hermitian matrix (but the output must 0023 % not be, as the output is the Hessian in the embedding Euclidean space.) 0024 % 0025 % By default, the retraction is only a first-order approximation of the 0026 % exponential. To force the use of a second-order approximation, call 0027 % M.retr = M.retr_polar after creating M. This switches from a QR-based 0028 % computation to an SVD-based computation. 0029 % 0030 % By default, k = 1. 0031 % 0032 % See also: stiefelcomplexfactory rotationsgroup stiefelfactory 0033 0034 % This file is part of Manopt: www.manopt.org. 0035 % Original author: Nicolas Boumal, June 18, 2019. 0036 % Contributors: 0037 % Change log: 0038 0039 0040 if ~exist('k', 'var') || isempty(k) 0041 k = 1; 0042 end 0043 0044 if k == 1 0045 M.name = @() sprintf('Unitary manifold U(%d)', n); 0046 elseif k > 1 0047 M.name = @() sprintf('Product unitary manifold U(%d)^%d', n, k); 0048 else 0049 error('k must be an integer no less than 1.'); 0050 end 0051 0052 M.dim = @() k*(n^2); 0053 0054 M.inner = @(x, d1, d2) real(d1(:)'*d2(:)); 0055 0056 M.norm = @(x, d) norm(d(:)); 0057 0058 M.typicaldist = @() pi*sqrt(n*k); 0059 0060 M.proj = @(X, H) multiskewh(multiprod(multihconj(X), H)); 0061 0062 M.tangent = @(X, H) multiskewh(H); 0063 0064 M.tangent2ambient_is_identity = false; 0065 M.tangent2ambient = @(X, U) multiprod(X, U); 0066 0067 M.egrad2rgrad = M.proj; 0068 0069 M.ehess2rhess = @ehess2rhess; 0070 function Rhess = ehess2rhess(X, Egrad, Ehess, H) 0071 % Reminder : H contains skew-Hermitian matrices. The actual 0072 % direction that the point X is moved along is X*H. 0073 Xt = multihconj(X); 0074 XtEgrad = multiprod(Xt, Egrad); 0075 symXtEgrad = multiherm(XtEgrad); 0076 XtEhess = multiprod(Xt, Ehess); 0077 Rhess = multiskewh( XtEhess - multiprod(H, symXtEgrad) ); 0078 end 0079 0080 % This QR-based retraction is only a first-order approximation 0081 % of the exponential map, not a second-order one. 0082 M.retr_qr = @retraction_qr; 0083 function Y = retraction_qr(X, U, t) 0084 % It is necessary to call qr_unique rather than simply qr to ensure 0085 % this is a retraction, to avoid spurious column sign flips. 0086 XU = multiprod(X, U); 0087 if nargin < 3 0088 Y = qr_unique(X + XU); 0089 else 0090 Y = qr_unique(X + t*XU); 0091 end 0092 end 0093 0094 % A second-order retraction is implemented here. To force its use, 0095 % after creating the factory M, execute M.retr = M.retr_polar. 0096 M.retr_polar = @retraction_polar; 0097 function Y = retraction_polar(X, U, t) 0098 if nargin == 3 0099 tU = t*U; 0100 else 0101 tU = U; 0102 end 0103 Y = X + multiprod(X, tU); 0104 for kk = 1 : k 0105 [Uk, ~, Vk] = svd(Y(:, :, kk)); 0106 Y(:, :, kk) = Uk*Vk'; 0107 end 0108 end 0109 0110 % By default, use QR retraction 0111 M.retr = M.retr_qr; 0112 0113 M.exp = @exponential; 0114 function Y = exponential(X, U, t) 0115 if nargin == 3 0116 exptU = t*U; 0117 else 0118 exptU = U; 0119 end 0120 for kk = 1 : k 0121 exptU(:, :, kk) = expm(exptU(:, :, kk)); 0122 end 0123 Y = multiprod(X, exptU); 0124 end 0125 0126 M.log = @logarithm; 0127 function U = logarithm(X, Y) 0128 U = multiprod(multihconj(X), Y); 0129 for kk = 1 : k 0130 U(:, :, kk) = logm(U(:, :, kk)); 0131 end 0132 % Ensure the tangent vector is in the Lie algebra. 0133 U = multiskewh(U); 0134 end 0135 0136 M.hash = @(X) ['z' hashmd5([real(X(:)) ; imag(X(:))])]; 0137 0138 M.rand = @() randunitary(n, k); 0139 0140 M.randvec = @randomvec; 0141 function U = randomvec(X) %#ok<INUSD> 0142 U = randskewh(n, k); 0143 nrmU = sqrt(U(:)'*U(:)); 0144 U = U / nrmU; 0145 end 0146 0147 M.lincomb = @matrixlincomb; 0148 0149 M.zerovec = @(x) zeros(n, n, k); 0150 0151 M.transp = @(x1, x2, d) d; 0152 M.isotransp = M.transp; % the transport is isometric 0153 0154 M.pairmean = @pairmean; 0155 function Y = pairmean(X1, X2) 0156 V = M.log(X1, X2); 0157 Y = M.exp(X1, .5*V); 0158 end 0159 0160 M.dist = @(x, y) M.norm(x, M.log(x, y)); 0161 0162 sz = n*n*k; 0163 M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))]; 0164 M.mat = @(x, u_vec) reshape(u_vec(1:sz), [n, n, k]) ... 0165 + 1i*reshape(u_vec((sz+1):end), [n, n, k]); 0166 M.vecmatareisometries = @() true; 0167 M.lie_identity = @() repmat(eye(n), [1, 1, k]); 0168 0169 end 0170