function M = fixedrankfactory_2factors_preconditioned(m, n, k)
% Manifold of m-by-n matrices of rank k with two factor quotient geometry.
%
% function M = fixedrankfactory_2factors_preconditioned(m, n, k)
%
% This geometry is tuned to least-squares problems such as low-rank matrix
% completion with ell-2 loss.
%
% A point X on the manifold is represented as a structure with two
% fields: L and R. The matrices L (m-by-k) and R (n-by-k) are 
% full column-rank matrices such that X = L*R'.
%
% Tangent vectors are represented as a structure with two fields: L, R.
% 
% Please cite the Manopt paper as well as the research paper:
%     @Techreport{mishra2012optimized,
%       Title   = {A {R}iemannian geometry for low-rank matrix completion},
%       Author  = {Mishra, B. and Adithya Apuroop, K. and Sepulchre, R.},
%       Journal = {Arxiv preprint arXiv:1211.1550},
%       Year    = {2012}
%     }
%
%
% See also: fixedrankembeddedfactory fixedrankfactory_2factors fixedrankfactory_3factors_preconditioned

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:
%
%    Apr. 4, 2015 (BM):
%        Cosmetic changes including avoiding storing the inverse of a
%        k-by-k matrix.
%
%    Sep. 6, 2018 (NB):
%        Removed M.exp() as it was not implemented.
  
    
    M.name = @() sprintf('LR''(tuned to least square problems) quotient manifold of %dx%d matrices of rank %d', m, n, k);
    
    M.dim = @() (m+n-k)*k;
    
    
    % Some precomputations at the point X to be used in the inner product 
    % (and pretty much everywhere else).
    function X = prepare(X)
        if ~all(isfield(X,{'LtL','RtR'}))
            L = X.L;
            R = X.R;
            X.LtL = L'*L;
            X.RtR = R'*R;
        end
    end
    
    
    % The choice of metric is motivated by symmetry and 
    % tuned to least-squares cost function.
    M.inner = @iproduct;
    function ip = iproduct(X, eta, zeta)
        X = prepare(X);
        ip = trace(X.RtR*(eta.L'*zeta.L)) + trace(X.LtL*(eta.R'*zeta.R)); % Scaled metric
    end
    
    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));
    
    M.dist = @(x, y) error('fixedrankfactory_2factors_preconditioned.dist not implemented yet.');
    
    M.typicaldist = @() 10*k;
    
    M.egrad2rgrad = @egrad2rgrad;
    function rgrad = egrad2rgrad(X, egrad)
        X = prepare(X);
        
        % Riemannian gradient
        rgrad.L = egrad.L/X.RtR;
        rgrad.R = egrad.R/X.LtL;
    end
    
    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)
        X = prepare(X);
        
        % Riemannian gradient.
        rgrad = egrad2rgrad(X, egrad);
        
        % Directional derivative of the Riemannian gradient.
        Hess.L = ehess.L/X.RtR - 2*egrad.L*(X.RtR \ (symm(eta.R'*X.R) / X.RtR));
        Hess.R = ehess.R/X.LtL - 2*egrad.R*(X.LtL \ (symm(eta.L'*X.L) / X.LtL));
        
        % We still need a correction factor for the non-constant metric.
        Hess.L = Hess.L + rgrad.L*(symm(eta.R'*X.R)/X.RtR) + eta.L*(symm(rgrad.R'*X.R)/X.RtR) - X.L*(symm(eta.R'*rgrad.R)/X.RtR);
        Hess.R = Hess.R + rgrad.R*(symm(eta.L'*X.L)/X.LtL) + eta.R*(symm(rgrad.L'*X.L)/X.LtL) - X.R*(symm(eta.L'*rgrad.L)/X.LtL);
        
        % Project on the horizontal space.
        Hess = M.proj(X, Hess);
    end
    
    M.proj = @projection;
    function etaproj = projection(X, eta)
        X = prepare(X);
        
        % Projection onto the horizontal space.
        Lambda = 0.5*((eta.R'*X.R)/X.RtR  -   X.LtL\(X.L'*eta.L));
        etaproj.L = eta.L + X.L*Lambda;
        etaproj.R = eta.R - X.R*Lambda';
    end
    
    M.tangent = M.proj;
    
    M.tangent2ambient = @(X, eta) eta;
    
    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Y.L = X.L + t*eta.L;
        Y.R = X.R + t*eta.R;
        
        % Numerical conditioning step: a simpler version.
        % We need to ensure that L and R are do not have very relative
        % skewed norms.
        
        scaling = norm(X.L, 'fro')/norm(X.R, 'fro');
        scaling = sqrt(scaling);
        Y.L = Y.L / scaling;
        Y.R = Y.R * scaling;
        
        % These are reused in the computations of gradient and Hessian.
        Y = prepare(Y);
    end
    
    
    M.hash = @(X) ['z' hashmd5([X.L(:) ; X.R(:)])];
    
    M.rand = @random;
    
    function X = random()
        X.L = randn(m, k);
        X.R = randn(n, k);
    end
    
    M.randvec = @randomvec;
    function eta = randomvec(X)
        eta.L = randn(m, k);
        eta.R = randn(n, k);
        eta = projection(X, eta);
        nrm = M.norm(X, eta);
        eta.L = eta.L / nrm;
        eta.R = eta.R / nrm;
    end
    
    M.lincomb = @lincomb;
    
    M.zerovec = @(X) struct('L', zeros(m, k),'R', zeros(n, k));
    
    M.transp = @(x1, x2, d) projection(x2, d);
    
    % vec and mat are not isometries, because of the scaled inner metric.
    M.vec = @(X, U) [U.L(:) ; U.R(:)];
    
    M.mat = @(X, u) struct('L', reshape(u(1:(m*k)), m, k), ...
        'R', reshape(u((m*k+1):end), n, k));
    
    M.vecmatareisometries = @() false;
    
    % Auxiliary functions
    symm = @(M) .5*(M+M');
end

% Linear combination of tangent vectors.
function d = lincomb(x, a1, d1, a2, d2) %#ok<INUSL>
    
    if nargin == 3
        d.L = a1*d1.L;
        d.R = a1*d1.R;
    elseif nargin == 5
        d.L = a1*d1.L + a2*d2.L;
        d.R = a1*d1.R + a2*d2.R;
    else
        error('Bad use of fixedrankfactory_2factors_preconditioned.lincomb.');
    end
    
end





