Home > manopt > manifolds > rotations > unitaryfactory.m

# unitaryfactory

## PURPOSE

Returns a manifold structure to optimize over unitary matrices.

## SYNOPSIS

function M = unitaryfactory(n, k)

## DESCRIPTION

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

## CROSS-REFERENCE INFORMATION

This function calls:
• randskewh Generates random skew Hermitian matrices with normal entries.
• randunitary Generates uniformly random unitary matrices.
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multihconj Hermitian-conjugate transpose the matrix slices of an N-D array
• multiherm Returns the Hermitian parts of the matrices in a 3D array
• multiprod Matrix multiply 2-D slices of N-D arrays
• multiskewh Returns the skew-Hermitian parts of the matrices in the 3D matrix X.
• qr_unique Thin QR factorization ensuring diagonal of R is real, positive if possible.
This function is called by:

## SOURCE CODE

```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 %
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
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);
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```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005