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:
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• hashmd5 Computes the MD5 hash of input data.
• lyapunov_symmetric Solves AX + XA = C when A = A', as a pseudo-inverse.
• 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 %   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
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
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.
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
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