Formal defintion

The most common definition of unfolding is the one defined and popularised by the seminal of Kolda and Bader (if you are interested, read their comprehensive paper "Tensor Decompositions and Applications"in SIAM REVIEW, 2009).

We will call that definition of unfolding, that of... Kolda&Bader! They define unfolding as follows:

Let \(\tilde X \in \mathbb{R}^{I_1 \times I_2 \times \cdots \times I_N}\) be a size \((I_1 \times I_2 \times \cdots \times I_N)\) tensor, then element \((i_1 \times i_2 \times \cdots \times i_N)\) of \(\tilde X\) maps to element \((i_n, j)\) of \(\tilde X^{Kolda}_{[n]}\), the n-mode unfolding of \(\tilde X\), with:

Note that if indexing starts at zero, the previous formula simplifies to:

For those used to manipulating arrays, you might recognise here the Fortran ordering of elements!

Example

Note that we index the modes (dimensions) of a tensor from 0, not 1.

Let's take for this example the tensor \(\tilde X\) of size \((3, 4, 2)\) defined by the frontal slices as follows:

and

Then, according to Kolda's definition, its unfolding are given by:

Properties

Given a tensor \(\tilde{X} \in \mathbb{R}^{I_1 \times I_2 \times \cdots \times I_N}\) and its Tucker decomposition \([\!\![~\tilde{G};~\mathbf{U}^{(1)}, \cdots, \mathbf{U}^{(N)}~]\!\!]\), we can express the mode-\(n\) unfolding of \(\tilde{X}\) as:

Notice the reverse order of the kronecker product of the matrix factors, that arises from the Fortran ordering of the elements when vectorising/folding/unfolding a tensor.

Formal definition

In TensorLy, given a tensor \(\tilde X \in \mathbb{R}^{I_1 \times I_2 \times \cdots \times I_N}\), we define the mode-n unfolding of \(\tilde X\) as a matrix \(\mathbf{X}^{TensorLy}_{[n]} \in \mathbb{R}^{I_n, I_M}\), with \(M = \prod_{\substack{k=1,\\k \neq n}}^N I_k\) and is defined by the mapping from element \((i_1, i_2, \cdots, i_N)\) to \((i_n, j)\), with

Example

Using the same tensor \(\tilde X\) as defined previously, the unfoldings according to this new definition are given by:

Properties

One of the nice things with this definition of the unfolding is that it offers nice properties.

Firstly, you might have recognised a C-ordering of the elements: as a result, this formulation leads to much more natural implementations and achieves better performance when using C-ordering of the elements (as numpy does by default). All that is now required is to move the required dimension into the first position and reshape the result into a matrix.

This definition also translates into more natural properties. For instance, given a tensor \(\tilde{X} \in \mathbb{R}^{I_1 \times I_2 \times \cdots \times I_N}\) and its Tucker decomposition \([\!\![~\tilde{G};~\mathbf{U}^{(1)}, \cdots, \mathbf{U}^{(N)}~]\!\!]\), we can express the mode-\(n\) unfolding of \(\tilde{X}\) as:

Example with code, using TensorLy

In TensorLy, we simply define tensors as numpy arrays (which are already multi-dimensional arrays, i.e. tensors):

X = np.array([[[ 0,  1],
               [ 2,  3],
               [ 4,  5],
               [ 6,  7]],

              [[ 8,  9],
               [10, 11],
               [12, 13],
               [14, 15]],

              [[16, 17],
               [18, 19],
               [20, 21],
               [22, 23]]])

You can quickly check the frontal slices as follows:

>>> X[..., 0]
array([[ 0,  2,  4,  6],
       [ 8, 10, 12, 14],
       [16, 18, 20, 22]])

>>> X[..., 1]
array([[ 1,  3,  5,  7],
       [ 9, 11, 13, 15],
       [17, 19, 21, 23]])

To unfold a tensor, simply use the unfold function from TensorLy:

>>> from tensorly import unfold
>>> unfold(X, 0)
array([[ 0,  1,  2,  3,  4,  5,  6,  7],
       [ 8,  9, 10, 11, 12, 13, 14, 15],
       [16, 17, 18, 19, 20, 21, 22, 23]])

Unfolding and ordering of the elements

Consider the tensor X defined above: you might be tempted to get the mode-0 unfolding as defined by Kolda and Bader by simply reshaping it into a matrix of size (3, 4*2). However, this is what happens if you do:

>>> X.reshape((3, -1))
array([[ 1, 13,  4, 16,  7, 19, 10, 22],
       [ 2, 14,  5, 17,  8, 20, 11, 23],
       [ 3, 15,  6, 18,  9, 21, 12, 24]])

What happened here?

Writing X.reshape((3, -1)) is equivalent to writing: X.reshape((3, -1), order='C'). In other words, the elements are read in C-order. Which means that the elements from the last axis are contiguous in memory and read first: you can see this by flattening the array in the default order (C order).

>>> X.ravel() # equivalent to X.ravel(order='C')
array([ 1, 13,  4, 16,  7, 19, 10, 22,  2, 14,  5, 17,  8, 20, 11, 23,  3,
       15,  6, 18,  9, 21, 12, 24])

What you really want is to read the elements column-wise, i.e. in Fortran order: indeed, if you reshape the array in Fortran order, you obtain the mode-0 unfolding:

>>> X.reshape((3, -1), order='F')
array([[ 1,  4,  7, 10, 13, 16, 19, 22],
       [ 2,  5,  8, 11, 14, 17, 20, 23],
       [ 3,  6,  9, 12, 15, 18, 21, 24]])

Unfolding using reshaping in F-order

Based on the above discussion, a simple function to unfold a tensor would be the following:

def f_unfold(tensor, mode=0):
    """Simple unfolding function
        Moves the `mode` axis to the beginning and reshapes in Fortran order
    """
    return np.reshape(np.moveaxis(tensor, mode), (tensor.shape[mode], -1), order='F')

Yes but...

1. we work mostly in C-order
2. moving the axis and reshaping in Fortran order is slow if you manipulate arrays in C-order, which they are by default in numpy which leads us to...
3. always try to use the defaults users expect. This is the Principle of least astonishment.

Unfolding using reshaping in C-order: playing with strides and shape

As we said, what we want is the mode-axis first, then the other dimensions in the correct order. Let's take the example of mode-0 unfolding: when we form X.reshape((mode, 4*2)), the elements of the last mode come first, while the first mode should. (One way to look at it is that the reshaping returns a tensor of size \((I_k, I_n \times \cdots \times I_{k+1} \times I_{k-1} \cdots \times I_1)\) instead of \((I_k, I_1 \times \cdots \times I_{k-1} \times I_{k+1} \cdots \times I_n)\)):

>>> X.reshape((3, -1))
array([[ 1, 13,  4, 16,  7, 19, 10, 22],
       [ 2, 14,  5, 17,  8, 20, 11, 23],
       [ 3, 15,  6, 18,  9, 21, 12, 24]])

To correct this, we need to re-arange the dimensions of the tensor before reshaping.

(0, 1, 2, 3, 4) -> (2, 4, 3, 1, 0)

def reorder(indices, mode):
    """Reorders the elements
    """
    indices = list(indices)
    element = indices.pop(mode)
    return ([element] + indices[::-1])

from numpy.lib.stride_tricks import as_strided

def c_unfold(tensor, mode=0):
    """Returns the mode-`mode` unfolding of `tensor`
        (modes start at `0`)
    """
    shape = list(tensor.shape)
    strides = list(tensor.strides)
    fibers_len = shape[mode]
    shape = reorder(shape, mode)
    strides = reorder(strides, mode)
    return as_strided(tensor, strides=strides, shape=shape).reshape((fibers_len, -1))

A shorter version

The previous code can be written equivalently in a much shorter way using numpy's built-in transpose function that takes care of changing the strides and the shape for us:

def unfold(tensor, mode=0):
    """Returns the mode-`mode` unfolding of `tensor`
    """
    return np.transpose(tensor, reorder(range(tensor.ndim), mode)).reshape((tensor.shape[mode], -1))

Timing all these candidates

This is all well and good but which of these functions is actually faster: we have mentioned this already, but let us see what it is in practice. I have run each method on a tensor of size \((100, 10, 15, 10, 100)\), using the timeit function from the timeit module.

With strides

We can use our knowledge on the ordering of the elements in the rows of the matrix to directly set the strides and the shape. Knowing the size of an element (in bytes), and the number of elements per dimension, we can easily determine the number of bytes to skip to go to the next dimension (strides):

def c_fold(unfolded, mode, shape):
    """Returns the folded tensor of shape `shape` from the `mode`-mode unfolding `unfolded`.
    """
    column_size, element_size = unfolded.strides
    column_dimensions = list(shape)
    # the mode-dimension is contained in the columns
    column_dimensions.pop(mode)
    # the elements of a row are ordered as the first dimension, then the 2nd, etc
    # to read the first, skip 0 elements
    # to read the second dimension, skip n_elements_first_dimension, etc.
    column_dimensions = [1] + np.cumprod(column_dimensions)[:-1].tolist()
    # number of bytes to skip to go from one element to another
    strides = [d * element_size for d in column_dimensions]
    # To go to the next element in dimension mode, we go to the next line of the unfolded tensor
    strides.insert(mode, column_size)
    return as_strided(unfolded, strides=strides, shape=shape)

Concise solution

Again, it is faster and clearer to let numpy do the work. We first reshape the unfolded tensor into a tensor after the transposition was done. Then we undo the transpose:

def fold(unfolded, mode, shape):
    """Returns the folded tensor of shape `shape` from the `mode`-mode unfolding `unfolded`.
    """
    unfolded_indices = reorder(range(len(shape)), mode)
    original_shape = [shape[i] for i in unfolded_indices]
    unfolded = unfolded.reshape(original_shape)

    folded_indices = list(range(len(shape)-1, 0, -1))
    folded_indices.insert(mode, 0)
    return np.transpose(unfolded, folded_indices)