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     if k < 2
0079         warning('manopt:elliptopefactory:lowk', ...
0080                 'k should be an integer >= 2. At k = 1, the set is discrete.');
0081     end
0082     
0083     
0084     M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with diagonal elements being 1', n, k);
0085     
0086     M.dim = @() n*(k-1) - k*(k-1)/2; % Extra -1 is because of the diagonal constraint that
0087     
0088     % Euclidean metric on the total space
0089     M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);
0090     
0091     M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));
0092     
0093     M.dist = @(Y, Z) error('elliptopefactory.dist not implemented yet.');
0094     
0095     M.typicaldist = @() 10*k;
0096     
0097     M.proj = @projection;
0098     
0099     M.tangent = M.proj;
0100     M.tangent2ambient = @(Y, eta) eta;
0101     
0102     M.retr = @retraction;
0103     
0104     M.egrad2rgrad = @egrad2rgrad;
0105     
0106     M.ehess2rhess = @ehess2rhess;
0107     
0108     M.exp = @exponential;
0109     
0110     % Notice that the hash of two equivalent points will be different...
0111     M.hash = @(Y) ['z' hashmd5(Y(:))];
0112     
0113     M.rand = @() random(n, k);
0114     
0115     M.randvec = @randomvec;
0116     
0117     M.lincomb = @matrixlincomb;
0118     
0119     M.zerovec = @(Y) zeros(n, k);
0120     
0121     M.transp = @(Y1, Y2, d) projection(Y2, d);
0122     
0123     M.vec = @(Y, u_mat) u_mat(:);
0124     M.mat = @(Y, u_vec) reshape(u_vec, [n, k]);
0125     M.vecmatareisometries = @() true;
0126     
0127 end
0128 
0129 % Given a matrix X, returns the same matrix but with each column scaled so
0130 % that they have unit 2-norm.
0131 % See obliquefactory.
0132 function X = normalize_rows(X)
0133     X = X';
0134     norms = sqrt(sum(X.^2, 1));
0135     X = bsxfun(@times, X, 1./norms);
0136     X = X';
0137 end
0138 
0139 % Orthogonal projection of each row of H to the tangent space at the
0140 % corresponding row of X, seen as a point on a sphere.
0141 % See obliquefactory.
0142 function PXH = project_rows(X, H)
0143     X = X';
0144     H = H';
0145     % Compute the inner product between each vector H(:, i) with its root
0146     % point X(:, i), that is, X(:, i).' * H(:, i). Returns a row vector.
0147     inners = sum(X.*H, 1);
0148     % Subtract from H the components of the H(:, i)'s that are parallel to
0149     % the root points X(:, i).
0150     PXH = H - bsxfun(@times, X, inners);
0151     PXH = PXH';
0152 end
0153 
0154 
0155 % Projection onto the tangent space, i.e., on the tangent space of
0156 % ||Y(i, :)|| = 1
0157 function etaproj = projection(Y, eta)
0158     [unused, k] = size(Y); %#ok<ASGLU>
0159     eta = project_rows(Y, eta);
0160 
0161     % Projection onto the horizontal space
0162     YtY = Y'*Y;
0163     SS = YtY;
0164     AS = Y'*eta - eta'*Y;
0165     try
0166         % This is supposed to work and indeed return a skew-symmetric
0167         % solution Omega.
0168         Omega = lyap(SS, -AS);
0169     catch up %#ok<NASGU>
0170         % It can happen though that SS will be rank deficient. The
0171         % Lyapunov equation we solve still has a unique skew-symmetric
0172         % solution, but solutions with a symmetric part now also exist,
0173         % and the lyap function doesn't like that. So we want to
0174         % extract the minimum norm solution. This is also useful if lyap is
0175         % not available (it is part of the control system toolbox).
0176         mat = @(x) reshape(x, [k k]);
0177         vec = @(X) X(:);
0178         is_octave = exist('OCTAVE_VERSION', 'builtin');
0179         if ~is_octave
0180             [vecomega, unused] = minres(@(x) vec(SS*mat(x) + mat(x)*SS), vec(AS)); %#ok<NASGU>
0181         else
0182             [vecomega, unused] = gmres(@(x) vec(SS*mat(x) + mat(x)*SS), vec(AS)); %#ok<NASGU>
0183         end
0184         Omega = mat(vecomega);
0185     end
0186     % % Make sure the result is skew-symmetric (does not seem necessary).
0187     % Omega = (Omega-Omega')/2;
0188     etaproj = eta - Y*Omega;
0189 end
0190 
0191 % Retraction
0192 function Ynew = retraction(Y, eta, t)
0193     if nargin < 3
0194         t = 1.0;
0195     end
0196     Ynew = Y + t*eta;
0197     Ynew = normalize_rows(Ynew);
0198 end
0199 
0200 % Exponential map
0201 function Ynew = exponential(Y, eta, t)
0202     if nargin < 3
0203         t = 1.0;
0204     end
0205 
0206     Ynew = retraction(Y, eta, t);
0207     warning('manopt:elliptopefactory:exp', ...
0208         ['Exponential for fixed rank spectrahedron ' ...
0209         'manifold not implemented yet. Used retraction instead.\n' ...
0210         'To disable this warning: warning(''off'', ''manopt:elliptopefactory:exp'')']);
0211 end
0212 
0213 % Euclidean gradient to Riemannian gradient conversion.
0214 % We only need the ambient space projection: the remainder of the
0215 % projection function is not necessary because the Euclidean gradient must
0216 % already be orthogonal to the vertical space.
0217 function rgrad = egrad2rgrad(Y, egrad)
0218     rgrad = project_rows(Y, egrad);
0219 end
0220 
0221 % Euclidean Hessian to Riemannian Hessian conversion.
0222 % TODO: speed this function up using bsxfun.
0223 function Hess = ehess2rhess(Y, egrad, ehess, eta)
0224     k = size(Y, 2);
0225 
0226     % Directional derivative of the Riemannian gradient
0227     scaling_grad = sum((egrad.*Y), 2); % column vector of size n
0228     scaling_grad_repeat = scaling_grad*ones(1, k);
0229 
0230     Hess = ehess - scaling_grad_repeat.*eta;
0231 
0232     scaling_hess = sum((eta.*egrad) + (Y.*ehess), 2);
0233     scaling_hess_repeat = scaling_hess*ones(1, k);
0234     % directional derivative of scaling_grad_repeat
0235     Hess = Hess - scaling_hess_repeat.*Y;
0236 
0237     % Project on the horizontal space
0238     Hess = projection(Y, Hess);
0239 end
0240 
0241 % Random point generation on the manifold
0242 function Y = random(n, k)
0243     Y = randn(n, k);
0244     Y = normalize_rows(Y);
0245 end
0246 
0247 % Random vector generation at Y
0248 function eta = randomvec(Y)
0249     eta = randn(size(Y));
0250     eta = projection(Y, eta);
0251     nrm = norm(eta, 'fro');
0252     eta = eta / nrm;
0253 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005