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}
 }


 See also: stiefelfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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 %
0043 % See also: stiefelfactory
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.
0086     M.egrad2rgrad = M.proj;
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 Fri 30-Sep-2022 13:18:25 by m2html © 2005