function M = SymPosSemDefinMatFixRank(n, k)
% Manifold of n-by-n symmetric positive semidefinite matrices of rank k.
%
% function M = SymPosSemDefinMatFixRank(n, k)
%
% A point X on the manifold is parameterized as YY^T = VDV^T where Y is a 
% matrix of size nxk, V is an orthogonal matrix of size nxk, and D is a 
% diagonal matrix of kxk. 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 V and D (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'. The metric is the
% canonical Euclidean metric on Y.
% 
% Since 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.
%
% Notice that this manifold is not complete: if optimization leads Y to be
% rank deficient, the geometry will break down. Hence, this geometry should
% only be used if it is expected that the points of interest will have rank
% exactly k. Reduce k if that is not the case.
% 
% An alternative, complete, geometry for positive semidefinite matrices of
% rank k is described in Bonnabel and Sepulchre 2009, "Riemannian Metric
% and Geometric Mean for Positive Semidefinite Matrices of Fixed Rank",
% SIAM Journal on Matrix Analysis and Applications.
%
%
% The geometry here implemented is the simplest case of 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
% The expressions for the distance and the logarithm come from the 2018
% preprint:
% Estelle Massart, P.-A. Absil, 
% "Quotient geometry with simple geodesics for the manifold of fixed-rank
% positive-semidefinite matrices".
% Paper link: https://sites.uclouvain.be/absil/2018-06/quotient_tech_report.pdf
% 
% 
% Please cite the Manopt paper as well as the research papers:
%     @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}
%     }
%
%     @TechReport{massart2018quotient,
%       author        = {Massart, Estelle and Absil, P.-A.},
%       title         = {Quotient geometry with simple geodesics for the manifold of fixed-rank positive-semidefinite matrices},
%       institution   = {UCLouvain},
%       year          = {2018},
%       number        = {UCL-INMA-2018.06-v2},
%       note          = {Preprint: \url{http://sites.uclouvain.be/absil/2018.06}},
%     }
% 

% See also: elliptopefactory spectrahedronfactory symfixedrankYYcomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors: Estelle Massart
% Change log:
%
%   July 10, 2013 (NB):
%       Added vec, mat, tangent, tangent2ambient ;
%       Correction for the dimension of the manifold.
%
%   Apr.  2, 2015 (NB):
%       Replaced trace(A'*B) by A(:)'*B(:) (equivalent but faster).
%
%   Apr. 17, 2018 (NB):
%       Removed dependence on lyap.
%
%   Sep.  6, 2018 (NB):
%       Removed M.exp() as it was not implemented.
%
%   June 7, 2019  (EM):
%       Added M.dist, M.exp, M.log and M.invretr.

    M.name = @() sprintf('%dx%d s.p.s.d. matrices of rank %d', n, k);

    M.dim = @() k*n - k*(k-1)/2;

    % Euclidean metric on the embedding space
    M.inner = @(x, eta, zeta) inner(x, eta, zeta);
    function output = inner(x, eta, zeta)
        tr1 = sum(sum(eta.S.*zeta.S));
        tr2 = sum(sum(eta.Z.*zeta.Z)); 
        output = double(tr1 + 2*tr2);
    end

    M.norm = @(x, eta) sqrt(M.inner(x, eta, eta));

%     M.dist = @(Y, Z) norm(logarithm(Y,Z),'fro');

    M.typicaldist = @() 10*k;

    M.proj = @projection;
    function output = projection(x, eta)
        % Projection onto the tangent space
        VZ = x.V' * eta.Z;
        output.S = sym(eta.S + VZ + VZ');
        output.Z = eta.Z - x.V * VZ;
    end

    M.tangent = M.proj;    

%     M.exp = @exponential;
%     function Ynew = exponential(Y, eta, t)
%         if nargin < 3
%             t = 1.0;
%         end
%         Ynew = Y + t*eta;
%     end

    M.retr = @retr;
        function x = retr(x, eta, t)
        if nargin == 3
            eta = veccomb(x, t, eta, 0, eta);
        end
        [Vp, R] = qr(eta.Z, 0);
        p = size(R,1);
        function [V1, D1] = pm(M)
           [V1, D1] = svd(sym(M));
           V1 = V1(:,1:p);
           D1 = D1(1:p,1:p);
        end
        [V1, D1] = pm([x.D + eta.S, R'; R, zeros(p, p)]);
        x.V = [x.V, Vp] * V1;
        x.D = D1;
    end
    
%     M.log = @logarithm;
%     function eta = logarithm(Y, Z)
%         YtZ = Y'*Z;
%         [U, ~, V] = svd(YtZ);
%         Qt = V*U';
%         eta = Z*Qt - Y;
%     end

%     M.invretr = M.log;

%     M.egrad2rgrad = @(Y, eta) eta;
%     M.ehess2rhess = @(Y, egrad, ehess, U) M.proj(Y, ehess);

    % Notice that the hash of two equivalent points will be different...
%     M.hash = @(Y) ['z' hashmd5(Y(:))];

    M.rand = @() random(n, k);
    function x = random(n, k)
        [x.V, ~] = qr(randn(n,k),0);
        x.D = diag(sort(rand(k,1), 'descend'));
    end

    M.randvec = @randomvec;
    function eta = randomvec(x)
        S = randn(size(x.D));   
        eta.S = sym(S);
        Z = randn(size(x.V));
        Z = Z - x.V*(x.V'*Z);
        eta.Z = Z;    
        neta = M.norm(x, eta);
        eta = veccomb(x, 1/neta, eta, 0, eta);
    end    

    M.lincomb = @veccomb;
    function output = veccomb(x, a, eta, b, zeta)
    % output = a*eta + b*zeta;
        output.S = a*eta.S + b*zeta.S;
        output.Z = a*eta.Z + b*zeta.Z;    
    end

    M.zerovec =@zerovec;
    function zero_eta = zerovec(x)
        zero_eta.S = zeros(size(x.D));
        zero_eta.Z = zeros(size(x.V));
    end

%     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 = @() false;

end

function x = sym(x)
    x  = 0.5 * (x + x');
end
