function low_rank_tensor_completion_embedded()
% Given partial observation of a low rank tensor (possibly including noise),
% attempts to complete it.
%
% function low_rank_tensor_completion_embedded()
%
% NOTE: Tensor Toolbox version 2.6 or higher is required for this factory:
% see https://www.tensortoolbox.org/ or https://gitlab.com/tensors/tensor_toolbox
%
% This example demonstrates how to use the geometry factory for the
% embedded submanifold of fixed-rank tensors in Tucker format:
% fixedranktensorembeddedfactory.
%
% This geometry is described in the article
% "A Riemannian trust-region method for low-rank tensor completion"
% Gennadij Heidel and Volker Schulz, doi:10.1002/nla.2175.
%
% This can be a starting point for many optimization problems of the form:
%
% minimize f(X) such that rank(X) = [r1 ... rd], size(X) = [n1 ... nd].
%
% Important: to keep this example short, the code for the cost function,
% gradient and Hessian do not properly exploit sparsity of the
% observations, which leads to significant slow-downs for large tensors.
% This example file should be considered a starting point for more
% sophisticated implementations.
%
% Input:  None. This example file generates random data with noise.
% 
% Output: None.
%
% Please cite the Manopt and Matlab Tensor Toolbox papers as well as the
% research paper:
%     @Article{heidel2018riemannian,
%       Title   = {A {R}iemannian trust-region method for low-rank tensor completion},
%       Author  = {G. Heidel and V. Schulz},
%       Journal = {Numerical Linear Algebra with Applications},
%       Year    = {2018},
%       Volume  = {23},
%       Number  = {6},
%       Pages   = {e1275},
%       Doi     = {10.1002/nla.2175}
%     }
%
% See also: fixedranktensorembeddedfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Gennadij Heidel, January 24, 2019.
% Contributors: 
% Change log:

    if ~exist('tenrand', 'file')
        fprintf('Tensor Toolbox version 2.6 or higher is required.\n');
        return;
    end

    % Random data generation with pseudo-random numbers from a 
    % uniform distribution on [0, 1].  
    tensor_dims = [60 40 20];
    core_dims = [8 6 5];
    total_entries = prod(tensor_dims);
    d = length(tensor_dims);
    
    % Standard deviation of normally distributed noise.
    % Set sigma to 0 to test the noise-free case.
    sigma = 0.1;
    
    % Generate a random tensor A of size n1-by-...-by-nd of rank (r1, ..., rd).
    U = cell(1, d);
    R = cell(1, d);
    for i = 1:d
        [U{i}, R{i}] = qr(randn(tensor_dims(i), core_dims(i)), 0);
    end

    Z.U = R;
    Z.G = tenrand(core_dims);
    Core = ttm(Z.G, Z.U);

    Y.U = U;
    Y.G = Core;
    A = ttm(Core, Y.U);
    
    % Add noise to low-rank tensor
    A = A + sigma*tensor(randn(tensor_dims));
    
    
    % Generate a random mask P for observed entries:
    % P(i, j, k) = 1 if the entry (i, j, k) of A is observed,
    %            = 0 otherwise.
    fraction = 0.1; % Fraction of observed entries.
    nr = round(fraction * total_entries);
    ind = randperm(total_entries);
    ind = ind(1 : nr);
    P = false(tensor_dims);
    P(ind) = true;
    % Hence, we observe the nonzero entries in PA:
    P = tensor(P);
    PA = P.*A; 
    % Note that an efficient implementation would require evaluating A as a
    % sparse tensor only at the indices of P.

    
    
    % Pick the submanifold of tensors of size n1-by-...-by-nd of
    % multilinear rank (r1, ..., rd).
    problem.M = fixedranktensorembeddedfactory(tensor_dims, core_dims);
    
    
    % Define the problem cost function.
    % The store structure is used to reduce full tensor evaluations.
    % Again: proper handling of sparse tensors would dramatically reduce
    % the computation time for large tensors. This file only serves as a
    % simple starting point. See help for the Tensor Toolbox regarding
    % sparse tensors. Same comment for gradient and Hessian below.
    problem.cost = @cost;
    function [f, store] = cost(X, store)
        if ~isfield(store, 'PXmPA')
            Xfull = full(X.X);
            store.PXmPA = P.*Xfull - PA;
        end
        f = .5*norm(store.PXmPA)^2;
    end

    % Define the Euclidean gradient of the cost function, that is, the
    % gradient of f(X) seen as a function of X without rank restrictions.
    problem.egrad =  @egrad;
    function [g, store] = egrad(X, store)
        if ~isfield(store, 'PXmPA')
            Xfull = full(X.X);
            store.PXmPA = P.*Xfull - PA;
        end
        g = store.PXmPA;
    end
    
    % Define the Euclidean Hessian of the cost at X along a vector eta.
    problem.ehess = @ehess;
    function H = ehess(X, eta)
        ambient_H = problem.M.tangent2ambient(X, eta);
        H = P.*ambient_H;
    end
    
    % Options
    X0 = problem.M.rand();
    options.maxiter = 3000;
    options.maxinner = 100;
    options.maxtime = inf;
    options.storedepth = 3;
    % Target gradient norm
    options.tolgradnorm = 1e-8*problem.M.norm(X0, getGradient(problem, X0));

    % Minimize the cost function using Riemannian trust-regions
    Xtr = trustregions(problem, X0, options);

    % Display some quality metrics for the computed solution
    Xtrfull = full(Xtr.X);
    fprintf('||X-A||_F / ||A||_F = %g\n', norm(Xtrfull - A)/norm(A));
    fprintf('||PX-PA||_F / ||PA||_F = %g\n', norm(P.*Xtrfull - PA)/norm(PA));
    
end
