Home > manopt > manifolds > symfixedrank > spectrahedronfactory.m

# spectrahedronfactory

## PURPOSE Manifold of n-by-n symmetric positive semidefinite matrices of rank k

## SYNOPSIS function M = spectrahedronfactory(n, k)

## DESCRIPTION ``` Manifold of n-by-n symmetric positive semidefinite matrices of rank k
with trace (sum of diagonal elements) equal to 1.

function M = spectrahedronfactory(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 trace(Xdot) == 0.
The metric is the canonical Euclidean metric on Y.

The trace constraint on X (trace(X) == 1) translates to a unit Frobenius
norm constraint on Y: trace(X) = norm(Y, 'fro')^2 == 1. The set of such
Y's forms the unit sphere in R^(nxk): see spherefactory. 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.
As an alternative, you may use the sphere manifold (it has larger
dimension (by 1), 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:

## SOURCE CODE ```0001 function M = spectrahedronfactory(n, k)
0002 % Manifold of n-by-n symmetric positive semidefinite matrices of rank k
0003 % with trace (sum of diagonal elements) equal to 1.
0004 %
0005 % function M = spectrahedronfactory(n, k)
0006 %
0007 % A point X on the manifold is parameterized as YY^T where Y is a matrix of
0008 % size nxk. As such, X is symmetric, positive semidefinite. We restrict to
0009 % full-rank Y's, such that X has rank exactly k. The point X is numerically
0010 % represented by Y (this is more efficient than working with X, which may
0011 % be big). Tangent vectors are represented as matrices of the same size as
0012 % Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y and trace(Xdot) == 0.
0013 % The metric is the canonical Euclidean metric on Y.
0014 %
0015 % The trace constraint on X (trace(X) == 1) translates to a unit Frobenius
0016 % norm constraint on Y: trace(X) = norm(Y, 'fro')^2 == 1. The set of such
0017 % Y's forms the unit sphere in R^(nxk): see spherefactory. But because for
0018 % any orthogonal Q of size k, it holds that (YQ)(YQ)' = YY', we "group" all
0019 % matrices of the form YQ in an equivalence class. The set of equivalence
0020 % classes is a Riemannian quotient manifold, implemented here.
0021 %
0022 %
0023 % Note that this geometry formally breaks down at rank-deficient Y's.
0024 % As an alternative, you may use the sphere manifold (it has larger
0025 % dimension (by 1), 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 11, 2013.
0050 % Contributors: Nicolas Boumal
0051 % Change log:
0052 %
0053 %   Apr. 2, 2015 (NB):
0054 %       Replaced trace(A'*B) by A(:)'*B(:) (equivalent but faster).
0055 %       Updated documentation.
0056 %
0057 %   Apr. 17, 2018 (NB):
0058 %       Removed dependence on lyap.
0059 %
0060 %   Sep.  6, 2018 (NB):
0061 %       Removed M.exp() as it was not implemented.
0062
0063
0064
0065     M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with trace 1', n, k);
0066
0067     M.dim = @() n*k - 1 - k*(k-1)/2;
0068
0069     % Euclidean metric on the total space
0070     M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);
0071
0072     M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));
0073
0074     M.dist = @(Y, Z) error('spectrahedronfactory.dist not implemented yet.');
0075
0076     M.typicaldist = @() 10*k;
0077
0078     M.proj = @projection;
0079     function etaproj = projection(Y, eta)
0080         % Projection onto the tangent space, i.e., on the tangent space of
0081         % ||Y|| = 1
0082
0083         eta = eta - (eta(:)'*Y(:))*Y;
0084
0085         % Projection onto the horizontal space
0086         YtY = Y'*Y;
0087         SS = YtY;
0088         AS = Y'*eta - eta'*Y;
0089         % Omega = lyap(SS, -AS);
0090         Omega = lyapunov_symmetric(SS, AS);
0091         etaproj = eta - Y*Omega;
0092     end
0093
0094     M.tangent = M.proj;
0095     M.tangent2ambient = @(Y, eta) eta;
0096
0097     M.retr = @retraction;
0098     function Ynew = retraction(Y, eta, t)
0099         if nargin < 3
0100             t = 1.0;
0101         end
0102         Ynew = Y + t*eta;
0103         Ynew = Ynew/norm(Ynew, 'fro');
0104     end
0105
0106
0108
0109     M.ehess2rhess = @ehess2rhess;
0110     function Hess = ehess2rhess(Y, egrad, ehess, eta)
0111
0112         % Directional derivative of the Riemannian gradient
0113         Hess = ehess - (egrad(:)'*Y(:))*eta - ( (ehess(:)'*Y(:)) + (eta(:)'*egrad(:)) )*Y;
0114         Hess = Hess - (Hess(:)'*Y(:))*Y;
0115
0116         % Project on the horizontal space
0117         Hess = M.proj(Y, Hess);
0118
0119     end
0120
0121
0122     % Notice that the hash of two equivalent points will be different...
0123     M.hash = @(Y) ['z' hashmd5(Y(:))];
0124
0125     M.rand = @random;
0126
0127     function Y = random()
0128         Y = randn(n, k);
0129         Y = Y/norm(Y,'fro');
0130     end
0131
0132     M.randvec = @randomvec;
0133     function eta = randomvec(Y)
0134         eta = randn(n, k);
0135         eta = projection(Y, eta);
0136         nrm = M.norm(Y, eta);
0137         eta = eta / nrm;
0138     end
0139
0140     M.lincomb = @matrixlincomb;
0141
0142     M.zerovec = @(Y) zeros(n, k);
0143
0144     M.transp = @(Y1, Y2, d) projection(Y2, d);
0145
0146     M.vec = @(Y, u_mat) u_mat(:);
0147     M.mat = @(Y, u_vec) reshape(u_vec, [n, k]);
0148     M.vecmatareisometries = @() true;
0149
0150 end```

Generated on Sun 05-Sep-2021 17:57:00 by m2html © 2005