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:
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• multihconj MULTIHCONJ Hermitian conjugating arrays of matrices.
• multiherm Returns the Hermitian parts of the matrices in the 3D matrix X
• multiprod Multiplying 1-D or 2-D subarrays contained in two N-D arrays.
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
0050     if ~exist('k', 'var') || isempty(k)
0051         k = 1;
0052     end
0053
0054     if k == 1
0055         M.name = @() sprintf('Complex Stiefel manifold St(%d, %d)', n, p);
0056     elseif k > 1
0057         M.name = @() sprintf('Product complex Stiefel manifold St(%d, %d)^%d', n, p, k);
0058     else
0059         error('k must be an integer no less than 1.');
0060     end
0061
0062     M.dim = @() k*(2*n*p - p^2); %! k*(n*p - .5*p*(p+1)) -> k*(2*n*p - p^2)
0063
0064     M.inner = @(x, d1, d2) real(d1(:)'*d2(:)); %! trace -> real-trace
0065
0066     M.norm = @(x, d) norm(d(:));
0067
0068     M.dist = @(x, y) error('stiefel.dist not implemented yet.');
0069
0070     M.typicaldist = @() sqrt(p*k);
0071
0072     M.proj = @projection;
0073     function Up = projection(X, U)
0074
0075         XHU = multiprod(multihconj(X), U); %! XtU -> XHU, multitransp -> multihconj
0076         herXHU = multiherm(XHU); %! symXtU -> herXHU, multisym -> multiherm
0077         Up = U - multiprod(X, herXHU); %! symXtU -> herXHU
0078
0079     end
0080
0081     M.tangent = M.proj;
0082
0083     % For Riemannian submanifolds, converting a Euclidean gradient into a
0084     % Riemannian gradient amounts to an orthogonal projection.
0086
0087     M.ehess2rhess = @ehess2rhess;
0088     function rhess = ehess2rhess(X, egrad, ehess, H)
0089         XHG = multiprod(multihconj(X), egrad); %! XtG -> XHG, multitransp -> multihconj
0090         herXHG = multiherm(XHG); %! symXtG -> herXHG, multisym(XtG) -> multiherm(XHG)
0091         HherXHG = multiprod(H, herXHG); %! HsymXtG -> HherXHG, symXtG -> herXHG
0092         rhess = projection(X, ehess - HherXHG); %! HsymXtG -> HherXHG
0093     end
0094
0095     M.retr = @retraction;
0096     function Y = retraction(X, U, t)
0097         if nargin < 3
0098             t = 1.0;
0099         end
0100         Y = X + t*U;
0101         for i = 1 : k
0102             [Q, R] = qr(Y(:, :, i), 0);
0103             % The instruction with R assures we are not flipping signs
0104             % of some columns, which should never happen in modern Matlab
0105             % versions but may be an issue with older versions.
0106             Y(:, :, i) = Q * diag(sign(sign(diag(R))+.5));
0107         end
0108     end
0109
0110     M.exp = @exponential;
0111     function Y = exponential(X, U, t)
0112         if nargin == 2
0113             t = 1;
0114         end
0115         tU = t*U;
0116         Y = zeros(size(X));
0117         for i = 1 : k
0118             % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
0119             Xi = X(:, :, i);
0120             Ui = tU(:, :, i);
0121             Y(:, :, i) = [Xi Ui] * ...
0122                          expm([Xi'*Ui , -Ui'*Ui ; eye(p) , Xi'*Ui]) * ...
0123                          [ expm(-Xi'*Ui) ; zeros(p) ];
0124         end
0125
0126     end
0127
0128     M.hash = @(X) ['z' hashmd5([real(X(:)) ; imag(X(:))])]; %! X(:) -> [real(X(:)) ; imag(X(:))]
0129
0130     M.rand = @random;
0131     function X = random()
0132         X = zeros(n, p, k);
0133         for i = 1 : k
0134             [Q, unused] = qr(randn(n, p) + 1i*randn(n,p), 0); %#ok<NASGU> %! Complex version
0135             X(:, :, i) = Q;
0136         end
0137     end
0138
0139     M.randvec = @randomvec;
0140     function U = randomvec(X)
0141         U = projection(X, randn(n, p, k) + 1i*randn(n, p, k)); %! Complex version
0142         U = U / norm(U(:));
0143     end
0144
0145     M.lincomb = @matrixlincomb;
0146
0147     M.zerovec = @(x) zeros(n, p, k);
0148
0149     M.transp = @(x1, x2, d) projection(x2, d);
0150
0151     M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))];
0152     M.mat = @(x, u_vec) reshape(u_vec(1:(n*p*k)) + 1i*u_vec((n*p*k+1):end), [n, p, k]);
0153     M.vecmatareisometries = @() true; % TODO : to check.
0154
0155 end```

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005