Home > manopt > manifolds > symfixedrank > elliptopefactory.m

elliptopefactory

PURPOSE ^

Manifold of n-by-n psd matrices of rank k with unit diagonal elements.

SYNOPSIS ^

function M = elliptopefactory(n, k)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005