function M = spectrahedronfactory(n, k)
% 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".
% 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: spherefactory elliptopefactory symfixedrankYYfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 11, 2013.
% Contributors: Nicolas Boumal
% Change log:
%
%   Apr. 2, 2015 (NB):
%       Replaced trace(A'*B) by A(:)'*B(:) (equivalent but faster).
%       Updated documentation.
%
%   Apr. 17, 2018 (NB):
%       Removed dependence on lyap.
%
%   Sep.  6, 2018 (NB):
%       Removed M.exp() as it was not implemented.
    
    
    
    M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with trace 1', n, k);
    
    M.dim = @() n*k - 1 - k*(k-1)/2;
    
    % Euclidean metric on the total space
    M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);
    
    M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));
    
    M.typicaldist = @() 10*k;
    
    M.proj = @projection;
    function etaproj = projection(Y, eta)
        % Projection onto the tangent space, i.e., on the tangent space of
        % ||Y|| = 1
        
        eta = eta - (eta(:)'*Y(:))*Y;
        
        % Projection onto the horizontal space
        YtY = Y'*Y;
        SS = YtY;
        AS = Y'*eta - eta'*Y;
        % Omega = lyap(SS, -AS);
        Omega = lyapunov_symmetric(SS, AS);
        etaproj = eta - Y*Omega;
    end
    
    M.tangent = M.proj;
    M.tangent2ambient = @(Y, eta) eta;
    
    M.retr = @retraction;
    function Ynew = retraction(Y, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Ynew = Y + t*eta;
        Ynew = Ynew/norm(Ynew, 'fro');
    end
    
    
    M.egrad2rgrad = @(Y, eta) eta - (eta(:)'*Y(:))*Y;
    
    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(Y, egrad, ehess, eta)
       
        % Directional derivative of the Riemannian gradient
        Hess = ehess - (egrad(:)'*Y(:))*eta - ( (ehess(:)'*Y(:)) + (eta(:)'*egrad(:)) )*Y;
        Hess = Hess - (Hess(:)'*Y(:))*Y;
        
        % Project on the horizontal space
        Hess = M.proj(Y, Hess);
        
    end
    
    
    % Notice that the hash of two equivalent points will be different...
    M.hash = @(Y) ['z' hashmd5(Y(:))];
    
    M.rand = @random;
    
    function Y = random()
        Y = randn(n, k);
        Y = Y/norm(Y,'fro');
    end
    
    M.randvec = @randomvec;
    function eta = randomvec(Y)
        eta = randn(n, k);
        eta = projection(Y, eta);
        nrm = M.norm(Y, eta);
        eta = eta / nrm;
    end
    
    M.lincomb = @matrixlincomb;
    
    M.zerovec = @(Y) zeros(n, k);
    
    M.transp = @(Y1, Y2, d) projection(Y2, d);
    
    M.vec = @(Y, u_mat) u_mat(:);
    M.mat = @(Y, u_vec) reshape(u_vec, [n, k]);
    M.vecmatareisometries = @() true;
    
end
