positivefactory.m

function M = positivefactory(m, n)
% Manifold of m-by-n matrices with positive entries; scale invariant metric
%
% function M = positivefactory(m)
% function M = positivefactory(m, n)
%
% A point X on the manifold M is represented as a matrix X of size mxn with
% all individual entries real, strictly positive. By default, n = 1.
%
% A tangent vector at X is represented as a matrix of the same size as X.
% Entries of tangent vectors are free (in particular, not necessarily
% positive.)
%
% The Riemannian metric for each individual entry is the bi-invariant
% metric for positive scalars, as a particular case of the bi-invariant
% metric for positive definite matrices studied in Chapter 6 of the book
%
%    "Positive definite matrices" by Rajendra Bhatia,
%    Princeton University Press, 2007.
%
% The Riemannian structure of M is obtained as the Cartesian product of the
% geometry for mxn positive real numbers.
%
% It should be stressed that matrices with one or more zero entries do not
% belong to this manifold: they appear to be infinitely far away as a
% result of the metric scaling like X.^(-1). Thus, if the solutions of an
% optimization problem have entries equal to zero, these solutions are not
% attainable on the manifold, which is likely to create serious numerical
% issues. This geometry is best used when the solutions of the optimization
% problem are indeed entry-wise positive, yet may have very different
% scales (with some entries being very small, and some entries being very
% large, relatively.)
%
% See also: sympositivedefinitefactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec 03, 2017.

    if ~exist('n', 'var') || isempty(n)
        n = 1;
    end
    
    M.name = @() sprintf('Element-wise positive %dx%d matrices', m, n);
    
    M.dim = @() m*n;
        
    % The metric is the scale invariant metric for scalars.
    M.inner = @myinner;
    function innerproduct = myinner(X, eta, zeta)
        innerproduct = (eta(:)./X(:))'*(zeta(:)./X(:));
    end
   
    M.norm = @(X, eta) sqrt(myinner(X, eta, eta));
    
    M.dist = @(X, Y) norm(log(Y./X), 'fro');
    
    M.typicaldist = @() sqrt(m*n);
    
    M.egrad2rgrad = @egrad2rgrad;
    function eta = egrad2rgrad(X, eta)
        eta = X.*(eta).*X;
    end
    
    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)
        % Directional derivatives of the Riemannian gradient
        Hess = X.*(ehess).*X + 2*(eta.*(egrad).*X);
        
        % Correction factor for the non-constant metric
        Hess = Hess - (eta.*(egrad).*X);
    end
    
    % Since this manifold is an open subset of R^(nxm), the tangent space
    % at any X on M is all of R^(nxm).
    M.proj = @(X, eta) eta;
    
    M.tangent = M.proj;
    M.tangent2ambient = @(X, eta) eta;
    
    M.retr = @exponential;
    
    M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        % It is unclear whether this is the numerically most stable way to
        % implement this operation. If you run into trouble with this
        % factory, please get in touch on the forum.
        Y = (X.*(exp((t*eta)./X)));
    end
    
    M.log = @logarithm;
    function H = logarithm(X, Y)
        % Same comment about numerical stability as for exp.
        H = (X.*(log(Y./X)));
    end
    
    M.hash = @(X) ['z' hashmd5(X(:))];
    
    % Generate a random element-wise positive matrix following a
    % certain distribution. The particular choice of a distribution is of
    % course arbitrary, and specific applications might require different
    % ones.
    M.rand = @random;
    function X = random()
        X = exp(randn(m, n));
    end
    
    % Generate a uniformly random unit-norm tangent vector at X.
    M.randvec = @randomvec;
    function eta = randomvec(X)
        eta = randn(m, n).*X;
        nrm = M.norm(X, eta);
        eta = eta / nrm;
    end
    
    M.lincomb = @matrixlincomb;
    
    M.zerovec = @(X) zeros(m, n);
    
    
    M.transp = @(X1, X2, eta) eta;
    
    % For reference, a proper transporter is given here, following
    % work by Sra and Hosseini: "Conic geometric optimisation on the
    % manifold of positive definite matrices".
    % This is not used by default. To force the use of this transport,
    % execute "M.transp = M.paralleltransp;" on your M returned by the
    % present factory.
    M.paralleltransp = @parallel_transport;
    M.isotransp = M.paralleltransp;
    function zeta = parallel_transport(X, Y, eta)
        zeta = eta.*Y./X;
    end
    
    % vec and mat are not isometries, because of the unusual inner metric.
    M.vec = @(X, U) U(:);
    M.mat = @(X, u) reshape(u, m, n);
    M.vecmatareisometries = @() true;
    
end