Manifold of n x n complex Hermitian pos. semidefinite matrices of rank k. function M = symfixedrankYYcomplexfactory(n, k) Manifold of n-by-n complex Hermitian positive semidefinite matrices of fixed rank k. This follows the quotient geometry described in Sarod Yatawatta's 2013 paper: "Radio interferometric calibration using a Riemannian manifold", ICASSP. Paper link: http://dx.doi.org/10.1109/ICASSP.2013.6638382. A point X on the manifold M is parameterized as YY^*, where Y is a complex matrix of size nxk. For any point Y on the manifold M, given any kxk complex unitary matrix U, we say Y*U is equivalent to Y, i.e., YY^* does not change. Therefore, M is the set of equivalence classes and is a Riemannian quotient manifold C^{nk}/SU(k). The metric is the usual real-trace inner product, that is, it is the usual metric for the complex plane identified with R^2. Notice that this manifold is not complete: if optimization leads Y to be rank-deficient, the geometry will break down. Hence, this geometry should only be used if it is expected that the points of interest will have rank exactly k. Reduce k if that is not the case. The geometry is based on the following papers (and references therein). Please cite the Manopt paper as well as the research papers: @INPROCEEDINGS{Yatawatta2013A, author={Yatawatta, S.}, booktitle={Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on}, title={Radio interferometric calibration using a {R}iemannian manifold}, year={2013}, month={May}, pages={3866--3870}, doi={10.1109/ICASSP.2013.6638382}, ISSN={1520-6149}, } @article{Yatawatta2013B, author = {Yatawatta, S.}, title = {On the interpolation of calibration solutions obtained in radio interferometry}, volume = {428}, number = {1}, pages = {828--833}, year = {2013}, doi = {10.1093/mnras/sts069}, journal = {Monthly Notices of the Royal Astronomical Society} } See also: symfixedrankYYfactory sympositivedefinitefactory
0001 function M = symfixedrankYYcomplexfactory(n, k) 0002 % Manifold of n x n complex Hermitian pos. semidefinite matrices of rank k. 0003 % 0004 % function M = symfixedrankYYcomplexfactory(n, k) 0005 % 0006 % Manifold of n-by-n complex Hermitian positive semidefinite matrices of 0007 % fixed rank k. This follows the quotient geometry described 0008 % in Sarod Yatawatta's 2013 paper: 0009 % "Radio interferometric calibration using a Riemannian manifold", ICASSP. 0010 % 0011 % Paper link: http://dx.doi.org/10.1109/ICASSP.2013.6638382. 0012 % 0013 % A point X on the manifold M is parameterized as YY^*, where 0014 % Y is a complex matrix of size nxk. For any point Y on the manifold M, 0015 % given any kxk complex unitary matrix U, we say Y*U is equivalent to Y, 0016 % i.e., YY^* does not change. Therefore, M is the set of equivalence 0017 % classes and is a Riemannian quotient manifold C^{nk}/SU(k). 0018 % The metric is the usual real-trace inner product, that is, 0019 % it is the usual metric for the complex plane identified with R^2. 0020 % 0021 % Notice that this manifold is not complete: if optimization leads Y to be 0022 % rank-deficient, the geometry will break down. Hence, this geometry should 0023 % only be used if it is expected that the points of interest will have rank 0024 % exactly k. Reduce k if that is not the case. 0025 % 0026 % The geometry is based on the following papers (and references therein). 0027 % Please cite the Manopt paper as well as the research papers: 0028 % 0029 % @INPROCEEDINGS{Yatawatta2013A, 0030 % author={Yatawatta, S.}, 0031 % booktitle={Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on}, 0032 % title={Radio interferometric calibration using a {R}iemannian manifold}, 0033 % year={2013}, 0034 % month={May}, 0035 % pages={3866--3870}, 0036 % doi={10.1109/ICASSP.2013.6638382}, 0037 % ISSN={1520-6149}, 0038 % } 0039 % 0040 % @article{Yatawatta2013B, 0041 % author = {Yatawatta, S.}, 0042 % title = {On the interpolation of calibration solutions obtained in radio interferometry}, 0043 % volume = {428}, 0044 % number = {1}, 0045 % pages = {828--833}, 0046 % year = {2013}, 0047 % doi = {10.1093/mnras/sts069}, 0048 % journal = {Monthly Notices of the Royal Astronomical Society} 0049 % } 0050 % 0051 % See also: symfixedrankYYfactory sympositivedefinitefactory 0052 0053 0054 % This file is part of Manopt: www.manopt.org. 0055 % Original author: Sarod Yatawatta, June 29, 2015. 0056 % Contributors: Bamdev Mishra. 0057 % Change log: 0058 % 0059 % June 28, 2016 (NB): 0060 % Metric scaled down by factor 2 to match the metric used in 0061 % euclideancomplexfactory. 0062 % 0063 % Apr. 18, 2018 (NB): 0064 % Removed lyap depedency. 0065 % 0066 % Sep. 6, 2018 (NB): 0067 % Removed M.exp() as it was not implemented. 0068 0069 M.name = @() sprintf('YY'' quotient manifold of Hermitian %dx%d complex matrices of rank %d.', n, n, k); 0070 0071 M.dim = @() 2*k*n - k*k; % SY: dim of ambient space (2*k*n) - dim of kxk unitary matrix (k^2). 0072 0073 % Euclidean metric on the total space. 0074 % BM: equivalent to real(trace(eta'*zeta)), but more efficient. 0075 M.inner = @(Y, eta, zeta) real(eta(:)'*zeta(:)); 0076 0077 M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta)); 0078 0079 % Find unitary U to minimize ||Y - Z*U||, 0080 % i.e., the Procrustes problem, with svd(Y'*Z). 0081 M.dist = @(Y, Z) distance; 0082 function distval = distance(Y, Z) 0083 [u, ignore, v] = svd(Z'*Y); %#ok<ASGLU> 0084 E = Y - Z*u*v'; % SY: checked. 0085 distval = real(E(:)'*E(:)); 0086 end 0087 0088 M.typicaldist = @() 10*k; % BM: To do. 0089 0090 M.proj = @projection; 0091 function etaproj = projection(Y, eta) 0092 % Projection onto the horizontal space 0093 xx = Y'*Y; 0094 rr = Y'*eta - eta'*Y; 0095 Omega = lyapunov_symmetric(xx, rr); 0096 etaproj = eta - Y*Omega; 0097 end 0098 0099 M.tangent = M.proj; 0100 M.tangent2ambient = @(Y, eta) eta; 0101 0102 M.retr = @retraction; 0103 function Ynew = retraction(Y, eta, t) 0104 if nargin < 3 0105 t = 1.0; 0106 end 0107 Ynew = Y + t*eta; 0108 end 0109 0110 0111 M.egrad2rgrad = @(Y, eta) eta; 0112 M.ehess2rhess = @(Y, egrad, ehess, U) M.proj(Y, ehess); 0113 0114 0115 % Notice that the hash of two equivalent points will be different... 0116 M.hash = @(Y) ['z' hashmd5([real(Y(:)); imag(Y(:))])]; 0117 0118 M.rand = @random; 0119 function Y = random() 0120 Y = randn(n, k) + 1i*randn(n,k); 0121 end 0122 0123 M.randvec = @randomvec; 0124 function eta = randomvec(Y) 0125 eta = randn(n, k) + 1i*randn(n,k); 0126 eta = projection(Y, eta); 0127 nrm = M.norm(Y, eta); 0128 eta = eta / nrm; 0129 end 0130 0131 M.lincomb = @matrixlincomb; 0132 0133 M.zerovec = @(Y) zeros(n, k); 0134 0135 M.transp = @(Y1, Y2, d) projection(Y2, d); 0136 0137 M.vec = @(Y, u_mat) [real(u_mat(:)); imag(u_mat(:))]; 0138 M.mat = @(Y, u_vec) reshape(u_vec(1 : n*k), [n, k]) + 1i*reshape(u_vec(n*k + 1: end), [n, k]); 0139 M.vecmatareisometries = @() true; 0140 0141 end