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 %   April 3, 2015 (NB):
0064 %       Replaced trace(A'*B) by A(:)'*B(:) : equivalent but faster.
0065 
0066 % TODO: modify normalize_rows and project_rows to work without transposes.
0067 % TODO: enhance ehess2rhess to also use bsxfun.
0068     
0069     
0070     if ~exist('lyap', 'file')
0071         warning('manopt:elliptopefactory:slowlyap', ...
0072                ['The function lyap to solve Lyapunov equations seems not to ' ...
0073                 'be available. This may slow down optimization over this ' ...
0074                 'manifold significantly. lyap is part of the control system ' ...
0075                 'toolbox.']);
0076     end
0077     
0078     
0079     M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with diagonal elements being 1', n, k);
0080     
0081     M.dim = @() n*(k-1) - k*(k-1)/2; % Extra -1 is because of the diagonal constraint that
0082     
0083     % Euclidean metric on the total space
0084     M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);
0085     
0086     M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));
0087     
0088     M.dist = @(Y, Z) error('elliptopefactory.dist not implemented yet.');
0089     
0090     M.typicaldist = @() 10*k;
0091     
0092     M.proj = @projection;
0093     
0094     M.tangent = M.proj;
0095     M.tangent2ambient = @(Y, eta) eta;
0096     
0097     M.retr = @retraction;
0098     
0099     M.egrad2rgrad = @egrad2rgrad;
0100     
0101     M.ehess2rhess = @ehess2rhess;
0102     
0103     M.exp = @exponential;
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     [unused, k] = size(Y); %#ok<ASGLU>
0154     eta = project_rows(Y, eta);
0155 
0156     % Projection onto the horizontal space
0157     YtY = Y'*Y;
0158     SS = YtY;
0159     AS = Y'*eta - eta'*Y;
0160     try
0161         % This is supposed to work and indeed return a skew-symmetric
0162         % solution Omega.
0163         Omega = lyap(SS, -AS);
0164     catch up %#ok<NASGU>
0165         % It can happen though that SS will be rank deficient. The
0166         % Lyapunov equation we solve still has a unique skew-symmetric
0167         % solution, but solutions with a symmetric part now also exist,
0168         % and the lyap function doesn't like that. So we want to
0169         % extract the minimum norm solution. This is also useful if lyap is
0170         % not available (it is part of the control system toolbox).
0171         mat = @(x) reshape(x, [k k]);
0172         vec = @(X) X(:);
0173         is_octave = exist('OCTAVE_VERSION', 'builtin');
0174         if ~is_octave
0175             [vecomega, unused] = minres(@(x) vec(SS*mat(x) + mat(x)*SS), vec(AS)); %#ok<NASGU>
0176         else
0177             [vecomega, unused] = gmres(@(x) vec(SS*mat(x) + mat(x)*SS), vec(AS)); %#ok<NASGU>
0178         end
0179         Omega = mat(vecomega);
0180     end
0181     % % Make sure the result is skew-symmetric (does not seem necessary).
0182     % Omega = (Omega-Omega')/2;
0183     etaproj = eta - Y*Omega;
0184 end
0185 
0186 % Retraction
0187 function Ynew = retraction(Y, eta, t)
0188     if nargin < 3
0189         t = 1.0;
0190     end
0191     Ynew = Y + t*eta;
0192     Ynew = normalize_rows(Ynew);
0193 end
0194 
0195 % Exponential map
0196 function Ynew = exponential(Y, eta, t)
0197     if nargin < 3
0198         t = 1.0;
0199     end
0200 
0201     Ynew = retraction(Y, eta, t);
0202     warning('manopt:elliptopefactory:exp', ...
0203         ['Exponential for fixed rank spectrahedron ' ...
0204         'manifold not implemented yet. Used retraction instead.\n' ...
0205         'To disable this warning: warning(''off'', ''manopt:elliptopefactory:exp'')']);
0206 end
0207 
0208 % Euclidean gradient to Riemannian gradient conversion.
0209 % We only need the ambient space projection: the remainder of the
0210 % projection function is not necessary because the Euclidean gradient must
0211 % already be orthogonal to the vertical space.
0212 function rgrad = egrad2rgrad(Y, egrad)
0213     rgrad = project_rows(Y, egrad);
0214 end
0215 
0216 % Euclidean Hessian to Riemannian Hessian conversion.
0217 % TODO: speed this function up using bsxfun.
0218 function Hess = ehess2rhess(Y, egrad, ehess, eta)
0219     k = size(Y, 2);
0220 
0221     % Directional derivative of the Riemannian gradient
0222     scaling_grad = sum((egrad.*Y), 2); % column vector of size n
0223     scaling_grad_repeat = scaling_grad*ones(1, k);
0224 
0225     Hess = ehess - scaling_grad_repeat.*eta;
0226 
0227     scaling_hess = sum((eta.*egrad) + (Y.*ehess), 2);
0228     scaling_hess_repeat = scaling_hess*ones(1, k);
0229     % directional derivative of scaling_grad_repeat
0230     Hess = Hess - scaling_hess_repeat.*Y;
0231 
0232     % Project on the horizontal space
0233     Hess = projection(Y, Hess);
0234 end
0235 
0236 % Random point generation on the manifold
0237 function Y = random(n, k)
0238     Y = randn(n, k);
0239     Y = normalize_rows(Y);
0240 end
0241 
0242 % Random vector generation at Y
0243 function eta = randomvec(Y)
0244     eta = randn(size(Y));
0245     eta = projection(Y, eta);
0246     nrm = norm(eta, 'fro');
0247     eta = eta / nrm;
0248 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005