Home > manopt > manifolds > stiefel > stiefelcomplexfactory.m

# stiefelcomplexfactory

## PURPOSE Returns a manifold struct. to optimize over complex orthonormal matrices.

## SYNOPSIS function M = stiefelcomplexfactory(n, p, k)

## DESCRIPTION ``` 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}
}

## CROSS-REFERENCE INFORMATION This function calls:
• 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
• qr_unique Thin QR factorization ensuring diagonal of R is real, positive if possible.
This function is called by:

## SOURCE CODE ```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 %
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.
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```

Generated on Sun 05-Sep-2021 17:57:00 by m2html © 2005