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".

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}
}

## CROSS-REFERENCE INFORMATION

This function calls:
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
This function is called by:
• maxcut Algorithm to (try to) compute a maximum cut of a graph, via SDP approach.

## 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".
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 %
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
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
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.
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
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