Manifold of n-by-n symmetric positive semidefinite matrices of rank k.
function M = symfixedrankYYfactory(n, k)
A point X on the manifold is parameterized as YY^T where Y is a matrix of
size nxk. As such, X is symmetric, positive semidefinite. We restrict to
full-rank Y's, such that X has rank exactly k. The point X is numerically
represented by Y (this is more efficient than working with X, which may
be big). Tangent vectors are represented as matrices of the same size as
Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y'. The metric is the
canonical Euclidean metric on Y.
Since for any orthogonal Q of size k, it holds that (YQ)(YQ)' = YY',
we "group" all matrices of the form YQ in an equivalence class. The set
of equivalence classes is a Riemannian quotient manifold, implemented
here.
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.
An alternative, complete, geometry for positive semidefinite matrices of
rank k is described in Bonnabel and Sepulchre 2009, "Riemannian Metric
and Geometric Mean for Positive Semidefinite Matrices of Fixed Rank",
SIAM Journal on Matrix Analysis and Applications.
The geometry here implemented is the simplest case of the 2010 paper:
M. Journee, P.-A. Absil, F. Bach and R. Sepulchre,
"Low-Rank Optimization on the Cone of Positive Semidefinite Matrices".
Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf
The expressions for the distance and the logarithm come from the 2018
preprint:
Estelle Massart, P.-A. Absil,
"Quotient geometry with simple geodesics for the manifold of fixed-rank
positive-semidefinite matrices".
Paper link: https://sites.uclouvain.be/absil/2018-06/quotient_tech_report.pdf
Please cite the Manopt paper as well as the research papers:
@Article{journee2010low,
Title = {Low-rank optimization on the cone of positive semidefinite matrices},
Author = {Journ{\'e}e, M. and Bach, F. and Absil, P.-A. and Sepulchre, R.},
Journal = {SIAM Journal on Optimization},
Year = {2010},
Number = {5},
Pages = {2327--2351},
Volume = {20},
Doi = {10.1137/080731359}
}
@TechReport{massart2018quotient,
author = {Massart, Estelle and Absil, P.-A.},
title = {Quotient geometry with simple geodesics for the manifold of fixed-rank positive-semidefinite matrices},
institution = {UCLouvain},
year = {2018},
number = {UCL-INMA-2018.06-v2},
note = {Preprint: \url{http://sites.uclouvain.be/absil/2018.06}},
}0001 function M = symfixedrankYYfactory(n, k) 0002 % Manifold of n-by-n symmetric positive semidefinite matrices of rank k. 0003 % 0004 % function M = symfixedrankYYfactory(n, k) 0005 % 0006 % A point X on the manifold is parameterized as YY^T where Y is a matrix of 0007 % size nxk. As such, X is symmetric, positive semidefinite. We restrict to 0008 % full-rank Y's, such that X has rank exactly k. The point X is numerically 0009 % represented by Y (this is more efficient than working with X, which may 0010 % be big). Tangent vectors are represented as matrices of the same size as 0011 % Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y'. The metric is the 0012 % canonical Euclidean metric on Y. 0013 % 0014 % Since for any orthogonal Q of size k, it holds that (YQ)(YQ)' = YY', 0015 % we "group" all matrices of the form YQ in an equivalence class. The set 0016 % of equivalence classes is a Riemannian quotient manifold, implemented 0017 % here. 0018 % 0019 % Notice that this manifold is not complete: if optimization leads Y to be 0020 % rank deficient, the geometry will break down. Hence, this geometry should 0021 % only be used if it is expected that the points of interest will have rank 0022 % exactly k. Reduce k if that is not the case. 0023 % 0024 % An alternative, complete, geometry for positive semidefinite matrices of 0025 % rank k is described in Bonnabel and Sepulchre 2009, "Riemannian Metric 0026 % and Geometric Mean for Positive Semidefinite Matrices of Fixed Rank", 0027 % SIAM Journal on Matrix Analysis and Applications. 0028 % 0029 % 0030 % The geometry here implemented is the simplest case of the 2010 paper: 0031 % M. Journee, P.-A. Absil, F. Bach and R. Sepulchre, 0032 % "Low-Rank Optimization on the Cone of Positive Semidefinite Matrices". 0033 % Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf 0034 % The expressions for the distance and the logarithm come from the 2018 0035 % preprint: 0036 % Estelle Massart, P.-A. Absil, 0037 % "Quotient geometry with simple geodesics for the manifold of fixed-rank 0038 % positive-semidefinite matrices". 0039 % Paper link: https://sites.uclouvain.be/absil/2018-06/quotient_tech_report.pdf 0040 % 0041 % 0042 % Please cite the Manopt paper as well as the research papers: 0043 % @Article{journee2010low, 0044 % Title = {Low-rank optimization on the cone of positive semidefinite matrices}, 0045 % Author = {Journ{\'e}e, M. and Bach, F. and Absil, P.-A. and Sepulchre, R.}, 0046 % Journal = {SIAM Journal on Optimization}, 0047 % Year = {2010}, 0048 % Number = {5}, 0049 % Pages = {2327--2351}, 0050 % Volume = {20}, 0051 % Doi = {10.1137/080731359} 0052 % } 0053 % 0054 % @TechReport{massart2018quotient, 0055 % author = {Massart, Estelle and Absil, P.-A.}, 0056 % title = {Quotient geometry with simple geodesics for the manifold of fixed-rank positive-semidefinite matrices}, 0057 % institution = {UCLouvain}, 0058 % year = {2018}, 0059 % number = {UCL-INMA-2018.06-v2}, 0060 % note = {Preprint: \url{http://sites.uclouvain.be/absil/2018.06}}, 0061 % } 0062 % 0063 0064 % See also: elliptopefactory spectrahedronfactory symfixedrankYYcomplexfactory 0065 0066 % This file is part of Manopt: www.manopt.org. 0067 % Original author: Bamdev Mishra, Dec. 30, 2012. 0068 % Contributors: Estelle Massart 0069 % Change log: 0070 % 0071 % July 10, 2013 (NB): 0072 % Added vec, mat, tangent, tangent2ambient ; 0073 % Correction for the dimension of the manifold. 0074 % 0075 % Apr. 2, 2015 (NB): 0076 % Replaced trace(A'*B) by A(:)'*B(:) (equivalent but faster). 0077 % 0078 % Apr. 17, 2018 (NB): 0079 % Removed dependence on lyap. 0080 % 0081 % Sep. 6, 2018 (NB): 0082 % Removed M.exp() as it was not implemented. 0083 % 0084 % June 7, 2019 (EM): 0085 % Added M.dist, M.exp, M.log and M.invretr. 0086 0087 M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d', n, k); 0088 0089 M.dim = @() k*n - k*(k-1)/2; 0090 0091 % Euclidean metric on the total space 0092 M.inner = @(Y, eta, zeta) eta(:)'*zeta(:); 0093 0094 M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta)); 0095 0096 M.dist = @(Y, Z) norm(logarithm(Y,Z),'fro'); 0097 0098 M.typicaldist = @() 10*k; 0099 0100 M.proj = @projection; 0101 function etaproj = projection(Y, eta) 0102 % Projection onto the horizontal space 0103 YtY = Y'*Y; 0104 SS = YtY; 0105 AS = Y'*eta - eta'*Y; 0106 % Omega = lyap(SS, -AS); 0107 Omega = lyapunov_symmetric(SS, AS); 0108 etaproj = eta - Y*Omega; 0109 end 0110 0111 M.tangent = M.proj; 0112 M.tangent2ambient = @(Y, eta) eta; 0113 0114 M.exp = @exponential; 0115 function Ynew = exponential(Y, eta, t) 0116 if nargin < 3 0117 t = 1.0; 0118 end 0119 Ynew = Y + t*eta; 0120 end 0121 0122 M.retr = M.exp; 0123 0124 M.log = @logarithm; 0125 function eta = logarithm(Y, Z) 0126 YtZ = Y'*Z; 0127 [U, ~, V] = svd(YtZ); 0128 Qt = V*U'; 0129 eta = Z*Qt - Y; 0130 end 0131 0132 M.invretr = M.log; 0133 0134 M.egrad2rgrad = @(Y, eta) eta; 0135 M.ehess2rhess = @(Y, egrad, ehess, U) M.proj(Y, ehess); 0136 0137 % Notice that the hash of two equivalent points will be different... 0138 M.hash = @(Y) ['z' hashmd5(Y(:))]; 0139 0140 M.rand = @random; 0141 function Y = random() 0142 Y = randn(n, k); 0143 end 0144 0145 M.randvec = @randomvec; 0146 function eta = randomvec(Y) 0147 eta = randn(n, k); 0148 eta = projection(Y, eta); 0149 nrm = M.norm(Y, eta); 0150 eta = eta / nrm; 0151 end 0152 0153 M.lincomb = @matrixlincomb; 0154 0155 M.zerovec = @(Y) zeros(n, k); 0156 0157 M.transp = @(Y1, Y2, d) projection(Y2, d); 0158 0159 M.vec = @(Y, u_mat) u_mat(:); 0160 M.mat = @(Y, u_vec) reshape(u_vec, [n, k]); 0161 M.vecmatareisometries = @() true; 0162 0163 end