symfixedrankYYfactory.m

function M = symfixedrankYYfactory(n, k)
% Manifold of n-by-n symmetric positive semidefinite matrices of rank k.
%
% function M = symfixedrankYYfactory(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'. 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.
%
% If you wish to optimize with bounded rank rather than fixed rank, see
% the tool manoptlift with burermonteirolift.
% 
% 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}
%     }
%
%     @article{MasAbs2020,
%         author = {Massart, Estelle and Absil, P.-A.},
%         title = {Quotient Geometry with Simple Geodesics for the Manifold of Fixed-Rank Positive-Semidefinite Matrices},
%         journal = {SIAM Journal on Matrix Analysis and Applications},
%         volume = {41},
%         number = {1},
%         pages = {171--198},
%         year = {2020},
%         doi = {10.1137/18M1231389}
%     }
% 

% See also: elliptopefactory spectrahedronfactory symfixedrankYYcomplexfactory
%           manoptlift burermonteirolift

% 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('YY'' quotient manifold of %dx%d psd matrices of rank %d', n, k);

    M.dim = @() k*n - 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.dist = @(Y, Z) norm(logarithm(Y,Z),'fro');

    M.typicaldist = @() 10*k;

    M.proj = @projection;
    function etaproj = projection(Y, eta)
        % 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.exp = @exponential;
    function Ynew = exponential(Y, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Ynew = Y + t*eta;
    end

    M.retr = M.exp;
    
    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;
    function Y = random()
        Y = randn(n, k);
    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