Home > examples > low_rank_tensor_completion_embedded.m

# low_rank_tensor_completion_embedded

## PURPOSE

Given partial observation of a low rank tensor (possibly including noise),

## SYNOPSIS

function low_rank_tensor_completion_embedded()

## DESCRIPTION

``` 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}
}

## CROSS-REFERENCE INFORMATION

This function calls:
• fixedranktensorembeddedfactory Manifold of tensors with fixed multilinear rank in Tucker format
• full FULL Convert TTeMPS tensor to full array
• norm NORM Norm of a TT/MPS tensor.
• round ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• full FULL Convert TTeMPS tensor to full array
• norm NORM Norm of a TT/MPS block-mu tensor.
• round ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• full FULL Convert TTeMPS_op operator to full array
• round ROUND Approximate TTeMPS operator within a prescribed tolerance.
• trustregions Riemannian trust-regions solver for optimization on manifolds.
This function is called by:

## SOURCE CODE

```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 %
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.
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;