Manifold of n-by-n symmetric positive semidefinite matrices of rank k
with trace (sum of diagonal elements) equal to 1.
function M = spectrahedronfactory(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 trace(Xdot) == 0.
The metric is the canonical Euclidean metric on Y.
The trace constraint on X (trace(X) == 1) translates to a unit Frobenius
norm constraint on Y: trace(X) = norm(Y, 'fro')^2 == 1. The set of such
Y's forms the unit sphere in R^(nxk): see spherefactory. 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.
As an alternative, you may use the sphere manifold (it has larger
dimension (by 1), 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: spherefactory elliptopefactory symfixedrankYYfactory0001 function M = spectrahedronfactory(n, k) 0002 % Manifold of n-by-n symmetric positive semidefinite matrices of rank k 0003 % with trace (sum of diagonal elements) equal to 1. 0004 % 0005 % function M = spectrahedronfactory(n, k) 0006 % 0007 % A point X on the manifold is parameterized as YY^T where Y is a matrix of 0008 % size nxk. As such, X is symmetric, positive semidefinite. We restrict to 0009 % full-rank Y's, such that X has rank exactly k. The point X is numerically 0010 % represented by Y (this is more efficient than working with X, which may 0011 % be big). Tangent vectors are represented as matrices of the same size as 0012 % Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y and trace(Xdot) == 0. 0013 % The metric is the canonical Euclidean metric on Y. 0014 % 0015 % The trace constraint on X (trace(X) == 1) translates to a unit Frobenius 0016 % norm constraint on Y: trace(X) = norm(Y, 'fro')^2 == 1. The set of such 0017 % Y's forms the unit sphere in R^(nxk): see spherefactory. But because for 0018 % any orthogonal Q of size k, it holds that (YQ)(YQ)' = YY', we "group" all 0019 % matrices of the form YQ in an equivalence class. The set of equivalence 0020 % classes is a Riemannian quotient manifold, implemented here. 0021 % 0022 % 0023 % Note that this geometry formally breaks down at rank-deficient Y's. 0024 % As an alternative, you may use the sphere manifold (it has larger 0025 % dimension (by 1), 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: spherefactory elliptopefactory symfixedrankYYfactory 0047 0048 % This file is part of Manopt: www.manopt.org. 0049 % Original author: Bamdev Mishra, July 11, 2013. 0050 % Contributors: Nicolas Boumal 0051 % Change log: 0052 % 0053 % Apr. 2, 2015 (NB): 0054 % Replaced trace(A'*B) by A(:)'*B(:) (equivalent but faster). 0055 % Updated documentation. 0056 % 0057 % Apr. 17, 2018 (NB): 0058 % Removed dependence on lyap. 0059 % 0060 % Sep. 6, 2018 (NB): 0061 % Removed M.exp() as it was not implemented. 0062 0063 0064 0065 M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with trace 1', n, k); 0066 0067 M.dim = @() n*k - 1 - k*(k-1)/2; 0068 0069 % Euclidean metric on the total space 0070 M.inner = @(Y, eta, zeta) eta(:)'*zeta(:); 0071 0072 M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta)); 0073 0074 M.dist = @(Y, Z) error('spectrahedronfactory.dist not implemented yet.'); 0075 0076 M.typicaldist = @() 10*k; 0077 0078 M.proj = @projection; 0079 function etaproj = projection(Y, eta) 0080 % Projection onto the tangent space, i.e., on the tangent space of 0081 % ||Y|| = 1 0082 0083 eta = eta - (eta(:)'*Y(:))*Y; 0084 0085 % Projection onto the horizontal space 0086 YtY = Y'*Y; 0087 SS = YtY; 0088 AS = Y'*eta - eta'*Y; 0089 % Omega = lyap(SS, -AS); 0090 Omega = lyapunov_symmetric(SS, AS); 0091 etaproj = eta - Y*Omega; 0092 end 0093 0094 M.tangent = M.proj; 0095 M.tangent2ambient = @(Y, eta) eta; 0096 0097 M.retr = @retraction; 0098 function Ynew = retraction(Y, eta, t) 0099 if nargin < 3 0100 t = 1.0; 0101 end 0102 Ynew = Y + t*eta; 0103 Ynew = Ynew/norm(Ynew, 'fro'); 0104 end 0105 0106 0107 M.egrad2rgrad = @(Y, eta) eta - (eta(:)'*Y(:))*Y; 0108 0109 M.ehess2rhess = @ehess2rhess; 0110 function Hess = ehess2rhess(Y, egrad, ehess, eta) 0111 0112 % Directional derivative of the Riemannian gradient 0113 Hess = ehess - (egrad(:)'*Y(:))*eta - ( (ehess(:)'*Y(:)) + (eta(:)'*egrad(:)) )*Y; 0114 Hess = Hess - (Hess(:)'*Y(:))*Y; 0115 0116 % Project on the horizontal space 0117 Hess = M.proj(Y, Hess); 0118 0119 end 0120 0121 0122 % Notice that the hash of two equivalent points will be different... 0123 M.hash = @(Y) ['z' hashmd5(Y(:))]; 0124 0125 M.rand = @random; 0126 0127 function Y = random() 0128 Y = randn(n, k); 0129 Y = Y/norm(Y,'fro'); 0130 end 0131 0132 M.randvec = @randomvec; 0133 function eta = randomvec(Y) 0134 eta = randn(n, k); 0135 eta = projection(Y, eta); 0136 nrm = M.norm(Y, eta); 0137 eta = eta / nrm; 0138 end 0139 0140 M.lincomb = @matrixlincomb; 0141 0142 M.zerovec = @(Y) zeros(n, k); 0143 0144 M.transp = @(Y1, Y2, d) projection(Y2, d); 0145 0146 M.vec = @(Y, u_mat) u_mat(:); 0147 M.mat = @(Y, u_vec) reshape(u_vec, [n, k]); 0148 M.vecmatareisometries = @() true; 0149 0150 end