Block term decomposition¶

Generating pseudorandom factor matrices and core tensor¶

In Tensorlab, a BTD with factor matrices \(\mat{U}^{(r,n)}\) and core tensors \(\ten{S}^{(r)}\) is stored as a cell of cells U. The \(r\)th entry U{r} contains the \(r\)th term in the BTD. Each term is a cell containing the \(N\) factor matrices and the core tensor for the \(r\)th term. This means U{r}{n} is the factor matrix \(\mat{U}^{(r,n)}\) of size \(I_n\times J_r^{(n)}\) for \(n = 1,\ldots,N\), and U{r}{N+1} is the core tensor \(\ten{S}^{(r)}\) of size \(J_r^{(1)} \times J_r^{(2)} \times \cdots \times J_r^{(N)}\). Hence, U = {U{1},U{2},...,U{R}} is a cell with \(R\) cells, each having \(N+1\) entries.

To generate pseudorandom factor matrices and core tensors corresponding to a BTD, btd_rnd can be used. The size of the tensor is given by size_tens, and the size of the \(R\) core tensors is given by the cell size_core, i.e., size_core{r} is the size of the \(r\)th core tensor. For example:

size_tens = [17 19 21];
size_core = {[3 5 7],[6 3 5],[4 3 4]};
U = btd_rnd(size_tens,size_core);

By default btd_rnd generates U{r}{n} using randn, after which the factor matrices U{r}{1:N} are orthogonalized. Other generators can also be provided with options given as a structure or as key-value pairs, e.g.:

% using a structure
options.Real = @rand;
options.Imag = @rand;
[U,S] = btd_rnd(size_tens,size_core,options);

% using key-value pairs
[U,S] = btd_rnd(size_tens, size_core, 'Real', @rand, 'Imag', @rand);

Each term is normalized such that the factors \(\mat{U}^{(r,n)}\) have orthonormal columns and such that \(\ten{S}^{(r)}\) is all-orthogonal and has ordered multilinear singular values. The normalization can be disabled by setting the Unitary and Normalize options to false.

Generating the associated full tensor¶

Given a cell array of block terms U = {U{1},U{2},...}, its associated full tensor T can be computed with

T = btdgen(U);

Using basic tensor operations, the tensor T can also be computed as

R = length(U);
N = length(U{1})-1;
size_tens = cellfun('size',U{1}(1:N),1);
T = zeros(size_tens);
for r = 1:R
     T = T + tmprod(U{r}{N+1},U{r}(1:N),1:N);
end

The ful method can also be used to generate the full tensor. Both the ful and the btdgen method allow the generation of subtensors:

U = btd_rnd([10 11 12], {[3 4 2], [2,2,1]});
T = btdgen(U); % or, equivalently
T = ful(U);


t1 = btdgen(U, [1 1 3 2 1 5 8]); % using linear indices
frob(t1 - T([1 1 3 2 1 5 8]))    % is zero up to machine precision

t2 = ful(U, ':', 9, 3);          % using subscripts (':' is a string!)
frob(t2 - T(:,9,3))              % is zero up to machine precision

With a specific algorithm¶

In contrast to the CPD, the LMLRA and the decomposition in multilinear rank-\((L_r,L_r,1)\) terms, no specialized high-level algorithm is available for a general BTD. This means that the user has to generate his or her own initialization and select a specific (family of) algorithm(s) such as btd_nls or btd_minf to compute the BTD given the initialization.

To compute a BTD of a dense, sparse, incomplete or structured tensor T given an initialization U0, the command btd_nls(T,U0) can be used. For example,

% Generate pseudorandom U and associated full tensor T.
size_tens = [17 19 21];
size_core = {[2 2 2],[2 2 2],[2 2 2]};
U = btd_rnd(size_tens,size_core);
T = btdgen(U);

% Generate an initialization Uinit and compute the BTD with nonlinear least squares.
U0   = noisy(U,20);
Uhat = btd_nls(T,U0);

generates a tensor T as the sum of three real rank-\((2,2,2)\) tensors and then computes its BTD.

The decomposition in multilinear rank-\((L_r,L_r,1)\) terms is a special case of a BTD for third-order tensors, in which each term \(r\) has multilinear rank \((L_r,L_r,1)\). For this kind of BTD, there exist efficient algorithms which are preferred over the generic BTD algorithms. See Decomposition in multilinear rank- terms for more information.

Setting the options¶

The selected algorithm may be customizable by supplying the method with additional options (see the relevant help for more information). These options can be given using an options structure, or using key-value pairs, e.g.,

options.Display = 5;     % Show convergence progress every 5 iterations.
options.MaxIter = 200;   % Set stop criteria.
options.TolFun  = 1e-12;
options.TolX    = 1e-12;
Uhat = btd_nls(T,U0,options);

Viewing the algorithm output¶

The selected method also returns output specific to the algorithm, such as the number of iterations and the algorithm’s objective function value. To obtain this information, capture the second output:

[Uhat,output] = btd_nls(T,U0,options);

and inspect its fields, for example by plotting the objective function value:

semilogy(0:output.iterations,sqrt(2*output.fval));
xlabel('iteration');
ylabel('frob(btdres(T,U))');
grid on;

Computing the error¶

The residual between a tensor \(\ten{T}\) and its BTD defined by the factors \(\mat{\hat{U}}^{(r,n)}\) and the core tensors \(\ten{\hat{S}}^{(r)}\) (stored in Uhat) can be computed with

res = btdres(T,Uhat);

If the tensor is dense, then the result is equivalent to btdgen(Uhat)-T. The relative error between the tensor and its BTD can be computed as

relerr = frob(btdres(T,Uhat))/frob(T); % or
relerr = frobbtdres(T,Uhat)/frob(T); % also works for structured tensors