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

0001 function low_rank_tensor_completion_embedded()
0002 % Given partial observation of a low rank tensor (possibly including noise),
0003 % attempts to complete it.
0004 %
0005 % function low_rank_tensor_completion_embedded()
0006 %
0007 % NOTE: Tensor Toolbox version 2.6 or higher is required for this factory:
0008 % see https://www.tensortoolbox.org/ or https://gitlab.com/tensors/tensor_toolbox
0009 %
0010 % This example demonstrates how to use the geometry factory for the
0011 % embedded submanifold of fixed-rank tensors in Tucker format:
0012 % fixedranktensorembeddedfactory.
0013 %
0014 % This geometry is described in the article
0015 % "A Riemannian trust-region method for low-rank tensor completion"
0016 % Gennadij Heidel and Volker Schulz, doi:10.1002/nla.2175.
0017 %
0018 % This can be a starting point for many optimization problems of the form:
0019 %
0020 % minimize f(X) such that rank(X) = [r1 ... rd], size(X) = [n1 ... nd].
0021 %
0022 % Important: to keep this example short, the code for the cost function,
0023 % gradient and Hessian do not properly exploit sparsity of the
0024 % observations, which leads to significant slow-downs for large tensors.
0025 % This example file should be considered a starting point for more
0026 % sophisticated implementations.
0027 %
0028 % Input:  None. This example file generates random data with noise.
0029 %
0030 % Output: None.
0031 %
0032 % Please cite the Manopt and Matlab Tensor Toolbox papers as well as the
0033 % research paper:
0034 %     @Article{heidel2018riemannian,
0035 %       Title   = {A {R}iemannian trust-region method for low-rank tensor completion},
0036 %       Author  = {G. Heidel and V. Schulz},
0037 %       Journal = {Numerical Linear Algebra with Applications},
0038 %       Year    = {2018},
0039 %       Volume  = {23},
0040 %       Number  = {6},
0041 %       Pages   = {e1275},
0042 %       Doi     = {10.1002/nla.2175}
0043 %     }
0044 %
0045 % See also: fixedranktensorembeddedfactory
0046 
0047 % This file is part of Manopt: www.manopt.org.
0048 % Original author: Gennadij Heidel, January 24, 2019.
0049 % Contributors:
0050 % Change log:
0051 
0052     if ~exist('tenrand', 'file')
0053         fprintf('Tensor Toolbox version 2.6 or higher is required.\n');
0054         return;
0055     end
0056 
0057     % Random data generation with pseudo-random numbers from a
0058     % uniform distribution on [0, 1].
0059     tensor_dims = [60 40 20];
0060     core_dims = [8 6 5];
0061     total_entries = prod(tensor_dims);
0062     d = length(tensor_dims);
0063     
0064     % Standard deviation of normally distributed noise.
0065     % Set sigma to 0 to test the noise-free case.
0066     sigma = 0.1;
0067     
0068     % Generate a random tensor A of size n1-by-...-by-nd of rank (r1, ..., rd).
0069     U = cell(1, d);
0070     R = cell(1, d);
0071     for i = 1:d
0072         [U{i}, R{i}] = qr(randn(tensor_dims(i), core_dims(i)), 0);
0073     end
0074 
0075     Z.U = R;
0076     Z.G = tenrand(core_dims);
0077     Core = ttm(Z.G, Z.U);
0078 
0079     Y.U = U;
0080     Y.G = Core;
0081     A = ttm(Core, Y.U);
0082     
0083     % Add noise to low-rank tensor
0084     A = A + sigma*tensor(randn(tensor_dims));
0085     
0086     
0087     % Generate a random mask P for observed entries:
0088     % P(i, j, k) = 1 if the entry (i, j, k) of A is observed,
0089     %            = 0 otherwise.
0090     fraction = 0.1; % Fraction of observed entries.
0091     nr = round(fraction * total_entries);
0092     ind = randperm(total_entries);
0093     ind = ind(1 : nr);
0094     P = false(tensor_dims);
0095     P(ind) = true;
0096     % Hence, we observe the nonzero entries in PA:
0097     P = tensor(P);
0098     PA = P.*A; 
0099     % Note that an efficient implementation would require evaluating A as a
0100     % sparse tensor only at the indices of P.
0101 
0102     
0103     
0104     % Pick the submanifold of tensors of size n1-by-...-by-nd of
0105     % multilinear rank (r1, ..., rd).
0106     problem.M = fixedranktensorembeddedfactory(tensor_dims, core_dims);
0107     
0108     
0109     % Define the problem cost function.
0110     % The store structure is used to reduce full tensor evaluations.
0111     % Again: proper handling of sparse tensors would dramatically reduce
0112     % the computation time for large tensors. This file only serves as a
0113     % simple starting point. See help for the Tensor Toolbox regarding
0114     % sparse tensors. Same comment for gradient and Hessian below.
0115     problem.cost = @cost;
0116     function [f, store] = cost(X, store)
0117         if ~isfield(store, 'PXmPA')
0118             Xfull = full(X.X);
0119             store.PXmPA = P.*Xfull - PA;
0120         end
0121         f = .5*norm(store.PXmPA)^2;
0122     end
0123 
0124     % Define the Euclidean gradient of the cost function, that is, the
0125     % gradient of f(X) seen as a function of X without rank restrictions.
0126     problem.egrad =  @egrad;
0127     function [g, store] = egrad(X, store)
0128         if ~isfield(store, 'PXmPA')
0129             Xfull = full(X.X);
0130             store.PXmPA = P.*Xfull - PA;
0131         end
0132         g = store.PXmPA;
0133     end
0134     
0135     % Define the Euclidean Hessian of the cost at X along a vector eta.
0136     problem.ehess = @ehess;
0137     function H = ehess(X, eta)
0138         ambient_H = problem.M.tangent2ambient(X, eta);
0139         H = P.*ambient_H;
0140     end
0141     
0142     % Options
0143     X0 = problem.M.rand();
0144     options.maxiter = 3000;
0145     options.maxinner = 100;
0146     options.maxtime = inf;
0147     options.storedepth = 3;
0148     % Target gradient norm
0149     options.tolgradnorm = 1e-8*problem.M.norm(X0, getGradient(problem, X0));
0150 
0151     % Minimize the cost function using Riemannian trust-regions
0152     Xtr = trustregions(problem, X0, options);
0153 
0154     % Display some quality metrics for the computed solution
0155     Xtrfull = full(Xtr.X);
0156     fprintf('||X-A||_F / ||A||_F = %g\n', norm(Xtrfull - A)/norm(A));
0157     fprintf('||PX-PA||_F / ||PA||_F = %g\n', norm(P.*Xtrfull - PA)/norm(PA));
0158     
0159 end