# Tensor PCA from basis in tensor space

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

###### Abstract

The aim of this paper is to present a mathematical framework for tensor PCA.
The proposed approach is able to overcome the limitations of previous methods
that extract a low dimensional subspace by iteratively solving an optimization
problem. The core of the proposed approach is the derivation of a basis in
tensor space from a real self-adjoint tensor operator, thus reducing the
problem of deriving a basis to an eigenvalue problem. Three different cases
have been studied to derive: i) a basis from a self-adjoint tensor operator;
ii) a rank-1 basis; iii) a basis in a subspace. In particular, the equivalence
between eigenvalue equation for a real self-adjoint tensor operator and
standard matrix eigenvalue equation has been proven. For all the three cases
considered, a subspace approach has been adopted to derive a tensor PCA.
Experiments on image datasets validate the proposed mathematical framework.

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

##  1 Introduction

Tensors are arrays indexed by multiple indices, a generalization of vectors
and matrices [1], that have been adopted in several branches of data analysis,
thanks to their ability to represent a wide range of real-world data. Tensors
are most popular in machine learning [2] [3], [4], [5], [6], [7], [8], [9],
and image processing [10], [11], [12], [13], [14], [15], [16], [17], but they
arise in numerous applications such as: sensor array processing [18], big data
representation [19], completion of big data [20], genomics [21], neuroimaging
data analysis [22], [23], tumor classification [24], computer vision [25],
finance data [26], human action and gesture recognition [27]. Tensorial data
are typically high dimensional involving high computational cost, thus
techniques for dimensionality reduction are required to make possible the
application of tensors to real data. A widely adopted technique for this
purpose is to represent high-dimensional tensors in a low-dimensional
subspace, by preserving useful underlying information or structures. In this
context several very powerful tensor decomposition approaches have been
developed in the past [28], [29], [30], [31]. Among these, Tucker
decomposition [32] and CANDECOMP/PARAFAC (CP) decomposition [33], [34] 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), tensor
rank-one decomposition (TROD), hierarchical Tucker PCA (HT-PCA) and tensor-
train PCA (TT-PCA).
MPCA [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], aims to find
low-rank orthogonal projection matrices forming a lower dimensional subspace,
such that the projection of a set of higher order data on this subspace
minimizes a mean square error. The solution is iterative and performs feature
extraction that captures most of the original tensorial input variations.
Similarly to MPCA, the TROD approach [45], [46], [33], [47], [48] formulates
the decomposition in terms of subspaces whose bases are rank-1 tensors. More
formally, TROD aims to find a set of vectors forming a subspace, such that the
projections of tensor data on this subspace minimizes a mean square error. An
efficient greedy least square minimization procedure is used to iteratively
compute each vector leaving all the others fixed.
Both HT and TT approaches fall within the class of tensor networks (TNs), that
represents a higher order tensor as a set of sparsely interconnected lower
order tensors, providing computational and storage benefits [49], [50], [51].
HT-PCA [52], [53], [54], which has been proposed to reduce the memory
requirements for Tucker decomposition, recursively splits the modes based on a
hierarchy and creates a binary tree $T$ combining a subset of modes at each
node. Thus, factor matrices are obtained from the singular value decomposition
(SVD) of the matricization of a tensor corresponding to each node. In current
applications a suboptimal algorithm is used [55], as estimating both tree and
subspace is an NP hard problem.
TT-PCA decomposition [56], [55], is a special case of HT in which all nodes of
the TN are connected in a cascade or train. Given a set of tensor data the
objective of TT-PCA is to find order-3 tensors such that a left unfolding
operator resulting in a matrix that represents a subspace is derived. In [57]
an approach based on successive SVD algorithm for computing TT decomposition
followed by thresholding the singular values, is proposed.
All the above methods for tensor decomposition try to find a low-dimensional
subspace by solving an optimization problem. A commonly used approach for this
purpose is the alternating least square (ALS) method [58], 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 [59].
In order to overcome these issues, this paper proposes an alternative approach
to tensor PCA that does not require solving an optimization problem. Tensor
PCA can be considered as multilinear extensions of principal component
analysis (PCA) [60], [61], [62], [63], [64], [65], one of the most popular and
efficient subspace learning techniques used in matrix analysis to reduce the
dimensionality of a dataset. More formally, PCA derives an orthogonal linear
transformation to extract a low-dimensional representation from high-
dimensional data. This transformation can be regarded as a projection of data
onto a coordinate system, termed as principal components (PCs), that generates
the maximum variance. There are two used approaches to perform the PCA: i) the
eigenvalue decomposition of data covariance matrix; and ii) the SVD of data
matrix. In both cases an orthogonal matrix $U$ is derived from the so-called
Gram-matrix $G=X^{T}X$, where $X$ is the data matrix [60], while $X^{T}$ its
transpose. The matrix $U$ represents the orthogonal linear transformation,
whose columns are the eigenvectors of Gram-matrix and the existence of matrix
$U$ is ensured by the symmetry of $G$. More specifically, the Gram-matrix can
be interpreted as a real nonnegative self-adjoint operator that admits a set
of real nonnegative eigenvalues, corresponding to a set of eigenvectors. Due
to orthogonality of eigenvectors, the matrix $U$ represents a basis in the
space of data vector $x$, whose observations are the rows of data matrix $X$.
With these considerations in mind, the PCA can be contextualized in the more
general framework of finite-dimensional linear spaces [66] and in this context
PCA can be considered as equivalent to derive a basis in vector space from a
real self-adjoint operator. As this concept is quite general can be profitably
used to extend PCA from vector-space to tensor-space. Thus in this paper some
results of the theory of self-adjoint operators will be used as the foundation
to derive tensor PCA.
The aim of this paper is to develop a mathematical framework for tensor PCA.
The core of the proposed approach is the derivation of a basis in tensor space
from a real self-adjoint operator. It will be shown that tensors with a
symmetrical structure belong to this important class of operators, can be
easily derived as tensor products from tensor data, similar to the way in
which Gram-matrix are derived from data matrix. In this way the problem of
deriving a basis is equivalent to an eigenvalue problem. As final objective,
tensors PCA is obtained as the projection of tensor data on the subspace
spanned by the coordinates corresponding to the highest eigenvalues. In the
paper three different problems will be carried out to derive: i) a basis from
a real self-adjoint operator; ii) a rank-1 basis from a set of tensors; iii) a
basis in a subset of tensors. For all these three cases a subspace approach
will be adopted to derive a tensor PCA. 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. This result represents an
advancement in the context of multilinear analysis, as it constitutes a
generalization of previous results limited only to the case of supersymmetric
tensor operators [67], [68], [69], [70], [71].
The paper is organized as follows. Section 2 reports some preliminaries and
notation. Section 3 deals with the finite linear space of tensors and some
fundamental concepts to derive a tensor basis. In Section 4 the mathematical
framework for tensor PCA is developed as projection of data on a subspace with
reference to three different problems. Section 5 validates the proposed
approach by experiments on image datasets.

##  2 Preliminaries and notation

A real order-$d$ tensor is a multidimensional array denoted by Euler script
calligraphic letters, e.g. $\tensor{\mathcal{X}}{}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{d}}$, where $\times$ represents the Cartesian
product. The number $d$ of dimensions, also known as modes, is the order of a
tensor. The generic element of an order-$d$ tensor $\tensor{\mathcal{X}}{}$
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.$
(1)

or

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

where the index range in the $k$th mode is $1\leq{i_{k}}\leq{I_{k}},\hskip
5.69046ptk=1:d$. The vector i, in bold font, is a subscript vector which
assumes $L=I_{1}I_{2}\ldots{I_{d}}$ distinct values in the Cartesian interval
$I=I_{1}\times I_{2}\times\ldots\times I_{d}$. $L$ represents the dimension of
the $\tensor{\mathcal{X}}{}$ domain that will denoted by $dim(I)$.

#### Tensor product

Given two tensors

$\tensor{\mathcal{X}}{}(\textbf{i})=x{(i_{1},\ldots,i_{d})},\hskip
14.22636pt\tensor{\mathcal{Y}}{}(\textbf{j})=y{(j_{1},\ldots,j_{f})}$ (3)

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})$ (4)

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}}{}(\textbf{l},\textbf{m},\textbf{p},\textbf{q})=\sum_{\textbf{k=1}}^{\textbf{K}}\tensor{\mathcal{X}}{}(\textbf{l},\textbf{k},\textbf{m})\tensor{\mathcal{Y}}{}(\textbf{p},\textbf{k},\textbf{q})$
(5)

Here a multi-index notation is used for summation, meaning that, if
$\textbf{k}={(k_{1},\ldots,k_{g})}$ is a length-g index vector, then

$\sum_{\textbf{k=1}}^{\textbf{K}}=\sum_{{k_{1}}=1}^{{K_{1}}}...\sum_{{k_{g}}=1}^{{K_{g}}}$
(6)

As a result the product (6) corresponds to a contraction of common index k.
Using the Einstein summation convention, that interprets repeated subscript as
summation over that index, (6) can be rewritten in a more compact form as

$\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}}$ (7)

that implies summation over the vector index k. In the following we will
extensively use the Einstein convention to simplify mathematical notation.

#### Inner product

The inner product of two tensors of the same size
$\tensor{\mathcal{X}}{},\,\tensor{\mathcal{Y}}{}\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}}$ (8)

and $\tensor{\mathcal{X}}{},\tensor{\mathcal{Y}}{}$ are said to be orthogonal
if $\left\langle\tensor{\mathcal{X}}{},\tensor{\mathcal{Y}}{}\right\rangle=0$.
From this definition 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}.$
(9)

#### Outer product

The outer product of tensor
$\tensor{\mathcal{X}}{}=\tensor{\mathcal{X}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{q}}$ with tensor
$\tensor{\mathcal{Y}}{}=\tensor{\mathcal{X}}{}_{\textbf{j}}\in\mathbb{R}^{J_{1}\times
J_{2}\times\ldots\times J_{p}}$ is the order-$(q+p)$ 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_{q})},\hskip
5.69046pt\textbf{j}={(j_{1},\ldots,j_{p})}$ (11)

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

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

#### Vector-to-linear index transformation

This operation turns 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)

to an order-$(p+1)$ tensor

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

transforming the vector index i to the linear index $m$, ranging from 1 to
$L$, (1:$L$ in Matlab notation) where $L=I_{1}I_{2}\ldots{I_{q}}$.

This is achieved by noticing that the indices $i_{1},\ldots,i_{q}$ can be
arranged in $L$ different ways (the number of permutations with repetition),
so that a correspondence $\alpha$ between $m$ and i

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

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{L},\hskip
14.22636pt1\leq{i_{k}}\leq{I_{k}},\hskip 5.69046ptk=1:q$ (16)

Among the possible choices for the correspondence, here we use the one
described by the following relationship

$m=\alpha(\textbf{i})=i_{1}+(i_{2}-1)I_{1}+(i_{3}-1)I_{1}I_{2}+\ldots+(i_{q}-1)I_{1}I_{2}\ldots{I_{q-1}}$
(17)

With reference to (17) the matrix $T$ whose $m$th row represents the
corresponding vector index $\textbf{i}(m)$ can be derived, so that the inverse
transformation of (17) for a value $m$ is formally given by

$\textbf{i}=T(m,:).$ (18)

An example of such a matrix for $\textbf{i}=(i_{1},i_{2},i_{3})$,
$I_{1}=3,I_{2}=2,I_{3}=2$ is the following

$T=\left(\begin{array}[]{c c c c c c c c c c c c}1&1&1&1&1&1&2&2&2&2&2&2\\\
1&1&1&2&2&2&1&1&1&2&2&2\\\ 1&2&3&1&2&3&1&2&3&1&2&3\\\ \end{array}\right)^{T}$

where $1\leq{m}\leq{12}$ represents the row-index.

##  3 The space of tensors

It is straightforward to show that the set of order-$d$ tensors
$\tensor{\mathcal{X}}{}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times
I_{d}}$ forms a linear space, since it is closed under addition and
multiplication by a scalar. With reference to the canonical basis
$e_{i_{1}}=(1,0,\ldots,0),\ldots,e_{i_{d}}=(0,0,\ldots,1)$ in vector space,
any order-$d$ tensor $\tensor{\mathcal{X}}{}$ can be decomposed as

$\tensor{\mathcal{X}}{}=\tensor{\mathcal{X}}{}_{i_{1}i_{2}\ldots
i_{d}}\;e_{i_{1}}e_{i_{2}}\ldots e_{i_{d}}$ (19)

The outer product $e_{i_{1}}e_{i_{2}}\ldots e_{i_{d}}$ is equivalent to a set
of $L=\prod_{j=1}^{d}I_{j}$ order-$d$ tensors, thus due to property (19) it
represents a canonical tensor basis. As this basis is of size $L$, the set of
order-$d$ tensors form a linear space of dimensionality $L$. As an example for
$d=2$ and $I_{1}=I_{2}=3$ we have $L=9$ and

$e_{1}e_{1}=\left(\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\
0&0&0\end{array}\right),\;$

$e_{1}e_{2}=\left(\begin{array}[]{ccc}0&1&0\\\ 0&0&0\\\
0&0&0\end{array}\right),\;$ $e_{1}e_{3}=\left(\begin{array}[]{ccc}0&0&1\\\
0&0&0\\\ 0&0&0\end{array}\right),\;\ldots$

#### Tensor operators

In the product (7) the tensor $\tensor{\mathcal{X}}{}_{\hskip
0.56917pt\textbf{l},\textbf{k},\textbf{m}}$ can be interpreted as a
transformation from space $R^{P\times{K}\times{Q}}$ to space
$R^{L\times{M}\times{P}\times{Q}}$, where $P=p_{1}\times\ldots\times{p_{p}}$,
$K=k_{1}\times\ldots\times{k_{k}}$ and so on, transforming tensor
$\tensor{\mathcal{Y}}{}_{\hskip 0.56917pt\textbf{p},\textbf{k},\textbf{q}}$ to
tensor $\tensor{\mathcal{Z}}{}_{\hskip
0.56917pt\textbf{l},\textbf{m},\textbf{p},\textbf{q}}$. Among the possible
transformations, we will interested in transformations from
$\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{d}}$ to itself. In this
particular case we have

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

where $\textbf{i}={(i_{1},\ldots,i_{d})}$ and
$\textbf{j}={(j_{1},\ldots,j_{d})}$, so that the tensors
$\tensor{\mathcal{Z}}{}_{\hskip 0.56917pt\textbf{i}}$ and
$\tensor{\mathcal{Y}}{}_{\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. Thus, denoting with $\tensor{\mathcal{L}}{}_{d}$
the order-$d$ tensor space $\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times
I_{d}}$ the operator $\tensor{\mathcal{A}}{}_{\textbf{i,j}}$ establishes a
transformation from $\tensor{\mathcal{L}}{}_{d}$ to
$\tensor{\mathcal{L}}{}_{d}$ and (20) can be rewritten without ambiguity as

$\tensor{\mathcal{Z}}{}=\tensor{\mathcal{A}}{}\tensor{\mathcal{Y}}{}$ (21)

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}}{}$
(22)

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

#### Real self-adjoint operators

An important class of operators is the class of real self-adjoint operators
$\tensor{\mathcal{A}}{}$, 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$
(23)

holds.

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})}$ (24)

such that

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

is a real 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.$
(26)

It is straightforward to show that the covariance of tensor
$\tensor{\mathcal{X}}{}_{\textbf{i}}$ defined as

$\tensor{\mathcal{R}}{}_{\textbf{i,j}}=E\\{\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{X}}{}_{\textbf{j}}\\},$
(27)

where $E\\{.\\}$ denotes the mean value, is a self-adjoint operator. In the
particular case of order-1 tensor space , achieved with $d=1$ in (24), the
self-adjoint operator (25) reduces to a symmetric matrix that operates on
space $R^{I_{1}}$.
On the basis of this property it follows that eigentensors of a self-adjoint
operator $\tensor{\mathcal{A}}{}$ belonging to distinct eigenvalues are
orthogonal. In fact 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$
(28)

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$
(29)

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}$.

The eigenvalues of a real self-adjoint operator are real, indeed combining
(22) and (23)
$\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}\right\rangle=\lambda{\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{V}}{}\right\rangle}=\lambda$
that yields a real 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}$
(30)

From (30) it follows that operators for which
$\left\langle\tensor{\mathcal{V}}{},\tensor{\mathcal{A}}{}\tensor{\mathcal{V}}{}\right\rangle\geq
0$ holds, ensure the eigenvalues are nonnegative. Operators belonging to a
such class are said to be nonnegative operators. A nonnegative operator can be
easily derived from a tensor $\tensor{\mathcal{X}}{}_{\textbf{i,j}}$ as
follows

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

For such an operator we have

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

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

${G}_{m,n}={X}_{l,m}{X}_{l,n}=({X}^{T}_{m,l}){X}_{l,n}$

that is $G=X^{T}X$. Similarly, operator
$\tensor{\mathcal{A^{\prime}}}{}_{\textbf{i,j}}$ reduces to the matrix
$G^{T}=XX^{T}$. The matrix $G$ represents the well known Gram-matrix, while
matrix $G^{T}$ its transpose. Thus, tensor
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}$ is the analog in tensor space of Gram-
matrix $G$.
Tensor eigenvalues and eigenvectors have received much attention recently in
the literature [72], [73], [74], following the seminal works of Qi [67], and
Lim [75], that introduced spectral theory of supersymmetric tensors. A tensor
$\tensor{\mathcal{A}}{}$ is called supersymmetric as its entries
$a_{i_{1},\dots i_{m}}$ are invariant under any permutations of their indices.
It is worth to note that tensor (25) has a symmetrical structure but it is not
supersymmeric, thus in order to distinguish it from the class of
supersymmetric tensors, simply referred as "symmetric" in some papers [68],
[69], [70], [71], in the following we will adopt the term self-adjoint for
tensor (25).

##  4 Tensor PCA from basis in tensor space

As discussed in the Introduction the core aspect of principal component
analysis (PCA) is the derivation of a basis from a self-adjoint operator, from
which a representation on a subspace that minimizes an error condition can be
achieved. The aim of this Section is to derive a basis in the order-$d$ tensor
space $\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{d}}$, that is a set
of $L=I_{1}I_{2}\ldots{I_{d}}$ orthogonal tensors

$\tensor{\mathcal{U}}{}_{\textbf{i},l},\hskip 5.69046ptl=1:L$ (33)

such that, for any tensor $\tensor{\mathcal{X}}{}_{\textbf{i}}$ in this space,
the following orthogonal expansion

$\tensor{\mathcal{X}}{}_{\textbf{i}}=\tensor{\mathcal{U}}{}_{\textbf{i},l}d_{l}$
(34)

holds. Using the orthogonality property
${\tensor{\mathcal{U}}{}_{\textbf{i},l}}{\tensor{\mathcal{U}}{}_{\textbf{i},l^{\prime}}}=\delta_{l,l^{\prime}}$
of basis ($\delta_{l,l^{\prime}}$ is the Dirac operator:
$\delta_{l,l^{\prime}}=1$ for $l=l^{\prime}$, while $\delta_{l,l^{\prime}}=0$
for $l\neq{l^{\prime}}$ ), the vector of coefficients $d_{l}$ in (34) can be
easily derived by the inner product of $\tensor{\mathcal{X}}{}_{\textbf{i}}$
and ${\tensor{\mathcal{U}}{}_{\textbf{i},l}}$

$d_{l}=\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{U}}{}_{\textbf{i},l}$
(35)

As final step tensor PCA is derived by extracting low-dimensional subspace
representation of dataset that ensures a minimum reconstruction error. In the
following results will be presented with reference to three different problems
to derive: i) a basis assuming a real self-adjoint tensor operator is given;
ii) a rank-1 basis from a set of tensors; iii) a basis in a subset of tensors.

### A. Basis from real self-adjoint operator

Let us given the order-$2d$ tensor $\tensor{\mathcal{A}}{}$

$\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})}$ (36)

with $1\leq{i_{k}}\leq{I_{k}},\hskip 5.69046pt1\leq{j_{k}}\leq{J_{k}},\hskip
5.69046ptk=1:d$, so that it results

$\tensor{\mathcal{A}}{}\in\mathbb{R}^{I\times J},\hskip
5.69046pt{I=I_{1}\times I_{2}\times\ldots\times I_{d}},\hskip
5.69046pt{J=J_{1}\times J_{2}\times\ldots\times J_{d}}$ (37)

and $dim(I)=dim(J)=I_{1}I_{2}\ldots{I_{d}}=L$. Assuming
$\tensor{\mathcal{A}}{}$ is self-adjoint, i.e.
$\tensor{\mathcal{A}}{}_{\textbf{i,j}}=\tensor{\mathcal{A}}{}_{\textbf{j,i}}$,
then we are interested in solving the following eigenvalue problem, that is to
determine the tensors $\tensor{\mathcal{U}}{}$ (the eigentensors) and the
scalars $\lambda$ (the eigenvalues) such that equation

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

is satisfied, where $\tensor{\mathcal{U}}{}\in\mathbb{R}^{I}$. To this end we
will prove the following proposition.

#### Proposition1

Eq. (38) is equivalent to the eigenvalue equation

$a(n,m)u(m)=\lambda{u(n)},\hskip 5.69046ptm,n=1:L$ (39)

where $a(m,n)$ is a symmetric matrix achieved using the vector-to-linear index
transformation (17) and the eigentensors $\tensor{\mathcal{U}}{}$ are obtained
by applying the inverse vector-to-linear transformation (18) to the
eigenvectors $u$ in (39).

#### Proof

Using index convention introduced in Section 2, Eq. (38) 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)$ (40)

which can be written as

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

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 (41) are independent on how the
terms in the equation are arranged as entries in the tensor
$\tensor{\mathcal{Y}}{}_{\textbf{i}}$, meaning that (41) 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,:)$, (41) 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$
(42)

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)}$ (43)

and the vector

$u(m)=\tensor{\mathcal{U}}{}_{\textbf{j}(m)},$ (44)

thus (42) can be re-arranged as

$(a(n,m))-\lambda\delta(n,m))u(m)=0$ (45)

or equivalently

$a(n,m)u(m)=\lambda{u(n)},\hskip 5.69046ptm,n=1:L$ (46)

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)$
(47)

that proves matrix $a(n,m)$ is symmetric (or self-adjoint), and this concludes
the proof.
Proposition1 shows that by a one-to-one vector-to-linear transformation the
eigenvalue tensor equation (38) reduces to an eigenvalue matrix equation. (39)
is a standard eigenvalue equation for matrix $a(n,m)$, that admits $L$
orthogonal eigenvectors due to the symmetry of $a(n,m)$. Once the eigenvectors
$u(n)$ are obtained by solving (39), they can be converted to tensors
$\tensor{\mathcal{U}}{}$ of proper size using the inverse vector-to-linear
index transformation (17)

$u(n(\textbf{i))}=\tensor{\mathcal{U(\textbf{i})}}{}=\tensor{\mathcal{U}}{}_{\textbf{i}}$
(48)

.

The statement proven with Proposition 1 shows that an eigenvalue equation for
real self-adjoint tensor operators is invariant under one-to-one
transformations of subscript vectors, thus reducing to a matrix eigenvalue
equation. This result represents an advancement in the general higher order
tensor eigenvalue problem, so far limited to the class of supersymmetric
tensor operators [67], [75][72],[73], [74].

Having derived a basis in the order-$d$ tensor space, a generic tensor
$\tensor{\mathcal{X}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{d}}$ in this space can be represented by the
orthogonal expansion (34). A pseudo-code of the procedure previously described
to derive a basis from a self-adjoint operator is described in Algorithm 1.

####  Tensor PCA

On the basis of previous results a tensor PCA can be carried out by deriving a
low-dimensional subspace representation of dataset that minimizes the
reconstruction error. The following proposition establishes this result.

#### Proposition 2

\- Let us refer to a generic basis
$\\{{\tensor{\mathcal{U}}{}_{\textbf{i},l},\hskip 5.69046ptl=1:L}\\}$ such
that for any $\tensor{\mathcal{X}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{d}}$ the representation

$\tensor{\mathcal{X}}{}_{\textbf{i}}=\tensor{\mathcal{U}}{}_{\textbf{i},l}d_{l}$
$d_{l}=\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{U}}{}_{\textbf{i},l}$
(49)

holds, and decompose $\tensor{\mathcal{X}}{}_{\textbf{i}}$ as

$\tensor{\mathcal{X}}{}_{\textbf{i}}={\hat{\tensor{\mathcal{X}}{}}_{\textbf{i}}+\tensor{\mathcal{N}}{}_{\mathbf{i}}}$
(50)

where
$\hat{\tensor{\mathcal{X}}{}}_{\textbf{i}}=\tensor{\mathcal{U}}{}_{\textbf{i},l}d_{l},\hskip
5.69046ptl=1:M$ with $M<L$, represents a subspace representation and
$\tensor{\mathcal{N}}{}_{\mathbf{i}}=\tensor{\mathcal{U}}{}_{\textbf{i},p}d_{p},\hskip
5.69046ptp=M+1:L$ the residual. Then the basis
$\\{{\tensor{\mathcal{U}}{}_{\textbf{i},l},\hskip 5.69046ptl=1:L}\\}$ that
minimizes the error

$\tensor{\mathcal{E}}{}=E{\\{\|\tensor{\mathcal{X}}{}_{\textbf{i}}-{\hat{\tensor{\mathcal{X}}{}}_{i}}\|^{2}}\\}$
(51)

is solution of the eigenvalue equation

$\tensor{\mathcal{R}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p}=\lambda_{p}{\tensor{\mathcal{U}}{}_{\textbf{i},p}}$
(52)

where
$\tensor{\mathcal{R}}{}_{\textbf{i,j}}=E\\{\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{X}}{}_{\textbf{j}}\\}$
is the covariance of $\tensor{\mathcal{X}}{}_{\textbf{i}}$.

A proof of this proposition is given as follows.

#### Proof

From
$\tensor{\mathcal{X}}{}_{\textbf{i}}={\hat{\tensor{\mathcal{X}}{}}_{\textbf{i}}+\tensor{\mathcal{N}}{}_{\mathbf{i}}}$
the mean square error $\tensor{\mathcal{E}}{}$ results as given by

$\tensor{\mathcal{E}}{}=E{\\{\tensor{\mathcal{N}}{}_{\textbf{i}}{{\tensor{\mathcal{N}}{}}_{\textbf{i}}}}\\}=E{\\{\tensor{\mathcal{U}}{}_{\textbf{i},p}d_{p}{{\tensor{\mathcal{U}}{}}_{\textbf{i},q}}d_{q}\\}}=$
$E{\\{\delta_{p,q}d_{p}d_{q}\\}}=E{\\{\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{U}}{}_{\textbf{i},p}{{\tensor{\mathcal{X}}{}_{\textbf{j}}\tensor{\mathcal{U}}{}}_{\textbf{j},q}}\\}}=$
$\tensor{\mathcal{U}}{}_{\textbf{i},p}E\\{\tensor{\mathcal{X}}{}_{\textbf{i}}\tensor{\mathcal{X}}{}_{\textbf{j}}\\}\tensor{\mathcal{U}}{}_{\textbf{j},q}=\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{R}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p},$
(53)

with $p=M+1:L$. Thus to find the basis that minimizes the error
$\tensor{\mathcal{E}}{}$ subject to the orthonormal condition for the basis
$\tensor{\mathcal{U}}{}_{\textbf{i},p}$, is equivalent to the problem

$\min\\{\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{R}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p}\\},\hskip
5.69046pts.t.\hskip
5.69046pt\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{U}}{}_{\textbf{i},p}=1$
(54)

Using the Lagrange multiplier method the problem reduces to the equation

$\frac{\partial}{\partial{\tensor{\mathcal{U}}{}_{\textbf{i},p}}}\\{\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{R}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p}+\lambda_{p}{(1-\tensor{\mathcal{U}}{}_{\textbf{i},p}\tensor{\mathcal{U}}{}_{\textbf{i},p})}\\}=0$
(55)

which once solved yields

$\tensor{\mathcal{R}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p}=\lambda_{p}\tensor{\mathcal{U}}{}_{\textbf{i},p}$
(56)

that is the eigenvalue equation of covariance
$\tensor{\mathcal{R}}{}_{\textbf{i,j}}$, and this concludes the proof.
As a corollary of this proposition, it follows that the basis elements in the
truncated representation
$\hat{\tensor{\mathcal{X}}{}}_{\textbf{i}}=\tensor{\mathcal{U}}{}_{\textbf{i},l}d_{l}$,
$l=1:M$ correspond to the highest eigenvalues in (52). In fact the error

$\min{\tensor{\mathcal{E}}{}}=\tensor{\mathcal{U}}{}_{\textbf{i},p}(\tensor{\mathcal{R}}{}_{\textbf{i,j}}\tensor{\mathcal{U}}{}_{\textbf{j},p})=\tensor{\mathcal{U}}{}_{\textbf{i},p}\lambda_{p}\tensor{\mathcal{U}}{}_{\textbf{i},p}=$
(57)
$\lambda_{p}\|\tensor{\mathcal{U}}{}_{\textbf{i},p}\|^{2}=\sum_{p=M+1}^{L}\lambda_{p}$
(58)

is minimum when the eigenvalues corresponding to the basis elements in
residual are the lowest or, that is the same, when the eigenvalues
corresponding to the basis elements in the subspace are the highest.
Proposition 2 proves that low-dimensional subspace representation
$\hat{\tensor{\mathcal{X}}{}_{\textbf{i}}}$ is optimal, since minimizes the
mean square reconstruction error.

Algorithm 1 Basis from real self-adjoint operator

INPUT: order-$2d$ self-adoint operator

$\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 matrix

$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)=\lambda{u(n)},\hskip 5.69046ptm,n=1:L$

4\. Convert eigenvector $u(n)$ to tensor by inverse vector-to-linear
transformation

$\tensor{\mathcal{U}}{}_{\textbf{i}}=u(n(\textbf{i}))$

OUTPUT: the basis

$\tensor{\mathcal{U}}{}=\\{\tensor{\mathcal{U}}{}_{\mathbf{i},l},\hskip
5.69046ptl=1:L\\}$

### B. Rank-1 basis

In this case, a set

$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},\hskip
5.69046ptn=1:N\\}$ (59)

of $N$ order-$d$ tensors
$\tensor{\mathcal{X}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{d}},$ $\hskip
5.69046pt\textbf{i}=(i_{1},\ldots,i_{d})$, is given. The problem to be solved
is to derive a set of $L=I_{1}I_{2}\ldots{I_{d}}$, rank-$1$ orthogonal tensors
$\tensor{\mathcal{U}}{}_{\textbf{i},\textbf{j}}$ with
${\textbf{j}=(j_{1},\ldots,j_{d})}$ having the same dimension of i, which can
be written as

$\tensor{\mathcal{U}}{}_{\textbf{i,j}}={U}_{i_{1},j_{1}}{U}_{i_{2},j_{2}}\ldots{U}_{i_{d},j_{d}}$
(60)

where the generic element ${U}_{i_{k},j_{k}}\in\mathbb{R}^{I_{k}\times J_{k}}$
is an orthonormal matrix $U^{k}$ that corresponds to the $k$-mode. In such a
way for a given value of j the resulting tensor

${U}_{\textbf{:j}}={U}_{:,j_{1}}{U}_{:,j_{2}}\ldots{U}_{:,j_{d}}$ (61)

is a rank-$1$ tensor being expressed as outer product of $d$ vectors, i.e. the
columns of matrices ${U}_{i_{k}j_{k}}$. The problem can be addressed by
deriving an eigenvector equation for each mode of the tensor
$\tensor{\mathcal{X_{\textbf{i}}}}{}$, so that once solved an orthogonal
matrix is obtained. This is equivalent to find a set of self-adjoint tensors
that operates along the $d$ modes. To this end let us refer to the $k$th mode.
The tensor

$\tensor{\mathcal{X}}{}(i_{1},\ldots,{i_{k}},\ldots,{i_{d}},n)\tensor{\mathcal{X}}{}(i_{1},\ldots,{j_{k}},\ldots,{i_{d}},n)=A^{k}(i_{k},j_{k})$
(62)

achieved by the contraction of all indices but $i_{k}$ and $j_{k}$, defines a
self-adjoint operator along $k$th mode as a symmetric matrix $A^{k}$ of
dimension $(I_{k}\times{J_{k}})$. From this matrix we can derive an orthogonal
matrix $U^{k}$ which is solution of the eigenvector equation

$A^{k}U^{k}=U^{k}\Sigma^{k}$ (63)

whose generic element is ${U}_{i_{k}j_{k}}$. By repeating this approach for
all the modes of tensor, we finally obtain the required $d$ orthogonal
matrices that can be used in (60) to derive a set of L orthogonal rank-$1$
tensors

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

Transforming vector index i to linear index by (18) yields the set

$\hat{\tensor{\mathcal{{U}}}{}}_{\textbf{i},m}=\tensor{\mathcal{U}}{}_{\textbf{i},\textbf{j}{(m)}},\hskip
5.69046pt,m=1:L$ (65)

that is a basis for the space of tensors
$\tensor{\mathcal{X}}{}_{\textbf{i}}\in\mathbb{R}^{I_{1}\times
I_{2}\times\ldots\times I_{d}}$. As a consequence each tensor of data set
$\tensor{\mathcal{X}}{}$ can be represented as

$\tensor{\mathcal{X}}{}_{\textbf{i},n}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},{m}}d_{n,m},\hskip
5.69046ptn=1:N,\hskip 5.69046ptm=1:L$ (66)

where the $(N\times{L})$ matrix ${d_{n,m}=[D]_{n,m}}$ is given by the inner
product

$d_{n,m}=\tensor{\mathcal{X}}{}_{\textbf{i},n}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},{m}}$
(67)

The matrix $D$ can be decomposed by the SVD as

$D=Y\Sigma{Z}^{T}$ (68)

where

$\Sigma=\left(\begin{array}[]{cc}S&0\\\ 0&0\end{array}\right),\;\hskip
5.69046ptS=diag(\sigma_{1},\ldots,\sigma_{r})$ (69)

and $r=rank(D)$. The generic term of $D$ is

$d_{n,m}=Y_{n,l}S_{l,p}Z_{m,p},\hskip 5.69046ptl,p=1:r,$ (70)

then, from (66)

$\tensor{\mathcal{X}}{}_{\textbf{i},n}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},m}Y_{n,l}S_{l,p}Z_{m,p}$
(71)

and $S_{l,p}=\delta_{l,p}\sigma_{p}$, it results

$\tensor{\mathcal{X}}{}_{\textbf{i},n}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},m}\sigma_{l}Y_{n,l}Z_{m,l},\hskip
5.69046ptl=1:r$ (72)

where $\sigma_{l}Y_{n,l}Z_{m,l}$ is the rank-1 decomposition of matrix $D$,
and $\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},m}\in{R^{I\times{L}}}$ is the
rank-1 basis.

A pseudo-code of the procedure previously described to derive a rank-1 basis
is described in Algorithm 2.

Algorithm 2 Rank-1 basis

INPUT: dataset of tensors

$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},\hskip
5.69046ptn=1:N\\},\hskip 5.69046pt\textbf{i}=(i_{1},\ldots,i_{d}),\hskip
5.69046pt1\leq{i_{k}}\leq{I_{k}}$

for k=1:d

1\. Derive the $(I_{k}\times{I_{k}})$ symmetric matrix $A^{k}$ from tensor
product along $k$-mode.

$A^{k}(i_{k},j_{k})=\tensor{\mathcal{X}}{}(i_{1},\ldots,{i_{k}},\ldots,{i_{d}},n)\tensor{\mathcal{X}}{}(i_{1},\ldots,{j_{k}},\ldots,{i_{d}},n)$

2\. Compute matrix $U^{k}=U_{i_{k},j_{k}}$ solving eigenvalue equation

$A^{k}=U^{k}\Sigma^{k}U^{k}$

end

3.The basis as outer product of matrices $U_{i_{k},j_{k}}$

$\tensor{\mathcal{U}}{}_{\textbf{i,j}}={U}_{i_{1},j_{1}}{U}_{i_{2},j_{2}}\ldots{U}_{i_{d},j_{d}}$

4\. Transform the index i to linear index $m$

$\hat{\tensor{\mathcal{{U}}}{}}_{\textbf{i},m}=\tensor{\mathcal{U}}{}_{\textbf{i},\textbf{j}{(m)}},\hskip
5.69046pt,m=1:L$

5\. Derive the $(N\times{L})$ matrix $(D)_{n,m}=d_{n,m}$ as inner product

$d_{n,m}=\tensor{\mathcal{X}}{}_{\textbf{i},n}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},m}$

6\. Decompose the matrix $d_{n,m}$ by SVD

$d_{n,m}=Y_{n,l}S_{l,p}Z_{m,p}\hskip 5.69046ptl,p=1:rank(D)$

OUTPUT: orthogonal expansion of tensor $\tensor{\mathcal{X}}{}_{\textbf{i},n}$

$\tensor{\mathcal{X}}{}_{\textbf{i},n}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},{m}}\sigma_{l}Y_{n,l}Z_{m,l},\hskip
5.69046ptl=1:r$

#### Tensor PCA

A principal component analysis can be carried out as follows. As the
eigenvalues are nonnegative, the elements of $S_{l,p}$ are put in decreasing
order and the corresponding eigenvectors in $Y_{n,l}$ and $Z_{m,p}$ are
ordered accordingly. Then the principal components analysis corresponds to the
truncated representation

$\hat{\tensor{\mathcal{X}}{}}_{\textbf{i},n}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},{l}}d_{n,l},\hskip
5.69046ptl=1:M$ (73)

with $M<r$, that minimizes the error

$\tensor{\mathcal{E}}{}=\|\tensor{\mathcal{X}}{}_{\textbf{i},n}-{\hat{\tensor{\mathcal{X}}{}}_{\textbf{i},n}\|^{2}}=\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},{m}}d_{n,m}\hat{\tensor{\mathcal{U}}{}}_{\textbf{i},{m}}d_{n,m}=d_{n,m}d_{n,m}$
(74)

with $m=M+1:r$. Using (70) and orthonormality of $Y_{n,l}$ and $Z_{m,l}$ we
have

$\tensor{\mathcal{E}}{}=Y_{n,l}Y_{n,l^{\prime}}\sigma_{l}\sigma_{l^{\prime}}Z_{m,l}Z_{m,l^{\prime}}=\sigma_{l}\sigma_{l},\hskip
5.69046ptl=M+1:r.$ (75)

The error is minimum when the singular values $\sigma_{l}$ corresponding to
the basis elements in residual are the lowest, or when the singular values
corresponding to elements in truncated representation are the highest. Thus
(73) is an optimal subspace representation of dataset
$\tensor{\mathcal{X}}{}$.

### C. Basis of a subspace

Let us suppose a set

$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},\hskip
5.69046ptn=1:N\\}$ (76)

of N order-$d$ tensors $\tensor{\mathcal{X}}{}_{\textbf{i}},$ is given, with
$N<L=I_{1}I_{2}\ldots{I_{d}}$. Assuming $r$ tensors within this set are
linearly independent, thus this subset spans a subspace in tensor space
$\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{d}}$, whose dimension is
$r\leq{N}$. The problem of deriving a basis from this subset is equivalent to
the QR problem [76], that is well known in matrix algebra, which aims to
derive an orthonormal set of vectors, from a set of independent vectors. In
tensor space the problem can be formulated as follows. Assume every tensors in
$\tensor{\mathcal{X}}{}$ can be represented as linear combination of $r$
orthogonal tensors $\tensor{\mathcal{Q}}{}_{\textbf{i},l},\hskip
2.84544ptl=1:r$. Formally we have

$\tensor{\mathcal{X}}{}_{\textbf{i},n}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}r_{n,l},\hskip
5.69046ptl=1:r$ (77)

where $r_{n,l}$ is the $(N\times{r})$ matrix of coefficients in the linear
combination (77).

The tensor product of $\tensor{\mathcal{X}}{}_{\textbf{i},n}$ with itself
along the vector index i, gives

$\tensor{\mathcal{X}}{}_{\textbf{i},n}\tensor{\mathcal{X}}{}_{\textbf{i},m}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}r_{n,l}\tensor{\mathcal{Q}}{}_{\textbf{i},l^{\prime}}r_{m,l^{\prime}}=r_{n,l}r_{m,l^{\prime}}\delta_{l,l^{\prime}}=r_{n,l}r_{m,l},\hskip
5.69046pt$ (78)

where the orthogonality of $\tensor{\mathcal{Q}}{}_{\textbf{i},l}$, i.e.
$\tensor{\mathcal{Q}}{}_{\textbf{i},l}\tensor{\mathcal{Q}}{}_{\textbf{i},l^{\prime}}=\delta_{l,l^{\prime}}$,
is used. The tensor product (78) represents a real nonnegative self-adjoint
operator in the space $R^{N}$, due to symmetry of the $(N\times{N})$ matrix
$r_{n,l}r_{m,l}$, thus ensuring a set of nonnegative eigenvalues exists.

It can be easily shown that the rank of matrix $[R]_{n,l}=r_{n,l}$ is $r$. In
fact due to the orthogonality of $\tensor{\mathcal{Q}}{}_{\textbf{i},l}$, it
results

$\tensor{\mathcal{Q}}{}_{\textbf{i},l^{\prime}}\tensor{\mathcal{X}}{}_{\textbf{i},n}=\tensor{\mathcal{Q}}{}_{\textbf{i},l^{\prime}}\tensor{\mathcal{Q}}{}_{\textbf{i},l}r_{n,l}=\delta_{l,l^{\prime}}r_{n,l}=r_{n,l}$
(79)

and

$r_{n,l}r_{n,l^{\prime}}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}\tensor{\mathcal{X}}{}_{\textbf{i},n}\tensor{\mathcal{Q}}{}_{\textbf{i},l^{\prime}}\tensor{\mathcal{X}}{}_{\textbf{i},n}=$
$\big{\|}\tensor{\mathcal{X}}{}\big{\|}^{2}\tensor{\mathcal{Q}}{}_{\textbf{i},l}\tensor{\mathcal{Q}}{}_{\textbf{i},l^{\prime}}=\big{\|}\tensor{\mathcal{X}}{}\big{\|}^{2}\delta_{l,l^{\prime}}$
(80)

(80) shows that two generic columns $r_{n,l}$ and $r_{n,l^{\prime}}$ of matrix
$R$ are orthogonal each other and, as a consequence $r=rank(R)$.
Due to symmetry, the $(N\times{N})$ matrix $r_{n,l}r_{m,l}$ can be decomposed
as

$r_{n,l}r_{m,l}=U_{n,p}\Sigma^{2}_{p,q}U^{T}_{q,m}=U_{n,p}\Sigma^{2}_{p,q}U_{m,q}=$
$U_{n,r}\Sigma_{p,q}\Sigma_{p,q}U_{m,q}=U_{n,p}\Sigma_{p,l}(U_{m,q}\Sigma_{q,l})$
(81)

The matrix on the left $r_{n,l}r_{m,l}=[RR^{T}]_{n,m}$ is the Gram-matrix of
$R^{T}$, and it is well known from matrix algebra that $rank(RR^{T})=rank(R)$.
This implies that the matrix $\Sigma_{p,q}$ is an
$(r^{\prime}\times{r^{\prime}})$ diagonal matrix
$\Sigma=diag(\sigma_{1}\ldots{\sigma_{r}})$, whose elements $\sigma_{j}$ are
the singular values, with $r^{\prime}=rank(R^{T}R)=rank(R)=r$. It is worth to
notice that any other matrix $RV$, whith $V$ orthogonal, satisfies (81), thus
$RR^{T}$ can be decomposed in a non unique way as

$r_{n,l}=U_{n,p}\Sigma_{p,l}$ (82)

Combining (77) and (82) gives

$\tensor{\mathcal{X}}{}_{\textbf{i},n}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}U_{n,p}\Sigma_{p,l}$
(83)

and the tensor product of $\tensor{\mathcal{X}}{}$ with matrix $U$ yields

$\tensor{\mathcal{X}}{}_{\textbf{i},n}U_{n,m}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}\Sigma_{p,l}U_{n,p}U_{n,m}=$
$\tensor{\mathcal{Q}}{}_{\textbf{i},l}\Sigma_{p,l}\delta_{p,m}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}\Sigma_{l,m}$
(84)

due to the orthogonality of $U$. The matrix $\Sigma_{p,l}$ is diagonal full
rank, thus $\Sigma^{-1}_{p,l}$ exists and from (84)

$\tensor{\mathcal{X}}{}_{\textbf{i},n}U_{n,m}\Sigma^{-1}_{m,l}=\tensor{\mathcal{Q}}{}_{\textbf{i},l}$
(85)

or

$\tensor{\mathcal{Q}}{}_{\textbf{i},l}=\tensor{\mathcal{X}}{}_{\textbf{i},n}b_{n,l},\hskip
5.69046ptl=1:r$ (86)

where

$b_{n,l}=U_{n,m}\Sigma^{-1}_{m,l}$ (87)

is an $(N\times{N})$ matrix. (86) and (87), together with the decomposition of
$RR^{T}$ give the requested relationships for the basis of subspace
$\tensor{\mathcal{X}}{}$.

A pseudo-code of the procedure previously described to derive a basis of a
subspace is described in Algorithm 3.

Algorithm 3 Basis of subspace

INPUT: dataset of tensors

$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},\hskip
5.69046ptn=1:N\\}$

for k=1:d

1\. Derive the symmetric $(N\times{N})$ matrix $(A)_{m,n}=a_{m,n}$ as tensor
product

$a_{m,n}=\tensor{\mathcal{X}}{}_{\textbf{i},m}\tensor{\mathcal{X}}{}_{\textbf{i},n}$

2\. Decompose the matrix $a_{m,n}$ by solving eigenvalue equation

$a_{m,n}=U_{n,p}\Sigma^{2}_{p,q}U_{m,q}\hskip 5.69046ptp,q=1:rank(A)$

3\. Compute the matrix $b_{n,l}$

$b_{n,l}=U_{n,m}\Sigma^{-1}_{m,l}$

OUTPUT: The basis

$\tensor{\mathcal{Q}}{}_{\textbf{i},l}=\tensor{\mathcal{X}}{}_{\textbf{i},n}b_{n,l},\hskip
5.69046ptl=1:rank(A)$

#### Tensor PCA

To derive a principal component analysis, as the eigenvalues are nonnegative,
we put the elements of $\Sigma_{p,l}$ in decreasing order and we order the
corresponding eigenvectors in $U_{n,p}$ accordingly. Then by defining the
truncated matrices

$\hat{U}_{n,p}=U_{n,p},\hskip 5.69046pt\hat{\Sigma}_{n,p}=\Sigma_{n,p},\hskip
5.69046ptp,l=1:M$ (88)

with $M<r$, principal components analysis corresponds to the subspace
representation

$\hat{\tensor{\mathcal{X}}{}}_{\textbf{i},n}={{Q}}_{\textbf{i},l}\hat{U}_{n,p}\hat{\Sigma}_{p,l}\hskip
5.69046ptn=1:N,\hskip 5.69046ptp,l=1:M$ (89)

that minimizes the error

$\tensor{\mathcal{E}}{}=\|\tensor{\mathcal{X}}{}_{\textbf{i},n}-{\hat{\tensor{\mathcal{X}}{}}_{\textbf{i},n}\|^{2}}={{Q}}_{\textbf{i},l}\hat{U}_{n,p}\hat{\Sigma}_{p,l}{{Q}}_{\textbf{i},l^{\prime}}\hat{U}_{n,p^{\prime}}\hat{\Sigma}_{p^{\prime},l^{\prime}}$
(90)

Due to orthonormality of ${{Q}}_{\textbf{i},l}$ and $\hat{U}_{n,p}$ it results

$\tensor{\mathcal{E}}{}=\hat{\Sigma}_{p,l}\hat{\Sigma}_{p,l}=\sigma_{p}\sigma_{p}$
(91)

with $p=M+1:r$. Thus the error is minimum when the singular values
corresponding to the basis elements in the residual are the lowest, or when
the singular values corresponding to elements in truncated representation are
the highest. Thus (89) is an optimal subspace repesentation of database
$\tensor{\mathcal{X}}{}$.

##  5 Experimental results

The aim of the experiments is to validate the three approaches previously
described in Section 4 for tensor PCA. To this end we refer to the space of
order-3 tensors, in which the generic element represents an RGB image. The
experiments have been conducted on three different datasets, namely CIFAR-10,
CALTECH-101 and FLOWERS, to show that tensor PCA gives a low-dimensional
subspace representation of datasets.

###  5.1 Experiments on CIFAR-10

The first experiment aims to validate the approach discussed in Subsection
4.A, in which a basis is derived from a self-adjoint operator. In this
experiment the CIFAR-10 dataset has been used[77, 78]. The dataset consists of
60000, $32\times 32$ colour images divided in 10 classes, with 6000 images per
class, of which 50000 are for training and 10000 for testing.

The experiment refers to the Class-0 alone of dataset and 5000 RGB images of
the training set for image reconstruction by representation (34) have been
used. For this purpose, the data has been organized in a order-$4$ tensor
$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},n=1:N\\}$,
with $N=5000,\;\tensor{\mathcal{X}}{}\in\mathbb{R}^{32\times 32\times 3\times
5000}$, so that the dimension of the basis is $L=32\cdot 32\cdot 3=3072$. A
self-adjoint operator can be derived from the dataset $\tensor{\mathcal{X}}{}$
as tensor product

$\tensor{\mathcal{A}}{}_{\textbf{i},\textbf{j}}=\tensor{\mathcal{X}}{}_{\textbf{i},n}\tensor{\mathcal{X}}{}_{\textbf{j},n}$

since it results in
$\tensor{\mathcal{A}}{}_{\textbf{i},\textbf{j}}=\tensor{\mathcal{A}}{}_{\textbf{j},\textbf{i}}$,
thus satisfying the requirements for Proposition 1. Solving equation (40) by
Algorithm 1 a set of $L$ eigenvalues and eigentensors are achieved.

Fig. 1 shows the behavior of eigenvalues as obtained from CIFAR-10 dataset,
using Algorithm 1, while some basis elements from CIFAR-10 dataset
corresponding to the first eigenvalues are reported in Fig. 2. Fig. 3 reports
several images taken from the dataset. For comparison the same images
reconstructed by (34), with the basis achieved with Algorithm 1, are reported
in Fig. 4. As you can see the two datasets match perfectly, thus proving the
validity of the approach. On the basis of the behavior of eigenvalues depicted
in Fig. 1 tensor PCA has been applied retaining only the principal components,
i.e. the ones corresponding to the most significant eigenvalues. Fig. 5 shows
the same images as Fig. 3, reconstructed with the 300 most significant
components of the basis, thus confirming the correctness of approximation by
tensor PCA.

Figure 1: The behavior of eigenvalues as obtained from CIFAR-10 dataset Figure
2: Some basis elements from CIFAR-10 dataset corresponding to the highest
eigenvalues Figure 3: Several images taken from the dataset CIFAR-10 Figure 4:
The same images of Fig. 3 reconstructed by (34) Figure 5: The same images of
Fig. 3 reconstructed with truncated basis

###  5.2 Experiments on CALTECH-101

The second experiment refers to the Algorithm 2 discussed in Subsection 4.B,
in which a rank-$1$ basis is derived from a dataset. For this purpose
CALTECH-101 dataset [79], that contains images from 101 object categories
(e.g., “helicopter”, “elephant” and “chair” etc.) and a background category
that contains the images not from the 101 object categories, has been used.
For each object category, there are about 40 to 800 images, while most classes
have about 50 images. The resolution of the image is roughly about 300×200
pixels. In this experiment 2990 RGB images reduced to a single size of
${64\times 64\times 3}$ have been used for image reconstruction by the
representation (34). For this purpose, the data has been organized in an
order-$4$ tensor
$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},n=1:N\\}$,
with $N=2990,\;\tensor{\mathcal{X}}{}\in\mathbb{R}^{64\times 64\times 3\times
2990}$, so that the dimension of the basis is $L=64\cdot 64\cdot 3=12288$.

Fig. 6 reports several images taken from the dataset. For comparison the same
images reconstructed by (34), with the basis achieved with Algorithm 2, are
reported in Fig. 7. Also in this case these results show a perfect
reconstruction of dataset, thus proving the validity of the algorithm.

Fig. 8 shows the behavior of singular values as obtained from CALTECH-101
dataset, while some basis elements from CALTECH-101 dataset corresponding to
the first singular values are reported in Fig. 9.

A truncated representation of dataset has been achieved from (72) by retaining
only the components corresponding to the most significant singular values in
matrix $S$.

Fig. 10 shows the same images as Fig. 6, reconstructed with the 500 most
significant components of the basis.

Figure 6: Several images taken from the dataset CALTECH-101 Figure 7: The same
images of Fig. 6 reconstructed by (34) Figure 8: The behavior of singular
values as obtained from CALTECH-101 dataset Figure 9: some basis elements from
CALTECH-101 dataset corresponding to the first eigenvalues Figure 10: The same
images of Fig. 6 reconstructed with truncated basis

###  5.3 Experiments on FLOWERS

The objective of the last experiment is to validate the approach reported in
Subsection 4.C to derive a basis in a subspace. In this experiment the FLOWERS
dataset [75] that contains 4242 colour images of flowers has been used. The
pictures are divided into five classes: chamomile, tulip, rose, sunflower,
dandelion. For each class there are about 800 photos. Photos, which are not
high resolution (about 320x240 pixels), are reduced to a single size of
${128\times 128\times 3}$. In this experiment we use 3660 RGB images of the
set for image reconstruction by the representation (34), with the basis
achieved solving the eigenvalue problem discussed in Sect.4.C. For this
purpose, the data has been organized as a order-$4$ tensor
$\tensor{\mathcal{X}}{}=\\{\tensor{\mathcal{X}}{}_{\textbf{i},n},n=1:N\\}$,
with $N=3660,\;\tensor{\mathcal{X}}{}\in\mathbb{R}^{128\times 128\times
3\times 3660}$, so that the dimension of the basis is $L=128\cdot 128\cdot
3=49152$.

Fig. 11 shows the behavior of eigenvalues as obtained from FLOWERS dataset,
using Algorithm 3, while some basis elements from FLOWERS dataset
corresponding to the first eigenvalues are reported in Fig. 12. Fig. 13
reports several images taken from the dataset. For comparison the same images
reconstructed by (34), with the basis achieved with Algorithm 3, are reported
in Fig. 14, showing a perfect match. These results confirm the validity of the
Algorithm 3 to derive a basis of subspace. A truncated representation of
dataset has been achieved from (89), by retaining only the components
corresponding to the most significant eigenvalues in matrix $\hat{\Sigma}$.
Fig. 15 shows the same images as Fig. 13, reconstructed with the 300 most
significant components of the basis.

Figure 11: The behavior of eigenvalues as obtained from FLOWERS dataset Figure
12: Some basis elements from FLOWERS dataset corresponding to the first
eigenvalues Figure 13: Several images taken from the dataset FLOWERS Figure
14: The same images of Fig. 13reconstructed by (34) Figure 15: The same images
of Fig. 13 reconstructed with truncated basis

##  6 Conclusion

Tensor PCA is a multilinear extension of PCA, one of the most popular and
efficient subspace techniques used in matrix analysis to reduce the
dimensionality of a dataset. In this paper a mathematical framework for tensor
PCA has been developed by deriving a basis in tensor space and by projecting
data on a subspace. In particular the problem of deriving a basis from a self-
adjoint tensor operator has been studied and the equivalence between
eigenvalue equation for this operator and standard matrix eigenvalue equation
has been proven.

## References

  * [1] 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.
  * [2] 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.
  * [3] D. Tao, X. Li, W. Hu, S. Maybank, and X. Wu, “Supervised tensor learning,” in _Fifth IEEE International Conference on Data Mining (ICDM’05)_. IEEE, 2005, pp. 8–pp.
  * [4] 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.
  * [5] M. Hou and B. Chaib-draa, “Fast recursive low-rank tensor learning for regression.” in _IJCAI_ , 2017, pp. 1851–1857.
  * [6] S. Zafeiriou, “Discriminant nonnegative tensor factorization algorithms,” _IEEE Transactions on Neural Networks_ , vol. 20, no. 2, pp. 217–235, 2009\.
  * [7] 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.
  * [8] 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.
  * [9] 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\.
  * [10] S. Aja-Fernández, R. de Luis Garcia, D. Tao, and X. Li, _Tensors in image processing and computer vision_. Springer Science & Business Media, 2009.
  * [11] 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.
  * [12] 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.
  * [13] 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.
  * [14] 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.
  * [15] 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.
  * [16] 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.
  * [17] 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.
  * [18] 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.
  * [19] 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.
  * [20] Q. Song, H. Ge, J. Caverlee, and X. Hu, “Tensor completion algorithms in big data analytics,” _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , vol. 13, no. 1, pp. 1–48, 2019.
  * [21] 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.
  * [22] 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.
  * [23] W. W. Sun and L. Li, “Store: sparse tensor response regression and neuroimaging analysis,” _The Journal of Machine Learning Research_ , vol. 18, no. 1, pp. 4908–4944, 2017.
  * [24] 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.
  * [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.-K. Kim, S.-F. Wong, and R. Cipolla, “Tensor canonical correlation analysis for action classification,” in _2007 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2007, pp. 1–8.
  * [28] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM review_ , vol. 51, no. 3, pp. 455–500, 2009.
  * [29] 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.
  * [30] 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.
  * [31] 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.
  * [32] L. Tucker, “Implications of factor analysis of three way matrices formeasurement of change,” _Problems inmeasuring change. Madison, Wis.: University of WisconsinPress_ , 1963.
  * [33] R. A. Harshman _et al._ , “Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis,” 1970.
  * [34] 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.
  * [35] K. Inoue, “Generalized tensor pca and its applications to image analysis,” _Applied Matrix and Tensor Variate Data Analysis_ , pp. 51–71, 2016.
  * [36] 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.
  * [37] 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.
  * [38] 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.
  * [39] 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.
  * [40] 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.
  * [41] 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.
  * [42] 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.
  * [43] J. Ye, “Generalized low rank approximations of matrices,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 112.
  * [44] 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.
  * [45] R. Bro, “Parafac. tutorial and applications,” _Chemometrics and intelligent laboratory systems_ , vol. 38, no. 2, pp. 149–171, 1997.
  * [46] 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.
  * [47] 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\.
  * [48] 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.
  * [49] A. Cichocki, “Era of big data processing: A new approach via tensor networks and tensor decompositions,” _arXiv preprint arXiv:1403.2048_ , 2014.
  * [50] A. Cichocki, A.-H. Phan, Q. Zhao, N. Lee, I. Oseledets, M. Sugiyama, D. P. Mandic _et al._ , “Tensor networks for dimensionality reduction and large-scale optimization: Part 2 applications and future perspectives,” _Foundations and Trends® in Machine Learning_ , vol. 9, no. 6, pp. 431–673, 2017.
  * [51] R. Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” _Annals of physics_ , vol. 349, pp. 117–158, 2014.
  * [52] L. Grasedyck, “Hierarchical singular value decomposition of tensors,” _SIAM journal on matrix analysis and applications_ , vol. 31, no. 4, pp. 2029–2054, 2010.
  * [53] 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.
  * [54] 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.
  * [55] M. Chaghazardi and S. Aeron, “Sample, computation vs storage tradeoffs for classification using tensor subspace models,” _arXiv preprint arXiv:1706.05599_ , 2017.
  * [56] 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.
  * [57] W. Wang, V. Aggarwal, and S. Aeron, “Principal component analysis with tensor train subspace,” _Pattern Recognition Letters_ , vol. 122, pp. 86–91, 2019\.
  * [58] 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.
  * [59] 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.
  * [60] I. T. Jolliffe and J. Cadima, “Principal component analysis: a review and recent developments,” _Philosophical transactions of the royal society A: Mathematical, Physical and Engineering Sciences_ , vol. 374, no. 2065, p. 20150202, 2016.
  * [61] J. P. Cunningham and Z. Ghahramani, “Linear dimensionality reduction: Survey, insights, and generalizations,” _The Journal of Machine Learning Research_ , vol. 16, no. 1, pp. 2859–2900, 2015.
  * [62] K. I. Diamantaras and S. Y. Kung, _Principal component neural networks: theory and applications_. John Wiley & Sons, Inc., 1996.
  * [63] J. Shlens, “A tutorial on principal component analysis,” _arXiv preprint arXiv:1404.1100_ , 2014.
  * [64] H. Abdi and L. J. Williams, “Principal component analysis,” _Wiley interdisciplinary reviews: computational statistics_ , vol. 2, no. 4, pp. 433–459, 2010.
  * [65] P. Sanguansat, _Principal component analysis_. BoD–Books on Demand, 2012.
  * [66] P. R. Halmos, _Finite-dimensional vector spaces_. Courier Dover Publications, 2017.
  * [67] L. Qi, “Eigenvalues of a real supersymmetric tensor,” _Journal of Symbolic Computation_ , vol. 40, no. 6, pp. 1302–1324, 2005.
  * [68] J. Brachat, P. Comon, B. Mourrain, and E. Tsigaridas, “Symmetric tensor decomposition,” _Linear Algebra and its Applications_ , vol. 433, no. 11-12, pp. 1851–1872, 2010.
  * [69] C.-F. Cui, Y.-H. Dai, and J. Nie, “All real eigenvalues of symmetric tensors,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 35, no. 4, pp. 1582–1601, 2014.
  * [70] K.-C. Chang, K. Pearson, and T. Zhang, “On eigenvalue problems of real symmetric tensors,” _Journal of Mathematical Analysis and Applications_ , vol. 350, no. 1, pp. 416–422, 2009.
  * [71] H. Chen, M. Li, H. Yan, and G. Zhou, “Further results on eigenvalues of symmetric decomposable tensors from multilinear dynamical systems,” _Applied Mathematics Letters_ , vol. 129, p. 107980, 2022.
  * [72] L. Qi, “Eigenvalues and invariants of tensors,” _Journal of Mathematical Analysis and Applications_ , vol. 325, no. 2, pp. 1363–1377, 2007.
  * [73] ——, “The spectral theory of tensors (rough version),” _arXiv preprint arXiv:1201.3424_ , 2012.
  * [74] K. Chang, L. Qi, and T. Zhang, “A survey on the spectral theory of nonnegative tensors,” _Numerical Linear Algebra with Applications_ , vol. 20, no. 6, pp. 891–912, 2013.
  * [75] L.-H. Lim, “Singular values and eigenvalues of tensors: a variational approach,” in _1st IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2005._ IEEE, 2005, pp. 129–132.
  * [76] G. H. Golub and C. F. Van Loan, _Matrix computations_. JHU press, 2013.
  * [77] A. Krizhevsky, G. Hinton _et al._ , “Learning multiple layers of features from tiny images,” 2009.
  * [78] “The CIFAR-10 dataset,” [Online]. Available: https://www.cs.toronto.edu/~kriz/cifar.html.
  * [79] “CALTECH-101 dataset,” [Online]. Available: https://data.caltech.edu/records/mzrjq-6wc02.