Basic use of Tensorlab¶

Construction of a tensor with random entries

A random tensor can be constructed by using the basic randn and rand commands in Matlab. For example, we can construct a tensor \(\ten{T}\in\R^{10\times 20\times 30}\) with entries randomly sampled from a Gaussian distribution as follows:

Construction of a tensor as a sum of R rank-1 terms

To construct a third-order tensor \(\ten{T}\in\R^{I\times J \times K}\) which can be written as a sum of \(R\) rank-1 terms, we start with the construction of the three factor matrices \(\mat{U}^{(1)}\in\R^{I\times R}\), \(\mat{U}^{(2)}\in\R^{J\times R}\) and \(\mat{U}^{(3)}\in\R^{K\times R}\), corresponding to the three different modes (see Canonical polyadic decomposition). This can be done with the command cpd_rnd, returning a cell which collects these matrices:

Construction of a tensor obtained by a multilinear transformation of a core tensor

We start with the construction of the core tensor and the transformation matrices using the lmlra_rnd command. The tensor can be obtained through the tensor-matrix product of the core tensor with the different matrices, for which one can use the lmlragen command.

Perturbation of a tensor by noise

A tensor can be perturbed by Gaussian noise with a user-defined signal-to-noise-ratio using the command noisy:

Visualization

Tensorlab offers four commands for tensor visualization: slice3, surf3, voxel3 and visualize . The plots of the first two commands contain sliders to go through the different slices. The voxel3 command plots the elements of a tensor as voxels whose color and opacity are proportial to the value of the elements. The plot has two sliders for the threshold parameter and the degree parameter. Figures 55, 56 and 57 visualize the tensor constructed in the following block of code:

t = linspace(0,1,60).';
T = cpdgen({sin(2*pi*t+pi/4),sin(4*pi*t+pi/4),sin(6*pi*t+pi/4)});
figure(); slice3(T);
figure(); surf3(T);
figure(); voxel3(T);

Figure 55: Visualization of a tensor using slice3.

Figure 56: Visualization of a tensor using surf3.

Figure 57: Visualization of a tensor using voxel3.

Finally, the visualize method can be used to ‘walk through’ the data and to compare tensor a decomposed tensor with its decomposition. For example:

t = linspace(0,1,60).';
U = {sin(2*pi*t+pi/4),sin(4*pi*t+pi/4),sin(6*pi*t+pi/4)};
T = noisy(cpdgen(U), 30);
visualize(U, 'original', T);

The result can be seen in Figure 58. The check boxes can be used to select the dimensions to plot. The sliders or the text fields allow a user to select which slices are plotted. More examples involving the visualize method can be found in the user guide.

Figure 58: Visualization of a tensor using visualize.

Rank estimation

If the rank of a tensor is not known or cannot be obtained from prior knowledge, the command rankest from Tensorlab may be useful. It estimates the corner of the L-curve obtained by sweeping over a number of rank-1 terms. This corner is taken as an estimate of the rank of the tensor. If there is no output argument, an L-curve as in Figure 59 is plotted.

figure();
U = cpd_rnd([10 20 30],4);
T = cpdgen(U);
rankest(T);

Figure 59: The plot generated by the rankest command.

The basic CPD command

The basic cpd algorithm requires only a tensor and a value for the rank:

R = 4;
Uest = cpd(T,R)

An extra output argument can be provided with [U,output] = cpd(T,R) to obtain information about the convergence of the algorithm, among other things.

Error measures

The relative error on the estimates of the factor matrices can be calculated with the command cpderr. It takes into account the permutation and scaling indeterminacies. The command returns the relative errors (for each factor matrix), the permutation matrix, the scaling matrices and the permuted and scaled factor matrices:

[relerrU,P,D,Uest] = cpderr(U,Uest)
% D is a cell collecting the different scaling matrices

The relative error on the estimate of the tensor can be determined using frob and cpdres, with the latter calculating the residuals:

relerrT = frob(cpdres(T,Uest))/frob(T)

Compression

In the compression step, the tensor is compressed to a smaller tensor using the mlsvd_rsi command by default, which applies a randomized SVD algorithm in combination with randomized subspace iterations. Compression is only performed if it is expected to lead to a significant reduction in computational complexity. The user can alter the compression method with the Compression option by providing a compression technique:

Uest = cpd(T,R,'Compression',@lmlra_aca);

In this example, a low multilinear rank approximation with adaptive cross-approximation is used. The user can also disable the compression step ( note that there will not be a refinement step then):

Uest = cpd(T,R,'Compression',false);

Initialization

To initialize the computation of the CPD of the compressed tensor, cpd uses a method based on the generalized eigenvalue decomposition (implemented in cpd_gevd) if the rank does not exceed the second largest dimension of the tensor. Otherwise, a random initialization is used. The initialization method can be altered with the Initialization option:

Uest = cpd(T,R,'Initialization',@cpd_rnd);

In this example, a random initialization is chosen. The user can also provide particular factor matrices Uinit for the initialization:

Uest = cpd(T,Uinit);

When specific factor matrices are provided, the compression step is skipped.

Core algorithms

Tensorlab provides several core algorithms for the computation of a CPD. First, there are optimization-based methods such as cpd_als (alternating least squares), cpd_minf (unconstrained nonlinear optimization) and cpd_nls (nonlinear least squares). Second, there are (semi-)algebraic methods such as cpd3_sd and cpd3_sgsd. The former reduces the computation to the simultaneous diagonalization of symmetric matrices, also in cases where only one factor matrix has full column rank, while the latter reduces the computation to the estimation of orthogonal factors in a simultaneous generalized Schur decomposition.

By default, cpd_nls is used for the CPD of the compressed tensor and for the refinement. The user can alter the choice of algorithm with the Algorithm and Refinement options. By default, Refinement is equal to Algorithm. The following example uses an alternating least squares approach for the CPD of the compressed tensor and for the refinement:

Uest = cpd(T,R,'Algorithm',@cpd_als);

Core algorithm parameters

Parameters can be passed to the core algorithm with the AlgorithmOptions option for cpd. Tensorlab provides default choices for all parameters.

A first parameter is Display, which determines after how many iterations intermediate results are printed. The first column of the printed results shows the objective function values. The second and third column give the relative change in the objective function values and the step sizes, respectively.

There are three parameters that affect the stopping criterion. The parameter Maxiter determines the maximum number of iterations. The parameters TolX and TolFun are tolerances for the step size and for the relative change in objective function value, respectively. For the NLS algorithm, the CGMaxIter option determines the maximum number of conjugate gradient iterations in each NLS iteration. The following piece of code illustrates how these parameters can be changed for the cpd_nls and cpd_als algorithms:

R = 4;
T = cpdgen(cpd_rnd([10 20 30],R));
Uinit = cpd_rnd([10 20 30],R);

options = struct;
options.Compression = false;
options.Algorithm = @cpd_nls;
options.AlgorithmOptions.Display = 1;
options.AlgorithmOptions.MaxIter = 100;      % Default 200
options.AlgorithmOptions.TolFun = eps^2;     % Default 1e-12
options.AlgorithmOptions.TolX = eps;         % Default 1e-6
options.AlgorithmOptions.CGMaxIter = 20;     % Default 15
[Uest_nls,output_nls] = cpd(T,Uinit,options);

options = struct;
options.Compression = false;
options.Algorithm = @cpd_als;
options.AlgorithmOptions.Display = 1;
options.AlgorithmOptions.MaxIter = 100;      % Default 200
options.AlgorithmOptions.TolFun = eps^2;     % Default 1e-12
options.AlgorithmOptions.TolX = eps;         % Default 1e-6
[Uest_als,output_als] = cpd(T,Uinit,options);

A convergence plot can be obtained from the link in the command terminal output, or by plotting the objective function values from output.Algorithm (and output.Refinement):

figure();
semilogy(output_nls.Algorithm.fval); hold all;
semilogy(output_als.Algorithm.fval);
ylabel('Objective function'); xlabel('Iteration');
title('Convergence plot'); legend('cpd\_nls','cpd\_als')

Figure 60 shows the output of the code block above and illustrates first-order convergence for the ALS algorithm and second-order convergence for the NLS algorithm.

Figure 60: Convergence plot for the cpd_nls and cpd_als algorithm.