Manifold of n-by-n psd matrices of rank k with unit diagonal elements. function M = elliptopefactory(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 and diag(Xdot) == 0. The metric is the canonical Euclidean metric on Y. The diagonal constraints on X (X(i, i) == 1 for all i) translate to unit-norm constraints on the rows of Y: norm(Y(i, :)) == 1 for all i. The set of such Y's forms the oblique manifold. But because 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. Note that this geometry formally breaks down at rank-deficient Y's. This does not appear to be a major issue in practice when optimization algorithms converge to rank-deficient Y's, but convergence theorems no longer hold. As an alternative, you may use the oblique manifold (it has larger dimension, but does not break down at rank drop.) The geometry is taken from 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 Please cite the Manopt paper as well as the research paper: @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} } See also: obliquefactory symfixedrankYYfactory spectrahedronfactory
0001 function M = elliptopefactory(n, k) 0002 % Manifold of n-by-n psd matrices of rank k with unit diagonal elements. 0003 % 0004 % function M = elliptopefactory(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 and diag(Xdot) == 0. 0012 % The metric is the canonical Euclidean metric on Y. 0013 % 0014 % The diagonal constraints on X (X(i, i) == 1 for all i) translate to 0015 % unit-norm constraints on the rows of Y: norm(Y(i, :)) == 1 for all i. 0016 % The set of such Y's forms the oblique manifold. But because for any 0017 % orthogonal Q of size k it holds that (YQ)(YQ)' = YY', we "group" all 0018 % matrices of the form YQ in an equivalence class. The set of equivalence 0019 % classes is a Riemannian quotient manifold, implemented here. 0020 % 0021 % Note that this geometry formally breaks down at rank-deficient Y's. 0022 % This does not appear to be a major issue in practice when optimization 0023 % algorithms converge to rank-deficient Y's, but convergence theorems no 0024 % longer hold. As an alternative, you may use the oblique manifold (it has 0025 % larger dimension, but does not break down at rank drop.) 0026 % 0027 % The geometry is taken from the 2010 paper: 0028 % M. Journee, P.-A. Absil, F. Bach and R. Sepulchre, 0029 % "Low-Rank Optimization on the Cone of Positive Semidefinite Matrices". 0030 % Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf 0031 % 0032 % 0033 % Please cite the Manopt paper as well as the research paper: 0034 % @Article{journee2010low, 0035 % Title = {Low-rank optimization on the cone of positive semidefinite matrices}, 0036 % Author = {Journ{\'e}e, M. and Bach, F. and Absil, P.-A. and Sepulchre, R.}, 0037 % Journal = {SIAM Journal on Optimization}, 0038 % Year = {2010}, 0039 % Number = {5}, 0040 % Pages = {2327--2351}, 0041 % Volume = {20}, 0042 % Doi = {10.1137/080731359} 0043 % } 0044 % 0045 % 0046 % See also: obliquefactory symfixedrankYYfactory spectrahedronfactory 0047 0048 % This file is part of Manopt: www.manopt.org. 0049 % Original author: Bamdev Mishra, July 12, 2013. 0050 % Contributors: 0051 % Change log: 0052 % July 18, 2013 (NB): 0053 % Fixed projection operator for rank-deficient Y'Y. 0054 % 0055 % Aug. 8, 2013 (NB): 0056 % No longer using nested functions, to aim at Octave compatibility. 0057 % Sign error in right hand side of the call to minres corrected. 0058 % 0059 % June 24, 2014 (NB): 0060 % Used code snippets from obliquefactory to speed up projection, 0061 % retraction, egrad2rgrad and rand: the code now uses bsxfun for this. 0062 % 0063 % Apr. 3, 2015 (NB): 0064 % Replaced trace(A'*B) by A(:)'*B(:) : equivalent but faster. 0065 % 0066 % Apr. 17, 2018 (NB): 0067 % Removed dependency on lyap entirely. 0068 % 0069 % Sep. 6, 2018 (NB): 0070 % Removed M.exp() as it was not implemented. 0071 0072 % TODO: modify normalize_rows and project_rows to work without transposes. 0073 % TODO: enhance ehess2rhess to also use bsxfun. 0074 0075 if k < 2 0076 warning('manopt:elliptopefactory:lowk', ... 0077 'k should be an integer >= 2. At k = 1, the set is discrete.'); 0078 end 0079 0080 0081 M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with diagonal elements being 1', n, k); 0082 0083 M.dim = @() n*(k-1) - k*(k-1)/2; % Extra -1 is because of the diagonal constraint that 0084 0085 % Euclidean metric on the total space 0086 M.inner = @(Y, eta, zeta) eta(:)'*zeta(:); 0087 0088 M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta)); 0089 0090 M.dist = @(Y, Z) error('elliptopefactory.dist not implemented yet.'); 0091 0092 M.typicaldist = @() 10*k; 0093 0094 M.proj = @projection; 0095 0096 M.tangent = M.proj; 0097 M.tangent2ambient = @(Y, eta) eta; 0098 0099 M.retr = @retraction; 0100 0101 M.egrad2rgrad = @egrad2rgrad; 0102 0103 M.ehess2rhess = @ehess2rhess; 0104 0105 % Notice that the hash of two equivalent points will be different... 0106 M.hash = @(Y) ['z' hashmd5(Y(:))]; 0107 0108 M.rand = @() random(n, k); 0109 0110 M.randvec = @randomvec; 0111 0112 M.lincomb = @matrixlincomb; 0113 0114 M.zerovec = @(Y) zeros(n, k); 0115 0116 M.transp = @(Y1, Y2, d) projection(Y2, d); 0117 0118 M.vec = @(Y, u_mat) u_mat(:); 0119 M.mat = @(Y, u_vec) reshape(u_vec, [n, k]); 0120 M.vecmatareisometries = @() true; 0121 0122 end 0123 0124 % Given a matrix X, returns the same matrix but with each column scaled so 0125 % that they have unit 2-norm. 0126 % See obliquefactory. 0127 function X = normalize_rows(X) 0128 X = X'; 0129 norms = sqrt(sum(X.^2, 1)); 0130 X = bsxfun(@times, X, 1./norms); 0131 X = X'; 0132 end 0133 0134 % Orthogonal projection of each row of H to the tangent space at the 0135 % corresponding row of X, seen as a point on a sphere. 0136 % See obliquefactory. 0137 function PXH = project_rows(X, H) 0138 X = X'; 0139 H = H'; 0140 % Compute the inner product between each vector H(:, i) with its root 0141 % point X(:, i), that is, X(:, i).' * H(:, i). Returns a row vector. 0142 inners = sum(X.*H, 1); 0143 % Subtract from H the components of the H(:, i)'s that are parallel to 0144 % the root points X(:, i). 0145 PXH = H - bsxfun(@times, X, inners); 0146 PXH = PXH'; 0147 end 0148 0149 0150 % Projection onto the tangent space, i.e., on the tangent space of 0151 % ||Y(i, :)|| = 1 0152 function etaproj = projection(Y, eta) 0153 eta = project_rows(Y, eta); 0154 0155 % Projection onto the horizontal space 0156 YtY = Y'*Y; 0157 SS = YtY; 0158 AS = Y'*eta - eta'*Y; 0159 Omega = lyapunov_symmetric(SS, AS); 0160 0161 % It does not seem necessary to enforce skew-symmetry numerically. 0162 % Omega = (Omega-Omega')/2; 0163 0164 etaproj = eta - Y*Omega; 0165 end 0166 0167 % Retraction 0168 function Ynew = retraction(Y, eta, t) 0169 if nargin < 3 0170 t = 1.0; 0171 end 0172 Ynew = Y + t*eta; 0173 Ynew = normalize_rows(Ynew); 0174 end 0175 0176 % Euclidean gradient to Riemannian gradient conversion. 0177 % We only need the ambient space projection: the remainder of the 0178 % projection function is not necessary because the Euclidean gradient must 0179 % already be orthogonal to the vertical space. 0180 function rgrad = egrad2rgrad(Y, egrad) 0181 rgrad = project_rows(Y, egrad); 0182 end 0183 0184 % Euclidean Hessian to Riemannian Hessian conversion. 0185 % TODO: speed this function up using bsxfun. 0186 function Hess = ehess2rhess(Y, egrad, ehess, eta) 0187 k = size(Y, 2); 0188 0189 % Directional derivative of the Riemannian gradient 0190 scaling_grad = sum((egrad.*Y), 2); % column vector of size n 0191 scaling_grad_repeat = scaling_grad*ones(1, k); 0192 0193 Hess = ehess - scaling_grad_repeat.*eta; 0194 0195 scaling_hess = sum((eta.*egrad) + (Y.*ehess), 2); 0196 scaling_hess_repeat = scaling_hess*ones(1, k); 0197 % directional derivative of scaling_grad_repeat 0198 Hess = Hess - scaling_hess_repeat.*Y; 0199 0200 % Project on the horizontal space 0201 Hess = projection(Y, Hess); 0202 end 0203 0204 % Random point generation on the manifold 0205 function Y = random(n, k) 0206 Y = randn(n, k); 0207 Y = normalize_rows(Y); 0208 end 0209 0210 % Random vector generation at Y 0211 function eta = randomvec(Y) 0212 eta = randn(size(Y)); 0213 eta = projection(Y, eta); 0214 nrm = norm(eta, 'fro'); 0215 eta = eta / nrm; 0216 end