Tensorization techniques¶

Definition¶

Second-order Hankelization transforms a vector \(\vec{v}\in\C^{N}\) into a Hankel matrix \(\mat{H}\in\C^{I \times J}\) with \(N = I+J-1\), defined by:

Hence, each anti-diagonal (or skew-diagonal) of the matrix \(\mat{H}\) is constant, as illustrated in the following example.

In general, \(K\)th-order Hankelization transforms a vector \(\vec{v}\in\C^{N}\) into a higher-order Hankel tensor \(\mat{H}\in\C^{I_1 \times \cdots \times I_K}\) with \(N = I_1+\cdots+I_K-K+1\), defined by:

A Hankel tensor has constant anti-hyperplanes.

From the previous example, it is clear that the original vector \(\vec{v}\) appears in the Hankel tensor as follows. The vector \(\vec{v}\) is divided into \(K\) segments, of which consecutive segments overlap with one entry. The first segment appears as the first fiber along the first mode of the Hankel tensor, and the other segments appear as the last fibers along the other modes.

Hankelization in Tensorlab¶

Hankelization can be carried out with hankelize. We discuss some of its features. The Order, Dim, Full, FullLimit and PermToFirst options are used by the other tensorization commands as well. The Sizes and Ind options are only used by hankelize.

1. First, the user can specify the order of the Hankelization with the Order option. For a Hankelization of order \(K\), a vector is tensorized to a Hankel tensor of order \(K\). By default, a second-order Hankelization is carried out.

Todo

If Order = 2, then a vector \(\vec{v}\in\C^{N}\) is mapped to a matrix \(\mat{V}\in\C^{I\times J}\):

V = hankelize(v);  % By default, 'Order' is 2

If Order = 4, then a vector \(\vec{v}\in\C^{N}\) is mapped to a fourth-order tensor \(\ten{V}\in\C^{I_1\times I_2\times I_3 \times I_4}\):

V = hankelize(v,'Order',4);

2. The Sizes option can be used to set the values \(I_1,\ldots,I_{K-1}\). The size \(I_K\) is chosen as \(I_K = N-(I_1+\cdots+I_{K-1})+K-1\). By default, \(I_1,\ldots,I_K\) are chosen as close to each other as possible.

Instead of Sizes, the Ind option can be used as well. The vector \(\vec{v}\) is partitioned into \(K\) segments \(\vec{v}_{1:d_1},\vec{v}_{d_1:d_2},\ldots,\vec{v}_{d_{K-1}:N}\) indicating the first fiber along the first mode of the Hankel tensor and the last fibers along the other modes. The Ind option determines the indices \(d_1,\ldots,d_{K-1}\). The Ind and Sizes options are related by Ind = cumsum(Sizes)-(0:Order-1).

Todo

The third-order Hankel tensor \(\ten{H}\in\C^{4\times 3\times 2}\) from (1) can be obtained from a vector \(\vec{v}\in\C^{7}\) with

H = hankelize(v,'Sizes',[4 3]);

Equivalently, it can be obtained with

H = hankelize(v,'Ind',[4 6]);

as \(\vec{v}_{1:4}\) is used for the first fiber along the first mode, and \(\vec{v}_{4:6}\) and \(\vec{v}_{6:7}\) are used for the last fibers along the second and third modes, respectively.

3. When matrices or tensors are Hankelized, the user can specify along which mode the tensorization is carried out with the Dim option. All fibers along this mode are Hankelized. Consider an \(M\)th-order tensor \(\ten{T}\in\C^{N_1 \times \cdots \times N_M}\). If \(\ten{T}\) is tensorized along mode \(n\) with tensorization order \(K\), hankelize returns a tensor \(\ten{H}\) of size \(N_1 \times \cdots \times N_{n-1} \times I_{1} \times \cdots \times I_{K} \times N_{n+1} \times \cdots \times N_{M}\). By default, a tensorization is carried out along the first mode. Let us illustrate with some examples.

Todo

Let

\[\begin{split}\mat{M} = \begin{bmatrix} 1 & 6 \\ 2 & 7 \\ 3 & 8 \\ 4 & 9 \\ 5 & 10 \end{bmatrix}.\end{split}\]

If Dim = 1 and Order = 2, each column of \(\mat{M}\in\C^{5\times 2}\) is Hankelized. The different Hankel matrices are stacked along the third mode of a tensor \(\ten{Z}\in\C^{3\times 3 \times 2}\):

Z = hankelize(M);  % By default, 'Dim' is 1 and 'Order' is 2

\(\ten{Z}\) has slices

\[\begin{split}\mat{Z}_{::1} = \left[\begin{array}{ccc} 1 & 2 & 3 \\ 2 & 3 & 4 \\ 3 & 4 & 5\end{array}\right],\ \mat{Z}_{::2} = \left[\begin{array}{ccc} 6 & 7 & 8 \\ 7 & 8 & 9 \\ 8 & 9 & 10\end{array}\right].\end{split}\]

Todo

Consider a matrix \(\mat{M}\in\C^{21\times 3}\). If the tensorization order is 4, hankelize returns a fifth-order tensor \(\ten{Z}\in\C^{6\times 6 \times 6 \times 6 \times 3}\) containing Hankel tensors of size \(6\times 6 \times 6\times 6\):

M = randn(21,3);
Z = hankelize(M,'Order',4);
size(Z)

fweqrdfwsddfsafd¶

ans =
  [6   6   6   6   3]

Todo

Consider a tensor \(\ten{T}\in\C^{10 \times 20 \times 10}\) which is Hankelized along the second mode (Dim = 2) using third-order Hankelization (Order = 3). The hankelize command returns a fifth-order tensor \(\ten{Z}\in\C^{10\times 7 \times 8 \times 7 \times 10}\):

T = randn(10,20,10);
Z = hankelize(T,'Dim',2,'Order',3);
size(Z)

oijsodijdiojdi¶

ans =
  [10   7   8   7   10]

In this example, each mode-2 vector of length 20 is tensorized to a Hankel tensor of size \(7 \times 8 \times 7\).

4. The Hankelized data can become very large due to the curse of dimensionality. Instead of working with the full Hankelized tensor, an efficient representation can be used for calculations by exploiting the Hankel structure. This efficient representation can be obtained from the second output argument of hankelize. Alternatively, it can be retrieved from the first output argument if the Full option is set to false. Note that the latter avoids the construction and storage of the full Hankelized tensor alltogether. The CPD and decomposition in multilinear rank-\((L_r,L_r,1)\) terms are computed much faster when using efficient representations of Hankelized datasets than when using the full tensors.

Todo

The third-order Hankel tensor of a vector \(\vec{v}\) of size \(1000\times 1\) has size \(334\times 334 \times 334\) by default. Using double precision, this is a dataset of approximately 284 MB. Using

[~,H] = hankelize(v,'Order',3)
% or
H = hankelize(v,'Order',3,'full',false)

returns the following efficient representation of \(\ten{H}\), which only requires 9.8 MB:

 H =
       type: 'hankel'
        val: [1000x1 double]
        dim: 1
      order: 3
        ind: [334 667]
 ispermuted: 0
repermorder: [1 2 3]
       size: [334 334 334]
    subsize: [1x1 struct]

Note that hankelize automatically detects if the full tensor would become too large. When the storage limit FullLimit is exceeded, the full tensor is not constructed, but the efficient representation is returned. The FullLimit option is set to 1 (expressed in GB) by default, but can be adjusted by the user.

The full tensor can be obtained from the efficient representation with the ful command.

5. The PermToFirst option can be used to indicate whether the output modes of the tensorization should be the first modes of the resulting tensor. Consider an \(M\)th-order tensor \(\ten{T}\in\C^{N_1 \times \cdots \times N_M}\) which is tensorized along the \(n\)th mode with tensorization order \(K\). By default, PermToFirst is false, and hankelize returns a tensor of size \(N_1 \times \cdots \times N_{n-1} \times I_{1} \times \cdots \times I_{K} \times N_{n+1} \times \cdots \times N_{M}\). If PermToFirst is set to true, hankelize returns a tensor of size \(I_{1} \times \cdots \times I_{K} \times N_{1} \times \cdots \times N_{n-1} \times N_{n+1} \times \cdots \times N_{M}\).

Todo

Consider a tensor \(\ten{T}\in\C^{10 \times 20 \times 10}\), which is Hankelized along the second mode (Dim = 2) using third-order tensorization (Order = 3). By default, a tensor \(\ten{Z}\in\C^{10\times 7 \times 8 \times 7 \times 10}\) is obtained. However, when setting PermToFirst = true, hankelize returns a tensor \(\ten{Z}\in\C^{7 \times 8 \times 7 \times 10 \times 10}\):

Z = hankelize(T,'Dim',2,'Order',3,'PermToFirst',1);

This option can be useful for decompositions such as the decomposition in multilinear rank-\((L_r,L_r,1)\) terms, where \(L_r\) may correspond to the rank of a Hankel representation.

Dehankelization in Tensorlab¶

The original data can be extracted from a Hankel matrix or Hankel tensor with the dehankelize command, which provides a number of options:

1. The Order option indicates the order of the dehankelization. By default, Order = 2, and each second-order slice in the first two modes is dehankelized. This is consistent with the default use of second-order tensorization in hankelize. If a vector should be extracted from a \(K\)th-order Hankel tensor, Order should be set to \(K\).

Todo

Consider a Hankel tensor \(\ten{H}\in\C^{10 \times 10 \times 10}\). Dehankelization of \(\ten{H}\) with

M = dehankelize(H);

returns a matrix \(\mat{M}\) of size \(19\times 10\), while setting Order = 3 returns a vector of length 28:

v = dehankelize(H,'Order',3);

2. The Dims option can be used to indicate along which modes the dehankelization is carried out.

Todo

Consider a tensor \(\ten{H}\in\C^{5 \times 10 \times 5 \times 10 \times 5}\). The command

T = dehankelize(H,'Dims',[1 3 5]);

returns a tensor \(\ten{T}\) of size \(13 \times 10 \times 10\).

3. Instead of using the Dims option, a combination of the Dim options and Order can be used for consecutive modes along which the detensorization should take place. The result with Dim = d and Order = K is the same as with Dims = d:d+K-1.
4. The user can indicate how the data is extracted using the Method option. We discuss some possibilities:

• The first column and the last row of a Hankel matrix are retained if Method is set to fibers. Considering a \(K\)th-order Hankel tensor, the first fiber in the first mode, and the last fibers along the other modes are retained.

• The data can also be estimated by averaging along the anti-diagonals if Method is set to the string mean or the function handle @mean. This can be useful if the data is noisy and/or if the data does not exactly have Hankel structure. For tensors, the average along the anti-hyperplanes is computed.

  Todo

  Let

  \[\begin{split}\mat{H} = \begin{bmatrix} 1 & 2.1 & 3.1 \\ 1.9 & 3 & 4.1 \\ 2.9 & 4 & 5.1 \\ 3.9 & 4.9 & 6 \end{bmatrix}.\end{split}\]

  The command

  v = dehankelize(H,'Method','fibers');
  

  returns the concatenated first column and last row, \(\vec{v} = \transpose{[1,1.9,2.9,3.9,4.9,6]}\), while

  v = dehankelize(H,'Method','mean');   % or @mean
  

  returns the averaged anti-diagonals \(\vec{v} = \transpose{[1,2,3,4,5,6]}\).

• The user can provide an alternative for the averaging along the anti-diagonals or anti-hyperplanes by passing a function handle to Method, such as @median:

  v = dehankelize(H,'Method',@median);
  

By default, Method is set to @mean as the result is more robust than with the fibers method. Note that @median is computationally more demanding.

The dehankelize command works very efficiently on a polyadic representation of a Hankel tensor. Internally, it avoids the construction of the possibly large full tensor. Hence, while the results are the same,

v = dehankelize({U1,U2,U3});

is more time- and memory-efficient if the factor matrices are large than

T = cpdgen({U1,U2,U3});
v = dehankelize(T);

It is possible to dehankelize several polyadic representations at the same time by concatenating the factor matrices (provided that the dimensions agree) and indicating how the rank-1 terms are grouped with the L option.

W = cpd_rnd([10 10],5);       % Concatenated U and V
v = dehankelize(W);
size(v)

ans =
   [19   1]

W = cpd_rnd([10 10],5);       % Concatenated U and V
M = dehankelize(W,'L',[3 2]);
size(M)

ans =
   [19   2]

Application example¶

We consider blind signal separation for a mixture of exponentials to illustrate the use of Hankelization and dehankelization. Four exponential functions, sampled in 101 points \(\vec{t}=[0,0.01,0.02,\ldots,1]\), are mixed such that three observed signals are obtained. The goal is to determine the mixing coefficients and the exponentials using only the observed signals.

This can be achieved by Hankelizing the observed signals along the mode of the observations. The resulting tensor with Hankel slices is decomposed using a CPD. Estimates of the source signals can be obtained by dehankelizing the estimated source Hankel matrices. Note that, since there are more source signals than observed signals, the source signals cannot be estimated by inverting the estimated mixing matrix. More information on this application can for instance be found in [28].

% Signal construction
t = linspace(0,1,101).';
S = exp(t*[1 2 3 4]);
M = randn(4,3);
X = S*M;
% Separation
H = hankelize(X,'Dim',1);
U = cpd(H,4);
Sest = dehankelize(U(1:2),'L',[1 1 1 1]);
% Error calculation
cpderr(S,Sest)

ans =
   5.3055e-12

Definition¶

Consider a vector \(\vec{f}\in\C^{N}\) and a vector of evaluation points \(\vec{t}\in\C^{N}\). Let us partition \([1,2,\ldots,N]\) in two vectors \(\vec{x}\in\C^{I}\) and \(\vec{y}\in\C^{J}\) with \(N = I+J\). For example, an interleaved partitioning (\(\vec{x} = [1,3,\ldots]\) and \(\vec{y} = [2,4,\ldots]\)) or block partitioning (\(\vec{x} = [1,2,\ldots,I]\) and \(\vec{y} = [I+1,\ldots,N]\)) can be used. The corresponding Löwner matrix \(\mat{L}\in\C^{I\times J}\) is defined as follows:

Löwner matrices can be useful for dealing with rational functions, as there is a close relationship between the degree of the rational function and the rank of the Löwner matrix [41]. We do not discuss higher-order Löwner tensors in this section.

Löwnerization in Tensorlab¶

Löwnerization can be carried out with the loewnerize command. It can be applied using similar options as in hankelize:

1. The Order option indicates the order of the Löwnerization and is set to two by default. For example, L = loewnerize(f,'Order',3).
2. The Dim option indicates the mode along which the Löwnerization needs to be carried out. By default, Dim = 1. For example, L = loewnerize(f,'Dim',2).
3. Setting the Full option to false avoids the construction of the full Löwner matrix or tensor. If the partitioned evaluation points in \(\vec{x}\) and \(\vec{y}\) are equidistant, time- and memory-efficient calculations can be performed using the efficient representation. The full Löwnerized tensor can be obtained from the efficient representation with the ful command.
4. The PermToFirst option can be used to indicate whether the output modes of the tensorization should be the first modes of the result. For example, L = loewnerize(f,'Dim',2,'PermToFirst',true) can be used to make the Löwner matrices appear in the first and second mode. By default, PermToFirst = false.

Two additional options are provided. These are illustrated in the examples below.

5. The T option can be used to set the evaluation points. By default, T = 1:N.
6. The Ind option can be used to set the partitioning. It is assumed to be a cell array with \(K-1\) distinct sets for \(K\)th-order Löwnerization. The \(K\)th set consists of the remaining indices. By default, an interleaved partitioning is used. If the cell contains only a single array, the cell brackets can be omitted.

f = [6 7 8 9];
% By default, 'Order' is 2, 'T' is [1 2 3 4] and 'Ind' is [1,3].
L = loewnerize(f)

L =

   1     1
   1     1

f = [1 4 9 16];
L = loewnerize(f,'T',[2 4 6 8],'Ind',[1 2])

L =

   2     2.5
   2.5   3

Delöwnerization in Tensorlab¶

Delöwnerization extracts the original data from a Löwner matrix or tensor. This can be done with deloewnerize which solves a linear system [41]. The command works in a similar way as dehankelize. The user can set the Dims option to indicate the modes along which the detensorization should be carried out (Dims = [1 2] by default). Alternatively, a combination of Dim and Order can be used for consecutive modes. In the case of delöwnerization, the T option can be used to specify the evaluation points, and the Ind option can be used to specify the partitioning (interleaved by default).

The deloewnerize command also works on polyadic representations using factor matrices. The L option can be used to indicate separate delöwnerizations, similar to dehankelize.

The result of the delöwnerization of a Löwner matrix is determined up to a constant offset. The deloewnerize command fixes this indeterminacy by returning a vector with zero mean.

Application example¶

We consider blind signal separation for a mixture of rational functions to illustrate the use of Löwnerization and delöwnerization. Four rational functions, sampled in 101 points \(\vec{t}=[0,0.01,0.02,\ldots,1]\), are mixed such that three observed signals are obtained. The goal is to determine the mixing coefficients and the rational functions using only the observed signals.

This can be achieved by Löwnerizing the observed signals along the mode of the observations. The resulting tensor with Löwner slices is decomposed using a CPD. Estimates of the source signals can be obtained by delöwnerizing the estimated source Löwner matrices. Note that, since there are more source signals than observed signals, the source signals cannot be estimated by inverting the estimated mixing matrix. More information on this application can for instance be found in [41].

% Signal construction
t = linspace(0,1,100).';
% Each column of S consists of a rational function of the form 1/(t-p)
S = 1./bsxfun(@minus,t,[0.2 0.4 0.6 0.8]);
M = randn(4,3);
X = S*M;
% Separation
H = loewnerize(X,'Dim',1,'T',t);
U = cpd(H,4);
Sest = deloewnerize(U(1:2),'L',[1 1 1 1]);
% Error calculation: compare Sest to S with zero-mean columns
cpderr(bsxfun(@minus,S,mean(S,1)),Sest)

ans =
   3.1050e-15

Definition¶

Given a vector \(\vec{v}\in\C^{N}\), second-order segmentation without overlapping segments consists in reshaping \(\vec{v}\) into a matrix \(\mat{S}\in\C^{I\times J}\) with \(N=IJ\). Each column of \(\mat{S}\) is then called a segment of \(\vec{v}\) of length \(I\). For example, the segmentation of a vector \(\vec{v} = [1,2,\ldots,6]\) with \(I = 3\) and \(J = 2\) results in a matrix \(\mat{S}\) of size \(3\times 2\):

When the order of the segmentation is increased, each segment is segmented again, resulting in a tensor \(\ten{S}\in\C^{I_1\times I_2 \times\cdots\times I_K}\) for a \(K\)th-order segmentation of a vector \(\vec{v}\in\C^{N}\) with \(N=I_1I_2\cdots I_K\).

Tensorlab also provides the possibility of using overlapping segments, as discussed below.

Segmentation in Tensorlab¶

Segmentation can be obtained with the segmentize command. Similar to hankelize and loewnerize, the Order option can be used to define the order of the segmentation and the Dim option can be used to indicate the mode along which the tensorization is carried out. Setting the Full option to false avoids the construction of the full tensor, while the PermToFirst option allows the user to indicate whether the output modes of the tensorization should be the first modes modes of the resulting tensor.

A number of additional options are available for segmentize:

1. The segment length can be set to \(I\) for second-order segmentation with the Segsize option. The number of segments \(J\) is automatically chosen as \(J = \lfloor N/I\rfloor\).

In general, for \(K\)th-order segmentation, the Segsize option can be used to set the values \(I_1,I_2,\ldots,I_{K-1}\). The value \(I_K\) is automatically chosen as \(I_K = \lfloor N/(I_1I_2\cdots I_{K-1})\rfloor\).

Todo

Consider a matrix \(\mat{M}\in\C^{60\times 6}\) which is tensorized by third-order segmentation along the first mode. If \(I_1 = 3\) and \(I_2 = 4\), then the resulting tensor \(\ten{S}\) has size \(3\times 4 \times 5 \times 6\):

M = randn(60,6);
S = segmentize(M,'Segsize',[3 4]);
size(S)

qwerewdsfsdaf¶

ans =
  [3   4   5   6]

2. Equivalently, for second-order segmentation, the Nsegments option can be used to set the number of segments \(J\). The segment length \(I\) is then chosen as \(I = \lfloor N/J\rfloor\).

In general, for \(K\)th-order segmentation, Nsegments can be used to set the values \(I_2,I_3,\ldots,I_{K}\). The value \(I_1\) is automatically chosen as \(I_1 = \lfloor N/(I_2I_3\cdots I_K)\rfloor\).

Todo

The tensor \(\ten{S}\) from the previous example can also be constructed as follows:

M = randn(60,6);
S = segmentize(M,'Nsegments',[4 5]);
size(S)

wefqoijdsijijo¶

ans =
  [3   4   5   6]

3. The UseAllSamples option can be set to true if an error needs to be thrown if not all entries of the data are used (i.e., if \(N\neq I_1I_2\cdots I_K\) for \(K\)th-order segmentation). By default, UseAllSamples is set to false and no error is thrown.
4. The user can specify the shift between consecutive segments with the advanced Shift option. If the shift is smaller than the segment length for second-order segmentation, consecutive segments overlap. If the shift is larger than the segment length, some data entries are discarded. By default, Shift is equal to the segment length for second-order segmentation so that there is no overlap. The number of segments is chosen such that as many data values are used as possible.

Todo

Consider the segmentization of the vector \(\vec{v}=[1,2,\ldots,9]\) with a shift of 2:

v = 1:9;
S = segmentize(v,'Shift',2)

15sdfa48df81dfs4851d¶

S =
   1     3     5     7
   2     4     6     8
   3     5     7     9

For general \(K\)th-order segmentation, Shift is a vector with \(K-1\) entries. Recall the process of \(K\)th-order segmentation with recursive segmentation and without overlap: segments of length \(I_1I_2\cdots I_{K-r}\) are constructed in the \(r\)th segmentation step which are then further segmented in the next step. The shift in step \(r\) is indicated by the \(u\)th entry of the Shift vector with \(u=K-r\). Note that the segment length in segmentation step \(r\) depends not only on \(I_1,I_2,\ldots,I_{K-r}\) but on the first \(K-r\) shifts as well. By default, Shift is equal to the cumulative product of the values \(I_1,I_2,\ldots,I_{K-1}\) so that there is no overlap.

Todo

Consider the vector \(\vec{v}=[1,2,\ldots,18]\). By default, using \(I_1 = I_2 = 3\) and \(I_3 = 2\), the vector \(\vec{v}\) is first segmented into a matrix of size \(9\times 2\). Each column is then segmented to a \(3\times 3\) matrix. This results in a third-order tensor \(\ten{S}^{(1)}\in\C^{3\times 3 \times 2}\):

v = 1:18;
S1 = segmentize(v,'Segsize',[3 3]);

The tensor \(\ten{S}^{(1)}\) has slices

\[\begin{split}\mat{S}_{::1}^{(1)} = \left[\begin{array}{ccc} 1 & 4 & 7\\ 2 & 5 & 8\\ 3 & 6 & 9 \end{array}\right],\ \mat{S}_{::2}^{(1)} = \left[\begin{array}{ccc} 10 & 13 & 16 \\ 11 & 14 & 17 \\ 12 & 15 & 18 \end{array}\right].\end{split}\]

In this case, the Shift option is set by default to Shift = [3 9]. When setting Shift to [3 8], the segments of length 9 with starting indices 1 and 9 (\(\vec{v}_{1:9}\) and \(\vec{v}_{9:17}\)) are created first. This results in a matrix of size \(9\times 2\). Each column is segmented again into segments of length 3 without overlap as the first element of Shift is 3. segmentize returns a tensor \(\ten{S}^{(2)}\in\C^{3 \times 3 \times 2 }\):

v = 1:18;
S2 = segmentize(v,'Segsize',[3 3],'Shift',[3 8])

qwe54eweqr45eq41eq¶

 S2(:,:,1) =

   1     4     7
   2     5     8
   3     6     9

 S2(:,:,2) =

   9    12    15
  10    13    16
  11    14    17

Note that there is an overlapping element with value 9 along the third mode. Alternatively, with Segsize = [3 3] and Shift = [2 7], one obtains a tensor \(\ten{S}^{(3)}\) with slices

\[\begin{split}\mat{S}_{::1}^{(3)} = \left[\begin{array}{ccc} 1 & 3 & 5 \\ 2 & 4 & 6 \\ 3 & 5 & 7 \end{array}\right],\ \mat{S}_{::2}^{(3)} = \left[\begin{array}{ccc} 8 & 10 & 12 \\ 9 & 11 & 13 \\ 10 & 12 & 14 \end{array}\right].\end{split}\]

Note that there are overlapping elements along the second mode, but no overlap along the third mode. In this process, \(\vec{v}\) is first segmented into consecutive segments of length 7 to obtain a matrix of size \(7\times 2\). Each column is then segmented again in segments of length 3 with an overlap of 1 as the first shift element is 2.

If Segsize = [3 3] and Shift = [2 5], segments of length 7 are constructed first with an overlap of 2. Three such segments can be constructed from \(\vec{v}\) to obtain a matrix of size \(7\times 3\). Each column is then segmented again to segments of length 3 with an overlap of 1, to obtain a tensor \(\ten{S}^{(4)}\) of size \(3\times 3\times 3\):

v = 1:18;
S4 = segmentize(v,'Segsize',[3 3],'Shift',[2 5]);

\(\ten{S}^{(4)}\) has slices

\[\begin{split}\mat{S}_{::1}^{(4)} = \left[\begin{array}{ccc} 1 & 3 & 5 \\ 2 & 4 & 6 \\ 3 & 5 & 7 \end{array}\right],\ \mat{S}_{::2}^{(4)} = \left[\begin{array}{ccc} 6 & 8 & 10 \\ 7 & 9 & 11 \\ 8 & 10 & 12 \end{array}\right],\ \mat{S}_{::3}^{(4)} = \left[\begin{array}{ccc} 11 & 13 & 15 \\ 12 & 14 & 16 \\ 13 & 15 & 17 \end{array}\right].\end{split}\]

When setting Shift = [1,...,1] with \(K-1\) entries, the segments constructed in each segmentation step overlap almost entirely. The result is the same as that of \(K\)th-order Hankelization.

v = 1:6;
S = segmentize(v,'Segsize',[3 3],'Shift',[1 1])

S(:,:,1) =

   1     2     3
   2     3     4
   3     4     5

S(:,:,2) =

   2     3     4
   3     4     5
   4     5     6

v = 1:6;
H = hankelize(v,'sizes',[3 3]);

Hence, Hankelization with hankelize is a special case of segmentation with segmentize.

Note that if the Full option in segmentize is set to false, an efficient representation of the segmented tensor is returned. This can be useful to avoid the redundancy due to overlapping segments. Tensorlab does not yet exploit this structure in tensor computations and decompositions, however. The full tensor can be obtained from the efficient representation with the ful command.

Desegmentation in Tensorlab¶

Desegmentation extracts the original data from the segmentized results. The desegmentize command works in a similar way as dehankelize. The user is able to indicate the modes along which the detensorization is carried out (with Dims, or with a combination of Dim and Order), and to set the Shift option. The Method option indicates the method to be used for the extraction of the data, such as single fibers, means, medians, ...

v = 1:7;
S = segmentize(v,'Segsize',[3 3],'Shift',[1 2])

S(:,:,1) =

    1     2     3
    2     3     4
    3     4     5

S(:,:,2) =

    3     4     5
    4     5     6
    5     6     7

vest = desegmentize(S,'Shift',[1 2])

vest =
    1
    2
    3
    4
    5
    6
    7

Definition¶

Decimation is a tensorization technique similar to segmentation. Second-order decimation consists of the reshaping of a vector \(\vec{v}\in\C^{N}\) into a matrix \(\mat{D}\in\C^{I\times J}\) such that each column of \(\mat{D}\) is a subsampled version of \(\vec{v}\) with a subsampling factor \(J\).

Higher-order decimation consists in further reshaping the obtained vectors, resulting in a tensor \(\ten{D}\in\C^{I_1\times I_2 \times \cdots \times I_K}\) for the \(K\)th-order decimation of a vector \(\vec{v}\in\C^{N}\) with \(N = I_1I_2\cdots I_K\).

Decimation in Tensorlab¶

Decimation can be obtained with the decimate command. The command works in a similar way as segmentize. The subsampling factor \(J\) in second-order segmentation can be set with the subsample option, while the Nsamples option can be used to set the number of samples \(I\).

D1 = decimate(1:12,'subsample',2)

D1 =
   1     2
   3     4
   5     6
   7     8
   9    10
  11    12

D2 = decimate(1:12,'subsample',4)

D2 =
  1     2     3     4
  5     6     7     8
  9    10    11    12

For \(K\)th-order decimation, the values \(I_2, \ldots, I_K\) are set with Subsample while the values \(I_1, \ldots, I_{K-1}\) are set with Nsamples.

The decimate command supports the other options of segmentize such as Shift, Full and PermToFirst as well. We refer the interested reader to Segmentation.

Dedecimation in Tensorlab¶

The original data can be extracted from the decimated version with the dedecimate command. The command works in a similar way as desegmentize, supporting options such as Dims, Dim, Order, Shift and Method. We refer the interested reader to the discussion of desegmentation in Segmentation.

Computing covariance matrices along different modes¶

The Matlab command cov calculates the variance of a data vector or the covariance matrix of a data matrix, by assuming that each row represents an observation and each column represents a variable. This procedure maps a matrix \(\mat{X}\in\C^{N\times R}\) to a covariance matrix \(\mat{C}\in\C^{R\times R}\).

The dcov command generalizes the functionality of cov to tensor data by allowing the user to indicate the modes along which the tensorization is carried out. The Dims option is an array with two elements, of which the first entry indicates the mode of the observations and of which the second entry indicates the mode of the variables.

When applying dcov on a data tensor \(\ten{X}\in\C^{N\times I \times J}\) with Dims = [1 2], a tensor \(\ten{C}\in\C^{I \times I \times J}\) is returned. \(\ten{C}\) contains covariance matrices along the third mode. When applying dcov on a data tensor \(\ten{X}\in\C^{I\times J \times N \times R}\) with Dims = [3 4], a tensor \(\ten{C}\in\C^{I \times J \times R \times R}\) is returned.

Computing lagged covariance matrices¶

Lagged (or shifted) covariance matrices can be computed with scov. The first argument is the data matrix, and the second argument contains the lags. By definition, the lags should be positive. A third input argument indicates whether the data is prewhitened first. The lagged covariance matrices are stacked along the third mode of the resulting tensor.

X = randn(1000,10);
C = scov(X,[0 1 2]);  % No prewhitening by default

Computing third- and fourth-order cumulants¶

The cum3 and cum4 commands can be used to compute third-order and fourth-order cumulants, respectively. These commands also return the third-order and fourth-order moments as the second output arguments. The first input argument is the data matrix, in which each row represents an observation and in which each column represents a variable.

X = randn(1000,10);
C3 = cum3(X);
size(C3)

ans =
   [10  10  10]

C4 = cum4(X);
size(C4)

ans =
   [10  10  10  10]

A second input argument can be used to indicate whether the data is prewhitened first:

C4 = cum4(X,'prewhitening');

Computing fourth-order cross- and spatio-temporal cumulant¶

Given four data matrices, the xcum4 command computes the fourth-order cross-cumulant. For matrices \(\mat{A}\in\C^{N\times I}\), \(\mat{B}\in\C^{N\times J}\), \(\mat{C}\in\C^{N\times K}\) and \(\mat{D}\in\C^{N\times L}\), the obtained cumulant has size \(I \times J \times K \times L\).

The stcum4 command computes the fourth-order spatio-temporal cumulant. Consider a data matrix \(\mat{X}\in\C^{N\times I}\) and lag parameter \(L\). stcum4 returns a seventh-order tensor of size \(I\times I \times I \times I \times K \times K \times K\) with \(K = 2L+1\).

For a discussion of higher-order statistics, we refer to [42].