Multilinear singular value decomposition and low multilinear rank approximation¶

Tags: mlsvd, lmlra

The multilinear singular value decomposition is, as the term indicates, a multilinear generalization of the matrix singular value decomposition. In this demo, we provide some insight in this decomposition and the low multilinear rank approximation. A Matlab implementation of this demo is given in the demo_mlsvd.m file, which can be found in the demo folder of Tensorlab or downloaded here.

For matrices, the singular value decomposition \(\mat{M} = \mat{U}\mat{\Sigma}\transpose{\mat{V}}\) is well known. Using tensor-matrix products, this decomposition can be written as

\[\begin{aligned} \mat{M} = \mat{\Sigma}\tmprod_1\mat{U}\tmprod_2\mat{V}.\end{aligned}\]

The matrix \(\mat{\Sigma}\) is a diagonal matrix, and the matrices \(\mat{U}\) and \(\mat{V}\) are orthogonal matrices. A generalization of this SVD is the multilinear singular value decomposition (MLSVD). In the literature, one can also find the names higher-order SVD (HOSVD) and Tucker decomposition. The term Tucker decomposition has evolved over the years and is now often used in a more general sense. The MLSVD of a third-order tensor can be written as

\[\begin{aligned} \ten{T} = \ten{S}\tmprod_1\mat{U}^{(1)} \tmprod_2\mat{U}^{(2)} \tmprod_3\mat{U}^{(3)}. \label{eq:mlsvd}\end{aligned}\]

Like in the matrix case where \(\mat{U}\) and \(\mat{V}\) are orthonormal bases for the column and row space, respectively, the MLSVD computes \(N\) orthonormal bases \(\mat{U}^{(n)}\), \(n=1,\ldots,N\) for the \(N\) different subspaces of mode-\(n\) vectors. As is illustrated in Figure 74, these orthonormal bases have exactly the same interpretation in as the matrix case. Essentially, the MLSVD considers a given tensor as \(N\) sets of mode-\(n\) vectors, and computes the matrix SVD of these sets. The multilinear rank of the tensor \(\ten{T}\) corresponds to the \(N-\text{tuplet}\) consisting of the dimensions of the different subspaces.

Figure 74: Illustration of the multilinear singular value decomposition of a multilinear rank-\((R_1,R_2,R_3)\) tensor and the different spaces.

The noiseless case¶

First we study the properties of the MLSVD for the exact case. Random factor matrices and a random core tensor are used to construct a multilinear rank-\((2, 2, 3)\) tensor \(\ten{T}\). Next the MLSVD is computed using mlsvd.

S = randn(2,2,3);
U = {rand(4,2), rand(5,2), rand(6,3)};
T = lmlragen(U,S);
[Ue,Se,sve] = mlsvd(T)

We can check the properties of the MLSVD:

The first two, two and three columns of respectively Ue{1}, Ue{2}, and Ue{3} are a basis for the column, row and fiber subspaces of \(\ten{T}\), respectively. (subspace computes the angle between subspaces: a zero angle indicates that the spaces are identical.)
```
subspace(Ue{1}(:,1:2),U{1})
subspace(Ue{2}(:,1:2),U{2})
subspace(Ue{3}(:,1:3),U{3})
```
mlsvd-kladfq¶
```
 4.3356e-16
 4.9257e-16
 6.2976e-16
```

The remaining two, three and three columns of Ue{1}, Ue{2}, and Ue{3} are orthogonal bases for the complements of the column, row and fiber subspaces of \(\ten{T}\). These subspaces are orthogonal to the column, row, and fiber subspaces, respectively:

U{1}'*Ue{1}(:, 3:end)
U{2}'*Ue{2}(:, 3:end)
U{3}'*Ue{3}(:, 4:end)

nlsvdadlfj¶

 ans =
    1.0e-15 *
     0.2220    0.0555
     0.0694    0.4580
 ans =
    1.0e-15 *
    -0.0139    0.0230    0.0173
     0.1110    0.0607    0.1110
 ans =
    1.0e-15 *
    -0.1110    0.2637         0
    -0.0971    0.2463    0.0902
     0.0555   -0.2220   -0.1527

Similar to the matrix SVD, the multilinear singular values sve indicate the importance of a basis vector in a particular mode. In contrast to the SVD, the core tensor \(\ten{S}\) is not diagonal. The multilinear singular values are defined as the Frobenius norms of the slices of order \(N-1\) for the different modes, i.e. the norm of the subtensor created by fixing one index:
```
nrm = zeros(size(Se,1),1);
for i = 1:size(Se,1)
    nrm(i) = frob(Se(i,:,:));
end
[sve{1} nrm] % mode-1 singular values
```
lkasdls¶
```
 ans =
     2.0722    2.0722
     1.1272    1.1272
     0.0000    0.0000
     0.0000    0.0000
```
The multilinear rank of a tensor can be determined by looking at the multilinear singular values. As shown in Figure 75, the multilinear singular values \(\sigma_{r}^{(n)}\) become zero for \(r > R_n\), in which \(R_n\) is the mode-\(n\) rank.
```
for n = 1:3
    subplot(1,3,n);
    semilogy(sve{n},'.-');
    xlim([1 length(sve{n})])
end
```

Figure 75: Multilinear singular values in the noiseless case.

Finally, the all-orthogonality of the core tensor is investigated. This property means that all slices of order \(N-1\) are mutually orthogonal to each other. We test this for the second mode by computing the inner product between the slices. To do this efficiently, we unfold the tensor in the second mode.
```
M = tens2mat(Se, 2);
M*M'
```
lkasdsad¶
```
  4.7448         0    0.0000   -0.0000   -0.0000
       0    0.8197    0.0000    0.0000    0.0000
  0.0000    0.0000    0.0000    0.0000    0.0000
 -0.0000    0.0000    0.0000    0.0000    0.0000
 -0.0000    0.0000    0.0000    0.0000    0.0000
```
We indeed see that all the off-diagonal entries are zero. Note that the diagonal contains the squared multilinear singular values for the second mode.
```
[sve{2} sve{2}.^2]
```
lkasfadsfadsfasfdls¶
```
 2.1782    4.7448
 0.9054    0.8197
 0.0000    0.0000
 0.0000    0.0000
 0.0000    0.0000
```

Noisy case¶

The equality in Equation holds only for the exact case. When noise is added, the multilinear rank-\((2, 2, 3)\) decomposition only approximates the tensor. Here, the influence of the noise is investigated. In this example, noise is added such that the signal-to-noise ratio is 30dB.

Tn = noisy(T, 30);

Again, a MLSVD can be computed from this noisy tensor. The multilinear singular values are plotted in Figure 76. Estimating the multilinear rank is now a more difficult task.

[Uen,Sen,sven] = mlsvd(Tn);
for n = 1:3
    subplot(1,3,n);
    semilogy(sven{n},'.-');
    xlim([1 length(sve{n})])
end

Figure 76: Multilinear singular values in the noisy case.

From the construction, we know that the multilinear rank is \((2, 2, 3)\). We can only keep the respective vectors and part of the core tensor:

Uentrunc{1} = Uen{1}(:,1:2);
Uentrunc{2} = Uen{2}(:,1:2);
Uentrunc{3} = Uen{3}(:,1:3);
Sentrunc = Sen(1:2, 1:2, 1:3);

The approximation error will not be zero now. The relative error is:

format long
frob(lmlragen(Uentrunc,Sentrunc)-Tn)/frob(Tn)

lmlra-1¶

 0.049896088936784

In contrast to the matrix case, this approximation may not be optimal in terms of relative Frobenius norm error, however. The lmlra method attempts to find the tensor of multilinear rank \((2, 2, 3)\) that best approximates \(\ten{T}\) in a least squares sense.

[Uopt, Sopt] = lmlra(Tn, [2 2 3]);

We again compute the error:

frob(lmlragen(Uopt,Sopt)-Tn)/frob(Tn)

lmlra-2¶

 0.049884982735855

The error has indeed improved. (For other tensors the difference may be more substantial.) The computational cost has increased, however. This cost can be somewhat reduced by taking the MLSVD solution as initialization:

[Uopt, Sopt] = lmlra(Tn, Uentrunc, Sentrunc);

The relative approximation error can also be improved by keeping more multilinear singular values, and hence by taking larger subspaces.

Uentrunc2{1} = Uen{1}(:,1:3);
Uentrunc2{2} = Uen{2}(:,1:3);
Uentrunc2{3} = Uen{3}(:,1:3);
Sentrunc2 = Sen(1:3, 1:3, 1:3);

frob(lmlragen(Uentrunc2,Sentrunc2)-Tn)/frob(Tn)

lmlra-3¶

 0.042189469973458

While the error indeed decreases and the MLSVD is computationally cheaper than lmlra, taking larger subspaces may not be desirable in particular applications, e.g., in multidimensional harmonic retrieval.