Home > manopt > manifolds > symfixedrank > spectrahedronfactory.m

spectrahedronfactory

PURPOSE ^

Manifold of n-by-n symmetric positive semidefinite matrices of rank k

SYNOPSIS ^

function M = spectrahedronfactory(n, k)

DESCRIPTION ^

 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 symfixedrankYYfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 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

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005