unitaryfactory

 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