# Decomposition of linear tensor transformations

Claudio Turchetti
DII - Department of Information Engineering
Università Politecnica delle Marche
via Brecce Bianche,12, 60131 Ancona, Italy
<EMAIL_ADDRESS>

###### Abstract

One of the main issues in computing a tensor decomposition is how to choose
the number of rank-one components, since there is no finite algorithms for
determining the rank of a tensor. A commonly used approach for this purpose is
to find a low-dimensional subspace by solving an optimization problem and
assuming the number of components is fixed. However, even though this
algorithm is efficient and easy to implement, it often converges to poor local
minima and suffers from outliers and noise. The aim of this paper is to
develop a mathematical framework for exact tensor decomposition that is able
to represent a tensor as the sum of a finite number of low-rank tensors. In
the paper three different problems will be carried out to derive: i) the
decomposition of a non-negative self-adjoint tensor operator; ii) the
decomposition of a linear tensor transformation; iii) the decomposition of a
generic tensor.

_K_ eywords Machine learning, Computer vision, tensor PCA, tensor
decomposition, self-adjoint operator, tensor basis, eigentensors

##  1 Introduction

Tensors are the higher order generalization of vectors and matrices, and can
consequently be treated as arrays indexed by multiple indices [1],[2], [3].
Thanks to their ability to represent a wide range of real-world data, tensors
have expanded quickly from psychometrics [4], [5] and chemometrics [6], [1] to
image analysis [7], [8], [9], [10], [11], [12], [13], big data representation
[14], machine learning, [15], [16], [17], [18], [19], [20], sensor array
processing [21], and much more [22], [23], [24], [22],[25], [26]. As the order
of tensor increases, the number of entries in the array increases
exponentially thus involving prohibitively computational and storage costs. To
prevent these limitations, techniques for dimensionality reduction are
required to make possible the application of tensors to real data. A widely
adopted tecnique for this purpose is to represent high-dimensional tensors as
linear combination of low-dimensional tensors. If $\tensor{\mathcal{A}}{}$ is
an order-$d$ tensor, to reach this goal is essential to find low-order tensors
$\tensor{\mathcal{U}}{}_{i}$ so that

$\tensor{\mathcal{A}}{}=\sum_{{i}=1}^{{r}}\sigma_{i}\tensor{\mathcal{U}}{}_{i},$
(1)

where $\sigma_{i}$ are positive scalars. In this way a low-dimensional
representation can be derived by simply retaining the principal components of
such a decomposition, that is the ones corresponding to the highest values of
$\sigma_{i}$. In this context several very powerful tensor decomposition
approaches have been developed in the past [27], [28], [29], [30], in which
tensors $\tensor{\mathcal{U}}{}_{i}$ in the summation (1) are rank-one
tensors. Among these, Tucker decomposition [31] and CANDECOMP/PARAFAC (CP)
decomposition [32], [5] are the most popular and fundamental models. Tucker
model decomposes a tensor into a core tensor multiplied by a factor matrix
along each mode, while CP model factorizes a tensor into a weighted sum of
rank-1 tensors. Following these two seminal approaches, several extensions for
low rank tensor decomposition have been proposed in the past few decades. Some
of the most relevant developments in this field are: multilinear principal
component analysis (MPCA) [33], [34], [35], [36], [37], [38], [39], [40],
[41], [42], tensor rank-one decomposition (TROD) [43], [44], [32], [45], [46]
, hierarchical Tucker PCA (HT-PCA) [47], [48], [49] and tensor-train PCA (TT-
PCA) [50], [51].
One of the main issues in computing a tensor decomposition is how to choose
the number of rank-one components, since there is no finite algorithms for
determining the rank of a tensor [45], [52], defined as the smallest number of
rank-one tensors that generate $\tensor{\mathcal{A}}{}$ as their sum. As a
consequence, all the above methods for tensor decomposition try to find a low-
dimensional subspace by solving an optimization problem and assuming the
number of components is fixed. A commonly used approach for this purpose is
the alternating least square (ALS) method [53], which assumes the solution of
all but one mode is known and then estimating the unknown parameter set of the
remaining mode. However, even though this algorithm is efficient and easy to
implement, it often converges to poor local minima and suffers from outliers
and noise [54].
The aim of this paper is to develop a mathematical framework for exact tensor
decomposition that is able to represent a tensor as the sum of a finite number
of low-rank tensors. The core of the proposed approach is the derivation of a
decomposition for non-negative self-adjoint tensor operators. In particular it
will be proven that a correspondence exists between the spectral decomposition
of a symmetric matrix, a result that is the basis of PCA and SVD development
in matrix analysis, and the decomposition of a self-adjoint tensor operator.
In the paper three different problems will be carried out to derive: i) the
decomposition of a non-negative self-adjoint tensor operator; ii) the
decomposition of a linear tensor transformation; iii) the decomposition of a
generic tensor. In particular, with reference to the first issue, the
properties of a self-adjoint tensor operator will be studied and the
equivalence between eigenvalue equation for this operator and standard matrix
eigenvalue equation will be proven.
The paper is organized as follows. Section 2 deals with the finite linear
space of tensors and some fundamental concepts on linear tensor
transformations. In Section 3 the mathematical framework for exact tensor
decomposition is developed. Section 4 reports some numerical results to
validate the mathematical framework previously presented.

##  2 The linear space of tensors

Let us refer to the set $\tensor{\mathcal{L}}{}_{d}$ of real order-$d$ tensors
$\tensor{\mathcal{X}}{}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times{I_{d}}}$ whose generic element will be denoted as

$\tensor{\mathcal{X}}{}_{\textbf{i}}=\tensor{\mathcal{X}}{}_{i_{1},\ldots,i_{d}},\;\;\textbf{i}={(i_{1},\ldots,i_{d})}\quad.$
(2)

The vector i, in bold font, is a subscript vector where the index $i_{k}$ in
the $k$th mode ranges in the interval $1\leq{i_{k}}\leq{I_{k}},\hskip
5.69046ptk=1:d$. As $\tensor{\mathcal{L}}{}_{d}$ is closed under addition and
multiplication by scalar, this set forms a linear space, called the tensor
space. In the following we will refer to real tensors alone, thus for
simplicity the term real will be omitted being always implied.

#### Inner product

The inner product of two tensors of the same size
$\tensor{\mathcal{X}}{}_{\textbf{i}},\tensor{\mathcal{Y}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots I_{d}}$ is a real scalar defined as

$\left\langle\tensor{\mathcal{X}}{},\tensor{\mathcal{Y}}{}\right\rangle=\tensor{\mathcal{X}}{}_{\hskip
0.56917pt\textbf{i}}\tensor{\mathcal{Y}}{}_{\hskip 0.56917pt\textbf{i}}.$ (3)

Here the Einstein summation convention, that interprets repeated subscript as
summation over that index, has been used. Following this convention, (3) is
equivalent to

$\left\langle\tensor{\mathcal{X}}{},\tensor{\mathcal{Y}}{}\right\rangle=\sum_{{i_{1}}=1}^{{I_{1}}}...\sum_{{i_{d}}=1}^{{I_{d}}}\tensor{\mathcal{X}}{}_{i_{1},\ldots,i_{d}}\tensor{\mathcal{Y}}{}_{i_{1},\ldots,i_{d}}$
(4)

In the paper we will extensively use the Einstein convention to simplify
mathematical notation. From this definition
$\tensor{\mathcal{X}}{},\tensor{\mathcal{Y}}{}$ are said to be orthogonal if
$\left\langle\tensor{\mathcal{X}}{},\tensor{\mathcal{Y}}{}\right\rangle=0$,
and it follows that the norm of a tensor is given by

$\|{\tensor{\mathcal{X}}{}}\|=\sqrt{\left\langle\tensor{\mathcal{X}}{},\tensor{\mathcal{X}}{}\right\rangle}$
(5)

#### Subspace

A subspace $S$ of $\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times{I_{d}}}$ is
a subset that is also a tensor space. Given a collection of tensors
$\tensor{\mathcal{U}}{}_{1},\dots\tensor{\mathcal{U}}{}_{n}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times{I_{d}}}$, the set of all linear combinations of these
vectors is a subspace referred to as the
$span\\{\tensor{\mathcal{U}}{}_{1},\dots\tensor{\mathcal{U}}{}_{n}\\}$:

$span\\{\tensor{\mathcal{U}}{}_{1},\dots\tensor{\mathcal{U}}{}_{n}\\}=\\{\alpha_{i}\tensor{\mathcal{U}}{}_{i}:\hskip
5.69046pt\alpha_{i}\in\mathbb{R},\hskip 5.69046pti=1:n\\}.$ (6)

#### Contraction

Given two tensors

$\tensor{\mathcal{X}}{}_{\hskip 0.56917pt\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times{I_{d}}},\hskip
14.22636pt\tensor{\mathcal{Y}}{}_{\hskip
0.56917pt\textbf{j}}\in\mathbb{R}^{J_{1}\times
J_{2}\times\ldots\times{J_{e}}}$ (7)

and assuming a common vector index k exists such that i, j can be partitioned
as

$\textbf{i}=(\textbf{l},\textbf{k},\textbf{m}),\hskip
14.22636pt\textbf{j}=(\textbf{p},\textbf{k},\textbf{q})$ (8)

the tensor product of $\tensor{\mathcal{X}}{}$ and $\tensor{\mathcal{Y}}{}$
along the multi-index k combines the two tensors to give a third tensor
$\tensor{\mathcal{Z}}{}$, whose generic element is given by

$\tensor{\mathcal{Z}}{}_{\hskip
0.56917pt\textbf{l},\textbf{m},\textbf{p},\textbf{q}}=\tensor{\mathcal{X}}{}_{\hskip
0.56917pt\textbf{l},\textbf{k},\textbf{m}}\tensor{\mathcal{Y}}{}_{\hskip
0.56917pt\textbf{p},\textbf{k},\textbf{q}}$ (9)

that implies summation over the vector index k. We refer to this operation as
contraction of index k.

#### Outer product

The outer product of tensor
$\tensor{\mathcal{X}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{p}}$ with tensor
$\tensor{\mathcal{Y}}{}_{\textbf{j}}\in\mathbb{R}^{J_{1}\times
J_{2}\times\ldots\times J_{q}}$ is the order-$(p+q)$ tensor
$\tensor{\mathcal{Z}}{}$ defined as

$\tensor{\mathcal{Z}}{}=\tensor{\mathcal{X}}{}\circ\tensor{\mathcal{Y}}{},$
(10)

where the generic element of $\tensor{\mathcal{Z}}{}$ is

$\tensor{\mathcal{Z}}{}_{\textbf{i},\textbf{j}}=\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{Y}}{}_{\textbf{j}},\hskip
8.5359pt\textbf{i}={(i_{1},\ldots,i_{p})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{q})}$ (11)

In particular for two vectors $x,y\in{R^{I_{1}}}$ the outer product is the
matrix ${Z}=x\circ y$ whose generic element is

$z_{ij}=x_{i}\,y_{j}\quad.$ (12)

#### Vector-to-linear index transformation

With reference to an order-$(p+q)$ tensor

${\tensor{\mathcal{X}}{}}_{\textbf{j},\textbf{i}},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{p})},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{q})}$ (13)

it is worth to notice that the indices $i_{1},\ldots,i_{q}$ can be arranged in
$M=I_{1}I_{2}\dots I_{q}$ different ways (the number of permutations with
repetition), so that a correspondence $\alpha$

$m=\alpha(\textbf{i}),\hskip 5.69046pt1\leq{m}\leq{M},\hskip
5.69046pt1\leq{i_{k}}\leq{I_{k}},\hskip 5.69046ptk=1:q$ (14)

between $m$ and i holds.

The inverse transformation exists, since the correspondence is one-to-one, and
is denoted by

$\textbf{i}=\alpha^{-1}(m),\hskip 5.69046pt1\leq{m}\leq{M},\hskip
14.22636pt1\leq{i_{k}}\leq{I_{k}},\hskip 5.69046ptk=1:q$ (15)

By defining the $(M\times q)$ matrix $T$ whose $m$th row represents the
corresponding vector index $\textbf{i}(m)$, then the inverse transformation
(15) is formally given by

$\textbf{i}=T(m,:),\hskip 5.69046pt1\leq{m}\leq M$ (16)

As a result this operation, transforming the vector index i to the linear
index $m$, turns the order-$(p+q)$ tensor to the order-$(p+1)$ tensor

${\tensor{\mathcal{X}}{}}_{\textbf{j},{m}},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{p})},\hskip 5.69046ptm=1:M.$ (17)

The vector-to-linear-index transformation can be applied to more than one
index. As an example, by transforming the two vector indices of
$\tensor{\mathcal{X}}{}_{\textbf{j},\textbf{i}}$ gives rise to the matrix
$\tensor{\mathcal{X}}{}_{\textbf{j}(n),\textbf{i}(m)}=X_{n,m},\hskip
5.69046ptn=1:N,m=1:M$, with $N=J_{1}J_{2}\dots J_{p}$, $M=I_{1}I_{2}\dots
I_{q}$.

#### Linear transformations

A two-index tensor $\tensor{\mathcal{A}}{}_{\hskip
0.56917pt\textbf{i},\textbf{j}}\in{R^{I\times{J}}}$, where
$I=I_{1}\times\ldots\times{I_{d}}$, $J=J_{1}\times\ldots\times{J_{e}}$,
defines a linear transformation
$\tensor{\mathcal{A}}{}:\mathbb{R}^{J}\rightarrow{\mathbb{R}^{I}}$ from
$\mathbb{R}^{J}$ to $\mathbb{R}^{I}$ as the product

$\tensor{\mathcal{Y}}{}_{\hskip
0.56917pt\textbf{i}}=\tensor{\mathcal{A}}{}_{\hskip
0.56917pt\textbf{i},\textbf{j}}\tensor{\mathcal{X}}{}_{\hskip
0.56917pt\textbf{j}}$ (18)

that transforms a tensor $\tensor{\mathcal{X}}{}_{\hskip
0.56917pt\textbf{j}}\in{R^{J}}$ to the tensor $\tensor{\mathcal{Y}}{}_{\hskip
0.56917pt\textbf{i}}\in{R^{I}}$. Among the possible transformations a
particular role is assumed by transformations from $\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{d}}$ to itself. In this particular case the tensors
$\tensor{\mathcal{Y}}{}_{\hskip 0.56917pt\textbf{i}}$ and
$\tensor{\mathcal{X}}{}_{\hskip 0.56917pt\textbf{j}}$ are both in
$\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{d}}$ and the order-$2d$
tensor $\tensor{\mathcal{A}}{}_{\hskip 0.56917pt\textbf{i},\textbf{j}}$ is
said to be a tensor operator (or endomorphism). Thus, the operator
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}$ establishes a transformation from
$\tensor{\mathcal{L}}{}_{d}$ to $\tensor{\mathcal{L}}{}_{d}$ and (18) can be
rewritten without ambiguity as

$\tensor{\mathcal{Y}}{}=\tensor{\mathcal{A}}{}\tensor{\mathcal{X}}{}$ (19)

A tensor $\tensor{\mathcal{V}}{}$ in $\tensor{\mathcal{L}}{}_{d}$ is said to
be an eigentensor if $\tensor{\mathcal{V}}{}\neq{0}$ and if for some scalar
$\lambda$ the following equation

$\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}=\lambda\tensor{\mathcal{V}}{}$
(20)

is satisfied. The scalar $\lambda$ is known as the eigenvalue of
$\tensor{\mathcal{A}}{}$ associated with the eigentensor
$\tensor{\mathcal{V}}{}$,

#### Self-adjoint non-negative definite operators

A self-adjoint operator is an operator for which the following property

$\left\langle\tensor{\mathcal{A}}{}\tensor{\mathcal{Y}}{},\tensor{\mathcal{Z}}{}\right\rangle=\left\langle\tensor{\mathcal{Y}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{Z}}{}\right\rangle$
(21)

holds.

As an example, an order-$2d$ tensor

${\tensor{\mathcal{A}}{}}_{\hskip 0.56917pt\textbf{i},\textbf{j}},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{d})}$ (22)

such that

${\tensor{\mathcal{A}}{}}_{\hskip
0.56917pt\textbf{i},\textbf{j}}={\tensor{\mathcal{A}}{}}_{\hskip
0.56917pt\textbf{j},\textbf{i}}$ (23)

is a self-adjoint operator, since it results

$\left\langle\tensor{\mathcal{A}}{}\tensor{\mathcal{Y}}{},\tensor{\mathcal{Z}}{}\right\rangle={\tensor{\mathcal{A}}{}}_{\hskip
0.56917pt\textbf{i},\textbf{j}}{\tensor{\mathcal{Y}}{}}_{\hskip
0.56917pt\textbf{j}}{\tensor{\mathcal{Z}}{}}_{\hskip
0.56917pt\textbf{i}}={\tensor{\mathcal{Y}}{}}_{\hskip
0.56917pt\textbf{j}}{\tensor{\mathcal{A}}{}}_{\hskip
0.56917pt\textbf{j},\textbf{i}}{\tensor{\mathcal{Z}}{}}_{\hskip
0.56917pt\textbf{i}}=\left\langle\tensor{\mathcal{Y}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{Z}}{}\right\rangle.$
(24)

In the particular case the operator acts on the order-1 tensor space
$R^{I_{1}}$, that correspoinds to $d=1$ in (22), the self-adjoint operator
(23) reduces to a symmetric matrix $A_{i_{1},j_{1}}=A_{j_{1},i_{1}}$.
A non-negative definite operator is such that the following property

$\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}\right\rangle\geq
0$ (25)

holds for every tensor $\tensor{\mathcal{V}}{}$.

An important class of operators is the class of self-adjoint non-negative
definite (SA-NND) tensor operators that are both self-adjoint and non-
negative. For such tensors, on the basis of the properties stated before, it
follows that eigenvalues of an (SA-NND) tensor are non-negative and the
eigentensors belonging to distinct eigenvalues are orthogonal.

Indeed, combining (20) and (21) yields the scalar

$\lambda=\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}\right\rangle/{\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{V}}{}\right\rangle}$
(26)

that is non-negative due to property (25).

As for the eigentensor orthogonality, suppose that
$\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}_{1}=\lambda_{1}\tensor{\mathcal{V}}{}_{1}$
and
$\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}_{2}=\lambda_{2}\tensor{\mathcal{V}}{}_{2}$
for $\lambda_{1}\neq\lambda_{2}$, if $\tensor{\mathcal{A}}{}$ is self-adjoint
, then

$\left\langle\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{V}}{}_{2}\right\rangle=\lambda_{1}\left\langle\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{V}}{}_{2}\right\rangle$
(27)

and also

$\left\langle\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{V}}{}_{2}\right\rangle=\left\langle\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}_{2}\right\rangle=\lambda_{2}\left\langle\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{V}}{}_{2}\right\rangle$
(28)

Therefore,
$(\lambda_{1}-\lambda_{2})\left\langle\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{V}}{}_{2}\right\rangle=0$,
and hence
$\left\langle\tensor{\mathcal{V}}{}_{1},\tensor{\mathcal{V}}{}_{2}\right\rangle=0$,
since $\lambda_{1}\neq\lambda_{2}$.

An SA-NND operator can be easily derived from a tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}$ as follows

$\tensor{\mathcal{G}}{}_{\textbf{i,j}}=\tensor{\mathcal{A}}{}_{\textbf{k,i}}\tensor{\mathcal{A}}{}_{\textbf{k,j}}$
(29)

The operator $\tensor{\mathcal{G}}{}_{\textbf{i,j}}$ is self-adjoint, in fact
we have

$\tensor{\mathcal{G}}{}_{\textbf{i,j}}=\tensor{\mathcal{A}}{}_{\textbf{k,i}}\tensor{\mathcal{A}}{}_{\textbf{k,j}}=\tensor{\mathcal{A}}{}_{\textbf{k,j}}\tensor{\mathcal{A}}{}_{\textbf{k,i}}=\tensor{\mathcal{G}}{}_{\textbf{j,i}}$
(30)

The eigenvalues of $\tensor{\mathcal{G}}{}_{\textbf{i,j}}$ are non-negative,
since

$\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{G}}{}\tensor{\mathcal{V}}{}\right\rangle=\tensor{\mathcal{V}}{}_{\textbf{i}}\tensor{\mathcal{G}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j}}=\tensor{\mathcal{A}}{}_{\textbf{k,i}}\tensor{\mathcal{V}}{}_{\textbf{i}}(\tensor{\mathcal{A}}{}_{\textbf{k,j}}\tensor{\mathcal{V}}{}_{\textbf{j}})=\tensor{\mathcal{Z}}{}_{\textbf{k}}\tensor{\mathcal{Z}}{}_{\textbf{k}}\geq
0$ (31)

where
$\tensor{\mathcal{Z}}{}_{\textbf{k}}=\tensor{\mathcal{A}}{}_{\textbf{k,i}}\tensor{\mathcal{V}}{}_{\textbf{i}}$,
so that from (26) it results $\lambda\geq 0$. The operator
$\tensor{\mathcal{G^{\prime}}}{}_{\textbf{i,j}}=\tensor{\mathcal{A}}{}_{\textbf{i,k}}\tensor{\mathcal{A}}{}_{\textbf{j,k}}$
satisfies this property as well. In particular for linear indices
$\textbf{i}=i$, $\textbf{j}=j$, $\textbf{k}=k$, the operator reduces to the
symmetric matrix

${G}_{i,j}={A}_{k,i}{A}_{k,j}=({A}^{T}_{i,k}){A}_{k,j}$ (32)

which can be written as $G=A^{T}A$. Similarly, the operator
$\tensor{\mathcal{G^{\prime}}}{}_{\textbf{i,j}}$ reduces to the matrix
$G^{T}=AA^{T}$. The matrix $G$ represents the well known Gram-matrix, while
matrix $G^{T}$ its transpose. Thus, tensor
$\tensor{\mathcal{G}}{}_{\textbf{i,j}}$ is the analog in tensor space of Gram-
matrix $G$. Assuming $G$ has size $(L\times L)$ and rank $r$, i.e.
$r=rank(G)$, it is well known from geometry that, since $G$ is a symmetric
non-negative definite matrix, it admits $L$ eigenvalues $\\{\lambda_{i},\hskip
5.69046pti=1:L\\}$ such that
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{r}>{0}=\lambda_{r+1}=\dots=\lambda_{L}$.
Thus $G$ can be diagonalized as

$G_{n,m}U_{m,p}=U_{n,q}\Lambda_{q,p},\hskip 5.69046ptq,p=1:r$ (33)

where $U=(u_{1},\dots,u_{r})$ is the matrix whose columns are the eigenvectors
corresponding to non-zero eigenvalues in the matrix
$\Lambda=diag(\lambda_{1},\dots,\lambda_{r})$. Using orthogonality of matrix
$U$, i.e. $U_{m,p}U_{m,p^{\prime}}=\delta_{p,p^{\prime}}$ and
$\Lambda_{q,p}=\lambda_{p}\delta_{q,p}$, (33) becomes

$G_{n,m}=\lambda_{p}U_{n,p}U_{m,p}.$ (34)

This relationship is called spectral decomposition of matrix $G$ and
represents a fundamental result in matrix analysis, since both principal
component analysis (PCA) and singular value decomposition (SVD) are based on
this decomposition. Thus it is essential to ask whether a similar relationship
can be derived for the class of (SA-NND) tensors.

##  3 Tensor decomposition

Having proven there is a correspondence between self-adjoint tensor operators
and symmetric matrices, the aim of this section is to discover that a
correspondence to the spectral decomposition (34) exists in tensor space. To
take advantage of these considerations, tensor $\tensor{\mathcal{A}}{}$ to be
decomposed will be treated as a linear transformation from tensor space
$\mathbb{R}^{J}$ to tensor space $\mathbb{R}^{I}$, that is
$\tensor{\mathcal{A}}{}:\mathbb{R}^{J}\rightarrow{\mathbb{R}^{I}}$, which for
linear indices reduces to the matrix $A$. On the basis of this assumption, two
cases will be considered first: i) the tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times I}$ is real self-
adjoint operator and ii) the tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times J}$, with $I\neq
J$, is a linear transformation. Using these results, then the most general
case of a generic tensor will be treated at the end of the section.

### A. Decomposition of a self-adjoint non-negative tensor operator

The first case will be studied here deals with an order-$2d$ tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times I}$ that
represents a linear operator (or endomorphism)
$\tensor{\mathcal{A}}{}:\mathbb{R}^{I}\rightarrow{\mathbb{R}^{I}}$. Thus, with
reference to such a tensor, we state the following proposition.

#### Proposition1

Let

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times I},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{d})}$ (35)

be a non-negative self-adjoint order-$2d$ tensor operator with ${I=I_{1}\times
I_{2}\times\ldots\times I_{d}}$, there exists a set of scalars
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{r}>{0}$, and a set of
orthogonal tensors
$\tensor{\mathcal{U}}{}=\\{\tensor{\mathcal{U}}{}_{\textbf{i},p},\hskip
5.69046ptp=1:r\\}$, such that the following decomposition

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\lambda_{p}{\tensor{\mathcal{U}}{}_{\textbf{i},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p},\hskip
5.69046ptp=1:r$ (36)

holds.

#### Proof

First, we want to prove that the following eigenvalue equation

$\tensor{\mathcal{A}}{}\tensor{\mathcal{U}}{}=\lambda\tensor{\mathcal{U}}{},$
(37)

where $\tensor{\mathcal{A}}{}$ is a self-adjoint operator, admits solutions
with $\tensor{\mathcal{U}}{}\in\mathbb{R}^{I}$ and $\lambda>{0}$. Using index
convention introduced in Section 2, Eq. (37) is equivalent to the $L$
equations

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j}}=\lambda\tensor{\mathcal{U}}{}_{\textbf{i}},\hskip
8.5359pt\textbf{i}\in{I},\hskip 5.69046ptL=dim(I)=I_{1}I_{2}\dots I_{d}$ (38)

which can be written as

$(\tensor{\mathcal{A}}{}_{\textbf{i,j}}-\lambda\delta_{\textbf{i,j}})\tensor{\mathcal{U}}{}_{\textbf{j}}=0,$
(39)

where the term on the left
$\tensor{\mathcal{Y}}{}_{\textbf{i}}=(\tensor{\mathcal{A}}{}_{\textbf{i,j}}-\lambda\delta_{\textbf{i,j}})\tensor{\mathcal{U}}{}_{\textbf{j}}$
represents an order-d tensor. The solutions of (39) are independent on how the
terms in the equation are arranged as entries in the tensor
$\tensor{\mathcal{Y}}{}_{\textbf{i}}$, meaning that (39) is invariant to a
one-to-one transformation of subscript vectors. Assuming indices
$\textbf{i},\textbf{j}$ are obtained by the inverse vector-to-linear index
transformation (18) then, substituting the vector indices
$\textbf{i}(n)=T(n,:)$, $\textbf{j}(m)=T(m,:)$, (39) becomes

$(\tensor{\mathcal{A}}{}_{\textbf{i}(n),\textbf{j}(m)}-\lambda\delta_{\textbf{i}(n),\textbf{j}(m)})\tensor{\mathcal{U}}{}_{\textbf{j}(m)}=0$
(40)

By defining the $(L\times{L})$ matrices

$a_{n,m}=\tensor{\mathcal{A}}{}_{\textbf{i}(n),\textbf{j}(m)},\hskip
5.69046pt\delta_{n,m}=\delta_{\textbf{i}(n),\textbf{j}(m)}$ (41)

and the vector

$u_{m}=\tensor{\mathcal{U}}{}_{\textbf{j}(m)},$ (42)

thus (40) can be re-arranged as

$(a_{n,m}-\lambda\delta_{n,m})u_{m}=0$ (43)

or equivalently

$a_{n,m}u_{m}=\lambda{u_{n}},\hskip 5.69046ptm,n=1:L$ (44)

Equation (44) is equivalent to

$Au=\lambda{u}.$ (45)

where the matrix $[A]_{n,m}=a_{n,m}$ is a symmetric non-negative definite
($L\times{L}$) matrix. In fact, as $\tensor{\mathcal{A}}{}_{\textbf{i,j}}$ is
self-adjoint, then it results

$a_{n,m}=\tensor{\mathcal{A}}{}_{\textbf{i}(n),\textbf{j}(m)}=\tensor{\mathcal{A}}{}_{\textbf{j}(m),\textbf{i}(n)}=a_{m,n}$
(46)

that proves matrix $a_{n,m}$ is symmetric (or self-adjoint). In addition,
since
$\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}\right\rangle\geq
0$ for every tensor $\tensor{\mathcal{V}}{}$, it results

$0\leq\tensor{\mathcal{V}}{}_{\textbf{j}}\tensor{\mathcal{A}}{}_{\textbf{i},\textbf{j}}\tensor{\mathcal{V}}{}_{\textbf{j}}=\tensor{\mathcal{V}}{}_{\textbf{i}(m)}\tensor{\mathcal{A}}{}_{\textbf{i}(m),\textbf{j}(n)}\tensor{\mathcal{V}}{}_{\textbf{j}(n)}=v_{m}{A_{m,n}}v_{n}$
(47)

for every vector $[v]_{n}=v_{n}$, thus proving the matrix $A$ is non-negative.

Assuming $r=rank(A)$, a set of scalars
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{r}>{0}$ exist such that (45)
holds, thus (45) can be rewritten as

$A_{n,m}U_{m,p}=U_{n,q}\Lambda_{q,p}$ (48)

where $U=(u_{1},\dots,u_{r})$ and
$\Lambda=diag(\lambda_{1},\dots,\lambda_{r})$. Using the inverse vector-to-
linear index transformation (16) we have

$U_{n(\textbf{i}),q}=\tensor{\mathcal{U}}{}_{\textbf{i},q},\hskip
8.5359ptA_{n(\textbf{i}),m(\textbf{j})}=\tensor{\mathcal{A}}{}_{\textbf{i,j}}$
(49)

and (48) becomes

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p}=\tensor{\mathcal{U}}{}_{\textbf{i},q}\Lambda_{q,p},$
(50)

or that is the same

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}{\tensor{\mathcal{U}}{}_{\textbf{j},p}}=\lambda_{p}{\tensor{\mathcal{U}}{}_{\textbf{i},p}},\hskip
5.69046ptp=1:r$ (51)

This proves that (37) admits solutions with $\lambda>0$. To decompose
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}$, the tensor product of (51) and
$\tensor{\mathcal{U}}{}_{\textbf{j},p^{\prime}}$ yields

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}{\tensor{\mathcal{U}}{}_{\textbf{j},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p^{\prime}}=\lambda_{p}{\tensor{\mathcal{U}}{}_{\textbf{i},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p^{\prime}},\hskip
5.69046ptp=1:r$ (52)

and from orthogonality
${\tensor{\mathcal{U}}{}_{\textbf{j},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p^{\prime}}=\delta_{p,p^{\prime}}$,
it results

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\delta_{p,p^{\prime}}=\lambda_{l}{\tensor{\mathcal{U}}{}_{\textbf{i},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p^{\prime}},\hskip
5.69046ptp=1:r$ (53)

that is

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\lambda_{p}{\tensor{\mathcal{U}}{}_{\textbf{i},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p},\hskip
5.69046ptp=1:r$ (54)

and this concludes the proof.

A pseudo-code of the procedure previously described to derive a decomposition
of a self-adjoint non-negative definite operator is described in Algorithm 1.

Algorithm 1 Decomposition of self-adjoint non-negative definite tensor
operator

INPUT: order-$2d$ self-adjoint non-negative definite operator
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times I},L=dim(I)$

$\tensor{\mathcal{A}}{}_{\textbf{i,j}},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{d})}$

1\. Derive ($L\times{d}$) matrix $T$ defining vector-to-linear transformation

$\textbf{i}(n)=T(n,:),\hskip 5.69046pt\textbf{j}(m)=T(m,:),\hskip
5.69046ptm,n=1:L$

2\. Compute symmetric non-negative definite $(L\times L)$ matrix $A$

$A_{n,m}=\tensor{\mathcal{A}}{}_{\textbf{i}(n),\textbf{j}(m)},\hskip
5.69046ptm,n=1:L$

3\. Solve eigenvalue equation

$A_{n,m}U_{m,p}=U_{n,q}\Lambda_{q,p},\hskip 5.69046ptq,p=1:r,\hskip
5.69046pt\Lambda=diag(\lambda_{1},\dots,\lambda_{r}),\hskip
5.69046ptr=rank(A)$,

4\. Convert matrix $U_{n,q}$ to tensor $\tensor{\mathcal{U}}{}_{\textbf{i},q}$
by inverse vector-to-linear transformation

$U_{n(\textbf{i}),q}=\tensor{\mathcal{U}}{}_{\textbf{i},q}$

OUTPUT: the decomposition

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\lambda_{p}{\tensor{\mathcal{U}}{}_{\textbf{i},p}}\tensor{\mathcal{U}}{}_{\textbf{j},p},\hskip
5.69046ptp=1:r$

### B. Decomposition of a linear tensor transformation

The second case refers to the order-$(d+e)$ tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times J}$ that
represents a linear transformation
$\tensor{\mathcal{A}}{}:\mathbb{R}^{J}\rightarrow{\mathbb{R}^{I}}$. Thus, we
state the following proposition.

#### Proposition2

Given

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times J},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{e})}$ (55)

where ${I=I_{1}\times I_{2}\times\ldots\times I_{d}}$ and ${J=J_{1}\times
J_{2}\times\ldots\times J_{e}}$, there exist sets of orthogonal tensors
$\tensor{\mathcal{U}}{}=\\{\tensor{\mathcal{U}}{}_{\textbf{i},p}\in\mathbb{R}^{I},\hskip
5.69046ptp=1:r\\}$,
$\tensor{\mathcal{V}}{}=\\{\tensor{\mathcal{V}}{}_{\textbf{j},p}\in\mathbb{R}^{J},\hskip
5.69046ptp=1:r\\}$, and a set of scalars
$\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{r}>{0}$, such that the following
decomposition

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\sigma_{p}\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{V}}{}_{\textbf{j},p}\hskip
5.69046pt,p=1:r.$ (56)

holds.

#### Proof

The tensor

$\tensor{\mathcal{G}}{}_{\textbf{l,j}}=\tensor{\mathcal{A}}{}_{\textbf{i,l}}\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{J\times
J},\hskip 5.69046pt\textbf{l}=(l_{1},\dots,l_{d})$ (57)

is an SA-NND operator, thus the Proposition1 applies and the eigenvalue
equation

$\tensor{\mathcal{G}}{}_{\textbf{l,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}=\tensor{\mathcal{V}}{}_{\textbf{l},q}\Lambda_{q,p},\hskip
5.69046ptq,p=1:r$ (58)

can be derived from (36), where
$\Lambda_{q,p}=diag(\lambda_{1},\dots,\lambda_{r})$, and
$\tensor{\mathcal{V}}{}=\\{\tensor{\mathcal{V}}{}_{\textbf{j},p},\hskip
5.69046ptp=1:r\\}$ is a set of orthogonal tensors in $\mathbb{R}^{J}$. Since
$\lambda_{i}>0$ we can pose $\Lambda=S^{2}$ with
$S_{q,p}=diag(\sigma_{1},\dots,\sigma_{r})$, then multiplying each side of
(58) by $\tensor{\mathcal{V}}{}_{\textbf{l},q^{\prime}}$ yields

$\tensor{\mathcal{V}}{}_{\textbf{l},q^{\prime}}\tensor{\mathcal{G}}{}_{\textbf{l,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}=\tensor{\mathcal{V}}{}_{\textbf{l},q^{\prime}}\tensor{\mathcal{V}}{}_{\textbf{l},q}S^{2}_{q,p}$
(59)

and by orthogonality of $\tensor{\mathcal{V}}{}_{\textbf{l},q}$, that is
$\tensor{\mathcal{V}}{}_{\textbf{l},q^{\prime}}\tensor{\mathcal{V}}{}_{\textbf{l},q}=\delta_{q^{\prime},q}$,
we have

$\tensor{\mathcal{V}}{}_{\textbf{l},q}\tensor{\mathcal{A}}{}_{\textbf{i,l}}\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}=S^{2}_{q,p}$
(60)

.

The matrix $S_{q,p}$ is diagonal full rank, thus $S^{-1}_{q,p}$ exists and
(60) can be rewritten as

$\tensor{\mathcal{A}}{}_{\textbf{i,l}}\tensor{\mathcal{V}}{}_{\textbf{l},q}S^{-1}_{q,p}(\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}S^{-1}_{p,q})=\delta_{p,q}$
(61)

which represents the orthogonality condition for the tensors

$\tensor{\mathcal{U}}{}_{\textbf{i},q}=\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}S^{-1}_{p,q},\hskip
5.69046ptq=1:r.$ (62)

Thus

$\tensor{\mathcal{U}}{}=\\{\tensor{\mathcal{U}}{}_{\textbf{i},q},\hskip
5.69046ptq=1:r\\}$ (63)

is a set of orthogonal tensors in $\mathbb{R}^{I}$. The inner product of two
generic terms $\tensor{\mathcal{U}}{}_{\textbf{i},s}$ and
$\tensor{\mathcal{U}}{}_{\textbf{i},q}$ gives

$\tensor{\mathcal{U}}{}_{\textbf{i},s}\tensor{\mathcal{U}}{}_{\textbf{i},q}=\tensor{\mathcal{U}}{}_{\textbf{i},s}\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}S^{-1}_{p,q}$
(64)

and for orthogonality, that is
$\tensor{\mathcal{U}}{}_{\textbf{i},s}\tensor{\mathcal{U}}{}_{\textbf{i},q}=\delta_{s,q}$,
it results

$S_{q,p}=\tensor{\mathcal{U}}{}_{\textbf{i},q}\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}$
(65)

Using $S_{q,p}=\sigma_{p}\delta_{q,p}$ and the orthogonality of
$\tensor{\mathcal{U}}{}_{\textbf{i},q}$ and
$\tensor{\mathcal{V}}{}_{\textbf{j},p}$, (65) is equivalent to

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\sigma_{p}\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{V}}{}_{\textbf{j},p},\hskip
5.69046ptp=1:r$ (66)

and this concludes the proof.

A pseudo-code of the procedure previously described to derive a decomposition
of a linear tensor transformation is described in Algorithm 2.

Algorithm 2 Decomposition of a linear tensor transformation

INPUT: order-$(d+e)$ tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times J}$

$\tensor{\mathcal{A}}{}_{\textbf{i,j}},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{e})}$

1\. Derive non-negative self-adjoint operator

$\tensor{\mathcal{G}}{}_{\textbf{l,j}}=\tensor{\mathcal{A}}{}_{\textbf{i,l}}\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{J\times
J}$

2\. Solve the eigenvalue equation using Algorithm 1

$\tensor{\mathcal{G}}{}_{\textbf{l,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}=\tensor{\mathcal{V}}{}_{\textbf{l},q}S^{2}_{q,p},\hskip
5.69046ptq,p=1:r,\hskip 5.69046ptS_{q,p}=diag(\sigma_{1},\dots,\sigma_{r})$

3\. Derive the orthogonal tensors

$\tensor{\mathcal{U}}{}_{\textbf{i},q}=\tensor{\mathcal{A}}{}_{\textbf{i,j}}\tensor{\mathcal{V}}{}_{\textbf{j},p}S^{-1}_{p,q}$

OUTPUT: the decomposition

$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\sigma_{p}\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{V}}{}_{\textbf{j},p},\hskip
5.69046ptp=1:r$

### C. Decomposition of a generic tensor

The previous result can be generalized to an order-$(d+e+f)$ tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,l,k}}\in\mathbb{R}^{I\times J\times K}$. In
this case $\tensor{\mathcal{A}}{}_{\textbf{i,l,k}}$ does not unequivocally
represents a transformation. With reference to this tensor we state the
following proposition.

#### Proposition3

Given

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\in\mathbb{R}^{I\times J\times
K},\hskip 5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{e})},\hskip
5.69046pt\textbf{k}={(k_{1},\ldots,k_{f})}$ (67)

where ${I=I_{1}\times I_{2}\times\ldots\times I_{d}}$, ${J=J_{1}\times
J_{2}\times\ldots\times J_{e}}$ and ${K=K_{1}\times K_{2}\times\ldots\times
K_{f}}$, there exist sets of orthogonal tensors
$\tensor{\mathcal{U}}{}=\\{\tensor{\mathcal{U}}{}_{\textbf{i},m}\in\mathbb{R}^{I},\hskip
5.69046ptm=1:M\\}$,
$\tensor{\mathcal{Z}}{}=\\{\tensor{\mathcal{Z}}{}_{\textbf{j},m}\in\mathbb{R}^{J},\hskip
5.69046ptm=1:M\\}$,
$\tensor{\mathcal{W}}{}=\\{\tensor{\mathcal{W}}{}_{\textbf{k},m}\in\mathbb{R}^{K},\hskip
5.69046ptm=1:M\\}$, and a set of scalars
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{M}>{0}$ with $M=r_{1}r_{2}$,
such that $\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}$ can be decomposed as

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\lambda_{m}\tensor{\mathcal{U}}{}_{\textbf{i},m}\tensor{\mathcal{Z}}{}_{\textbf{j},m}\tensor{\mathcal{W}}{}_{\textbf{k},m},\hskip
5.69046ptm=1:M.$ (68)

#### Proof

The product
$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\tensor{\mathcal{A}}{}_{\textbf{l,j,k}}$
is a self-adjoint operator on $\mathbb{R}^{I}$, thus the eigenvalue equation

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\tensor{\mathcal{A}}{}_{\textbf{l,j,k}}\hat{\tensor{\mathcal{U}}{}}_{\textbf{l},p}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q}S^{2}_{q,p},\hskip
5.69046ptq,p=1:r_{1}$ (69)

admits a set of solutions
$\tensor{\mathcal{U}}{}=\\{\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q},\hskip
5.69046ptq=1:r_{1}\\}$ with $S_{q,p}=diag(\sigma_{1},\dots,\sigma_{r_{1}})$,
$\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{r_{1}}>{0}$. Since the matrix
$S_{q,p}$ is diagonal full-rank, thus $S^{-1}_{q,p}$ exists and multiplying by
$\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q}$ (69) can be rewritten as

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q}S^{-1}_{q,p}(\tensor{\mathcal{A}}{}_{\textbf{l,j,k}}\hat{\tensor{\mathcal{U}}{}}_{\textbf{l},p}S^{-1}_{q,p})=\delta_{p,q},$
(70)

that represents the orthogonal condition for the tensor

$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},q}=\tensor{\mathcal{A}}{}_{\textbf{l,j,k}}\hat{\tensor{\mathcal{U}}{}}_{\textbf{l},p}S^{-1}_{p,q}\in{\mathbb{R}^{J\times
K}}.$ (71)

Due to orthogonality of $\hat{\tensor{\mathcal{U}}{}}_{\textbf{l},p}$ we have

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q}S_{q,p}\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}$
(72)

that represents the decomposition of linear transformation
$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}$ from ${\mathbb{R}^{J\times K}}$ to
${\mathbb{R}^{I}}$ . (72) can be further decomposed since the product
$\tensor{\mathcal{V}}{}_{\textbf{j,k},p}\tensor{\mathcal{V}}{}_{\textbf{l,k},p}$
is a self-adjoint operator on $\mathbb{R}^{J}$, thus the eigenvalue equation

$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}\hat{\tensor{\mathcal{V}}{}}_{\textbf{l,k},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{l},r}=\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\Lambda^{2}_{s,r},\hskip
5.69046ptr,s=1:r_{2}$ (73)

admits a set of solutions
$\hat{\tensor{\mathcal{Z}}{}}=\\{\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s},\hskip
5.69046pts=1:r_{2}\\}$ with
$\Lambda_{r,s}=diag(\gamma_{1},\dots,\gamma_{r_{2}})$,
$\gamma_{1}\geq\gamma_{2}\geq\dots\geq\gamma_{r_{2}}>{0}$. The matrix
$\Lambda^{-1}_{s,r}$ exists and (73) can be rewritten as

$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\Lambda^{-1}_{s,r}(\hat{\tensor{\mathcal{V}}{}}_{\textbf{l,k},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{l},r}\Lambda^{-1}_{r,s})=\delta_{r,s}$
(74)

that represents the orthogonality condition for the tensor

$\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},p,r}=\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\Lambda^{-1}_{s,r}.$
(75)

By deriving $\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}$ we have

$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}=\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\Lambda_{s,r}\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},p,r}$
(76)

that represents the decomposition of linear transformation
$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}$ from ${\mathbb{R}^{K}}$ to
${\mathbb{R}^{J}}$. Finally combining (72) and (76) yields

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=S_{q,p}\Lambda_{s,r}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},p,r},$
(77)

and more concisely

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\sigma_{p}\gamma_{p}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},r}\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},p,r},$
(78)

since $S_{q,p}=\sigma_{p}\delta_{q,p}$ and
$\Lambda_{s,r}=\gamma_{r}\delta_{s,r}$. Now, let us derive the matrix $T$ that
establishes the transformation between the vector index
$\textbf{u}=(p,r),p=1:r_{1},r=1:r_{2}$ and the linear index $m=1:M$ with
$M=r_{1}{r_{2}}$, such that

$(p,r)=\textbf{u}(m)=(u_{1}(m),u_{2}(m)).$ (79)

Then (78) becomes

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\sigma_{u_{1}(m)}\gamma_{u_{2}(m)}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},u_{1}(m)}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},u_{2}(m)}\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},u_{1}(m),u_{2}(m)},$
(80)

By defining

$\lambda_{m}=\sigma_{u_{1}(m)}\gamma_{u_{2}(m)}$
$\tensor{\mathcal{U}}{}_{\textbf{i},m}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},u_{1}(m)}$
$\tensor{\mathcal{Z}}{}_{\textbf{j},m}=\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},u_{2}(m)}$
$\tensor{\mathcal{W}}{}_{\textbf{k},m}=\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},u_{1}(m),u_{2}(m)},$
(81)

from (79) we have

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\lambda_{m}\tensor{\mathcal{U}}{}_{\textbf{i},m}\tensor{\mathcal{Z}}{}_{\textbf{j},m}\tensor{\mathcal{W}}{}_{\textbf{k},m},$
(82)

and this concludes the proof.

A pseudo-code of the procedure previously described to derive a decomposition
of a generic tensor $\tensor{\mathcal{A}}{}_{\textbf{i,l,k}}$ is described in
Algorithm 3.

Algorithm 3 Decomposition of a generic tensor

INPUT: order-$(d+e+f)$ tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\in\mathbb{R}^{I\times J\times K}$

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}},\hskip
5.69046pt\textbf{i}={(i_{1},\ldots,i_{d})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{e})},\hskip
5.69046pt\textbf{k}={(k_{1},\ldots,k_{f})}$

1\. Solve eigenvalue equation

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\tensor{\mathcal{A}}{}_{\textbf{l,j,k}}\hat{\tensor{\mathcal{U}}{}}_{\textbf{l},p}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},q}S^{2}_{q,p},\hskip
5.69046ptq,p=1:r_{1}$

2\. Derive the orthogonal tensors

$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},q}=\tensor{\mathcal{A}}{}_{\textbf{l,j,k}}\hat{\tensor{\mathcal{U}}{}}_{\textbf{l},p}S^{-1}_{p,q}\in{\mathbb{R}^{J\times
K}}$

3\. Solve the eigenvalue equation

$\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}\hat{\tensor{\mathcal{V}}{}}_{\textbf{l,k},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{l},r}=\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\Lambda^{2}_{s,r},\hskip
5.69046ptr,s=1:r_{2}$

4\. Derive the orthogonal tensors

$\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},p,r}=\hat{\tensor{\mathcal{V}}{}}_{\textbf{j,k},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},s}\Lambda^{-1}_{s,r}$

5\. Decompose the tensor $\tensor{\mathcal{A}}{}$

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\sigma_{p}\gamma_{p}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},p}\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},r}\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},p,r}$

6\. Apply the inverse vector-to-linear transformation to the vector index u

$(p,r)=\textbf{u}(m)=(u_{1}(m),u_{2}(m))$

7\. Define

$\lambda_{m}=\sigma_{u_{1}(m)}\gamma_{u_{2}(m)}$

$\tensor{\mathcal{U}}{}_{\textbf{i},m}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},u_{1}(m)}$

$\tensor{\mathcal{Z}}{}_{\textbf{j},m}=\hat{\tensor{\mathcal{Z}}{}}_{\textbf{j},u_{2}(m)}$

$\tensor{\mathcal{W}}{}_{\textbf{k},m}=\hat{\tensor{\mathcal{W}}{}}_{\textbf{k},u_{1}(m),u_{2}(m)}$

OUTPUT: the decomposition

$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}=\lambda_{m}\tensor{\mathcal{U}}{}_{\textbf{i},m}\tensor{\mathcal{Z}}{}_{\textbf{j},m}\tensor{\mathcal{W}}{}_{\textbf{k},m}$

##  4 Numerical experiments

The aim of this section is to present some numerical experiments to validate
the theoretical results previously derived for decomposition of tensor
transformations. The first experiment refers to an order-2d tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times I}$ with $d=3$
and $I=[16,16,3]$, to validate the decomposition of an SA-NND tensor operator,
proven in Proposition1. The entries of tensor $\tensor{\mathcal{A}}{}$ have
been randomly generated and the behavior of eigenvalues as obtained by
Algorithm 1 are shown in Fig.1. As you can see all the eigenvalues are
positive and they decreases rapidly as the index $n$ increases. Tensor
$\tensor{\mathcal{A}}{}$ can be reconstructed by decomposition (36) to less
than a numerical error of $2.6485\times 10^{-10}$. The second experiment deals
with the decomposition of an order-$(1+2)$ tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}\in\mathbb{R}^{I\times J}$, with
$I=[64]$ and $J=[8,4]$. In this case the decomposition is obtained by applying
Algorithm 2. Fig.2 shows the behavior of eigenvalues and using the
decomposition (58) the tensor $\tensor{\mathcal{A}}{}$ can be reconstructed to
less than a numerical error of $2.4618\times 10^{-10}$. Finally the third
experiment refers to a tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j,k}}\in\mathbb{R}^{I\times J\times K}$
with $I=[64]$, $J=[16]$, $K=[3]$ and the decomposition described in Algorithm
3. The behavior of eigenvalues is shown in Fig.3 and the reconstruction of
tensor $\tensor{\mathcal{A}}{}$ is obtained with an error of $5.2411\times
10^{-14}$.

Figure 1: The behavior of eigenvalues as obtained from eq.(48) Figure 2: The
behavior of eigenvalues as obtained from eq.(58)  Figure 3: The behavior of
eigenvalues as obtained from eq.(69)

##  5 Conclusion

In this paper a mathematical framework for exact tensor decomposition, that is
able to represent a tensor as the sum of a finite number of low-rank tensors,
has been developed. The core of the proposed approach is the derivation of a
decomposition for non-negative self-adjoint tensor operators. In particular it
has been proven that a correspondence exists between the spectral
decomposition of a symmetric matrix and the decomposition of a self-adjoint
tensor operator.

## References

  * [1] A. K. Smilde, P. Geladi, and R. Bro, _Multi-way analysis: applications in the chemical sciences_. John Wiley & Sons, 2005.
  * [2] L. De Lathauwer and B. De Moor, “From matrix to tensor: Multilinear algebra and signal processing,” in _Institute of mathematics and its applications conference series_ , vol. 67. Citeseer, 1998, pp. 1–16.
  * [3] E. Acar and B. Yener, “Unsupervised multiway data analysis: A literature survey,” _IEEE transactions on knowledge and data engineering_ , vol. 21, no. 1, pp. 6–20, 2008.
  * [4] L. R. Tucker, “Some mathematical notes on three-mode factor analysis,” _Psychometrika_ , vol. 31, no. 3, pp. 279–311, 1966.
  * [5] J. D. Carroll and J.-J. Chang, “Analysis of individual differences in multidimensional scaling via an n-way generalization of “eckart-young” decomposition,” _Psychometrika_ , vol. 35, no. 3, pp. 283–319, 1970.
  * [6] C. J. Appellof and E. R. Davidson, “Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents,” _Analytical Chemistry_ , vol. 53, no. 13, pp. 2053–2056, 1981.
  * [7] S. E. Sofuoglu and S. Aviyente, “Multi-branch tensor network structure for tensor-train discriminant analysis,” _IEEE Transactions on Image Processing_ , vol. 30, pp. 8926–8938, 2021.
  * [8] W. Qin, H. Wang, F. Zhang, J. Wang, X. Luo, and T. Huang, “Low-rank high-order tensor completion with applications in visual data,” _IEEE Transactions on Image Processing_ , vol. 31, pp. 2433–2448, 2022.
  * [9] J.-L. Wang, T.-Z. Huang, X.-L. Zhao, T.-X. Jiang, and M. K. Ng, “Multi-dimensional visual data completion via low-rank tensor representation under coupled transform,” _IEEE Transactions on Image Processing_ , vol. 30, pp. 3581–3596, 2021.
  * [10] L. Sun, C. He, Y. Zheng, Z. Wu, and B. Jeon, “Tensor cascaded-rank minimization in subspace: a unified regime for hyperspectral image low-level vision,” _IEEE Transactions on Image Processing_ , vol. 32, pp. 100–115, 2022.
  * [11] Z. Long, C. Zhu, J. Liu, and Y. Liu, “Bayesian low rank tensor ring for image recovery,” _IEEE Transactions on Image Processing_ , vol. 30, pp. 3568–3580, 2021.
  * [12] M. Zhang, Y. Gao, C. Sun, and M. Blumenstein, “Robust tensor decomposition for image representation based on generalized correntropy,” _IEEE Transactions on Image Processing_ , vol. 30, pp. 150–162, 2020.
  * [13] X. Tian, K. Xie, and H. Zhang, “A low-rank tensor decomposition model with factors prior and total variation for impulsive noise removal,” _IEEE Transactions on Image Processing_ , vol. 31, pp. 4776–4789, 2022.
  * [14] D. Kaur, G. S. Aujla, N. Kumar, A. Y. Zomaya, C. Perera, and R. Ranjan, “Tensor-based big data management scheme for dimensionality reduction problem in smart grid systems: Sdn perspective,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 30, no. 10, pp. 1985–1998, 2018.
  * [15] Y. Ji, Q. Wang, X. Li, and J. Liu, “A survey on tensor techniques and applications in machine learning,” _IEEE Access_ , vol. 7, pp. 162 950–162 990, 2019.
  * [16] M. Signoretto, Q. Tran Dinh, L. De Lathauwer, and J. A. Suykens, “Learning with tensors: a framework based on convex optimization and spectral regularization,” _Machine Learning_ , vol. 94, pp. 303–351, 2014.
  * [17] M. Hou and B. Chaib-draa, “Fast recursive low-rank tensor learning for regression.” in _IJCAI_ , 2017, pp. 1851–1857.
  * [18] V. Lebedev, Y. Ganin, M. Rakhuba, I. Oseledets, and V. Lempitsky, “Speeding-up convolutional neural networks using fine-tuned cp-decomposition,” _arXiv preprint arXiv:1412.6553_ , 2014.
  * [19] Y.-D. Kim, E. Park, S. Yoo, T. Choi, L. Yang, and D. Shin, “Compression of deep convolutional neural networks for fast and low power mobile applications,” _arXiv preprint arXiv:1511.06530_ , 2015.
  * [20] S. Lin, R. Ji, C. Chen, D. Tao, and J. Luo, “Holistic cnn compression via low-rank decomposition with knowledge transfer,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 41, no. 12, pp. 2889–2905, 2018\.
  * [21] N. D. Sidiropoulos, R. Bro, and G. B. Giannakis, “Parallel factor analysis in sensor array processing,” _IEEE transactions on Signal Processing_ , vol. 48, no. 8, pp. 2377–2388, 2000.
  * [22] P. Sankaranarayanan, T. E. Schomay, K. A. Aiello, and O. Alter, “Tensor gsvd of patient-and platform-matched tumor and normal dna copy-number profiles uncovers chromosome arm-wide patterns of tumor-exclusive platform-consistent alterations encoding for cell transformation and predicting ovarian cancer survival,” _PloS one_ , vol. 10, no. 4, p. e0121396, 2015.
  * [23] L. Omberg, G. H. Golub, and O. Alter, “A tensor higher-order singular value decomposition for integrative analysis of dna microarray data from different studies,” _Proceedings of the National Academy of Sciences_ , vol. 104, no. 47, pp. 18 371–18 376, 2007.
  * [24] H. Zhou, L. Li, and H. Zhu, “Tensor regression with applications in neuroimaging data analysis,” _Journal of the American Statistical Association_ , vol. 108, no. 502, pp. 540–552, 2013.
  * [25] J. Liu, P. Musialski, P. Wonka, and J. Ye, “Tensor completion for estimating missing values in visual data,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 35, no. 1, pp. 208–220, 2012.
  * [26] Y. Han, R. Chen, and C.-H. Zhang, “Rank determination in tensor factor model,” _Electronic Journal of Statistics_ , vol. 16, no. 1, pp. 1726–1803, 2022.
  * [27] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM review_ , vol. 51, no. 3, pp. 455–500, 2009.
  * [28] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa, and H. A. Phan, “Tensor decompositions for signal processing applications: From two-way to multiway component analysis,” _IEEE signal processing magazine_ , vol. 32, no. 2, pp. 145–163, 2015.
  * [29] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, “Tensor decomposition for signal processing and machine learning,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 13, pp. 3551–3582, 2017.
  * [30] H. Chen, S. A. Vorobyov, H. C. So, F. Ahmad, and F. Porikli, “Introduction to the special issue on tensor decomposition for signal processing and machine learning,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, no. 3, pp. 433–437, 2021.
  * [31] L. Tucker, “Implications of factor analysis of three way matrices formeasurement of change,” _Problems inmeasuring change. Madison, Wis.: University of WisconsinPress_ , 1963.
  * [32] R. A. Harshman _et al._ , “Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis,” 1970.
  * [33] K. Inoue, “Generalized tensor pca and its applications to image analysis,” _Applied Matrix and Tensor Variate Data Analysis_ , pp. 51–71, 2016.
  * [34] J. Liu, S. Chen, Z.-H. Zhou, and X. Tan, “Generalized low-rank approximations of matrices revisited,” _IEEE Transactions on Neural Networks_ , vol. 21, no. 4, pp. 621–632, 2010.
  * [35] H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos, “Mpca: Multilinear principal component analysis of tensor objects,” _IEEE transactions on Neural Networks_ , vol. 19, no. 1, pp. 18–39, 2008.
  * [36] J. Mažgut, P. Tiňo, M. Bodén, and H. Yan, “Dimensionality reduction and topographic mapping of binary tensors,” _Pattern Analysis and Applications_ , vol. 17, pp. 497–515, 2014.
  * [37] Y. Panagakis, C. Kotropoulos, and G. R. Arce, “Non-negative multilinear principal component analysis of auditory temporal modulations for music genre classification,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 18, no. 3, pp. 576–588, 2009.
  * [38] J. Sun, D. Tao, S. Papadimitriou, P. S. Yu, and C. Faloutsos, “Incremental tensor analysis: Theory and applications,” _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , vol. 2, no. 3, pp. 1–37, 2008.
  * [39] D. Xu, S. Yan, L. Zhang, S. Lin, H.-J. Zhang, and T. S. Huang, “Reconstruction and recognition of tensor-based objects with concurrent subspaces analysis,” _IEEE Transactions on Circuits and Systems for Video Technology_ , vol. 18, no. 1, pp. 36–47, 2008.
  * [40] J. Yang, D. Zhang, A. F. Frangi, and J.-y. Yang, “Two-dimensional pca: a new approach to appearance-based face representation and recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 26, no. 1, pp. 131–137, 2004.
  * [41] J. Ye, “Generalized low rank approximations of matrices,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 112.
  * [42] J. Ye, R. Janardan, and Q. Li, “Gpca: An efficient dimension reduction scheme for image compression and retrieval,” in _Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining_ , 2004, pp. 354–363.
  * [43] R. Bro, “Parafac. tutorial and applications,” _Chemometrics and intelligent laboratory systems_ , vol. 38, no. 2, pp. 149–171, 1997.
  * [44] N. K. M. Faber, R. Bro, and P. K. Hopke, “Recent developments in candecomp/parafac algorithms: a critical review,” _Chemometrics and Intelligent Laboratory Systems_ , vol. 65, no. 1, pp. 119–137, 2003.
  * [45] J. B. Kruskal, “Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics,” _Linear algebra and its applications_ , vol. 18, no. 2, pp. 95–138, 1977\.
  * [46] A. Shashua and A. Levin, “Linear image coding for regression and classification using the tensor-rank principle,” in _Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001_ , vol. 1. IEEE, 2001, pp. I–I.
  * [47] L. Grasedyck, “Hierarchical singular value decomposition of tensors,” _SIAM journal on matrix analysis and applications_ , vol. 31, no. 4, pp. 2029–2054, 2010.
  * [48] L. Grasedyck, D. Kressner, and C. Tobler, “A literature survey of low-rank tensor approximation techniques,” _GAMM-Mitteilungen_ , vol. 36, no. 1, pp. 53–78, 2013.
  * [49] W. Hackbusch and S. Kühn, “A new scheme for the tensor representation,” _Journal of Fourier analysis and applications_ , vol. 15, no. 5, pp. 706–722, 2009.
  * [50] J. A. Bengua, P. N. Ho, H. D. Tuan, and M. N. Do, “Matrix product state for higher-order tensor compression and classification,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 15, pp. 4019–4030, 2017.
  * [51] M. Chaghazardi and S. Aeron, “Sample, computation vs storage tradeoffs for classification using tensor subspace models,” _arXiv preprint arXiv:1706.05599_ , 2017.
  * [52] J. Håstad, “Tensor rank is np-complete,” in _Automata, Languages and Programming: 16th International Colloquium Stresa, Italy, July 11–15, 1989 Proceedings 16_. Springer, 1989, pp. 451–460.
  * [53] A. Zare, A. Ozdemir, M. A. Iwen, and S. Aviyente, “Extension of pca to higher order data structures: An introduction to tensors, tensor decompositions, and tensor pca,” _Proceedings of the IEEE_ , vol. 106, no. 8, pp. 1341–1358, 2018.
  * [54] L. Cheng, Y.-C. Wu, and H. V. Poor, “Probabilistic tensor canonical polyadic decomposition with orthogonal factors,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 3, pp. 663–676, 2016.