1

Matrix Product State for Higher-Order Tensor
Compression and Classification

arXiv:1609.04541v1 [stat.ML] 15 Sep 2016

Johann A. Bengua1 , Ho N. Phien1 , Hoang D. Tuan1 and Minh N. Do2
Abstract—This paper introduces matrix product state (MPS)
decomposition as a new and systematic method to compress multidimensional data represented by higher-order tensors. It solves
two major bottlenecks in tensor compression: computation and
compression quality. Regardless of tensor order, MPS compresses
tensors to matrices of moderate dimension which can be used for
classification. Mainly based on a successive sequence of singular
value decompositions (SVD), MPS is quite simple to implement
and arrives at the global optimal matrix, bypassing local alternating optimization, which is not only computationally expensive
but cannot yield the global solution. Benchmark results show that
MPS can achieve better classification performance with favorable
computation cost compared to other tensor compression methods.
Index Terms—Higher-order tensor compression and classification, supervised learning, matrix product state (MPS), tensor
dimensionality reduction.

I. I NTRODUCTION

T

HERE is an increasing need to handle large multidimensional datasets that cannot efficiently be analyzed
or processed using modern day computers. Due to the curse
of dimensionality it is urgent to develop mathematical tools
which can evaluate information beyond the properties of large
matrices [1]. The essential goal is to reduce the dimensionality
of multidimensional data, represented by tensors, with a minimal information loss by compressing the original tensor space
to a lower-dimensional tensor space, also called the feature
space [1]. Tensor decomposition is the most natural tool to
enable such compressions [2].
Until recently, tensor compression is merely based on
Tucker decomposition (TD) [3], also known as higher-order
singular value decomposition (HOSVD) when orthogonality
constraints on factor matrices are imposed [4]. TD is also
an important tool for solving problems related to feature
extraction, feature selection and classification of large-scale
multidimensional datasets in various research fields. Its wellknown application in computer vision was introduced in [5] to
analyze some ensembles of facial images represented by fifthorder tensors. In data mining, the HOSVD was also applied
to identify handwritten digits [6]. In addition, the HOSVD
has been applied in neuroscience, pattern analysis, image
classification and signal processing [7], [8], [9]. The higherorder orthogonal iteration (HOOI) [10] is an alternating least
squares (ALS) for finding the TD approximation of a tensor.
1 Faculty of Engineering and Information Technology, University
of Technology Sydney, Ultimo, NSW 2007, Australia; Email:
johann.a.bengua@student.uts.edu.au,
ngocphien.ho@uts.edu.au,
tuan.hoang@uts.edu.au.
2 Department of Electrical and Computer Engineering and the Coordinated
Science Laboratory, University of Illinois at Urbana- Champaign, Urbana, IL
61801 USA; Email: minhdo@illinois.edu

Its application to independent component analysis (ICA) and
simultaneous matrix diagonalization was investigated in [11].
Another TD-based method is multilinear principal component
analysis (MPCA) [12], an extension of classical principal
component analysis (PCA), which is closely related to HOOI.
Meanwhile, TD suffers the following conceptual bottlenecks
in tensor compression:
• Computation. TD compresses an N th-order tensor in
tensor space RI1 ×I2 ×···IN of large dimension I =
Q
N
j=1 Ij to its N th-order core tensor in a tensor space
QN
∆1 ×∆2 ×···∆N
R
of smaller dimension Nf = j=1 ∆j by
using N factor matrices of size Ij × ∆j . Computation
of these N factor matrices is computationally intractable.
Instead, each factor matrix is alternatingly optimized with
all other N − 1 factor matrices held fixed, which is still
computationally expensive. Practical application of the
TD-based compression is normally limited to small-order
tensors.
• Compression quality. TD is an effective representation of
a tensor only when the dimension of its core Q
tensor is
N
fairly large [2]. Restricting dimension Nf = j=1 ∆j
to a moderate size for tensor classification results in
significant lossy compression, making TD-based compression a highly heuristic procedure for classification.
It is also almost impossible to tune ∆j ≤ Ij among
QN
j=1 ∆j ≤ N̄f for a prescribed N̄f to have a better
compression.
In this paper, we introduce the matrix product state (MPS)
decomposition [13], [14], [15], [16] as a new method to
compress tensors, which fundamentally circumvent all the
above bottlenecks of TD-based compression. Namely,
• Computation. The MPS decomposition is fundamentally
different from the TD in terms of its geometric structure as it is made up of local component tensors with
maximum order three. Consequently, using the MPS
decomposition for large higher-order tensors can potentially avoid the computational bottleneck of the TD and
related algorithms. Computation for orthogonal common
factors in MPS is based on successive SVDs without
any recursive local optimization procedure and is very
efficient with low-cost.
• Compression quality. MPS compresses N th-order tensors
to their core matrices of size RN1 ×N2 . The dimension
Nf = N1 N2 can be easily tuned to a moderate size with
minimum information loss by pre-positioning the core
matrix in the MPS decomposition.
MPS has been proposed and applied to study quantum

2

many-body systems with great success, prior to its introduction
to the mathematics community under the name tensor-train
(TT) decomposition [17]. However, to the best of our knowledge its application to machine learning and pattern analysis
has not been proposed.
Our main contribution is summarized as follows:
•

•

•

Propose MPS decomposition as a new and systematic
method for compressing tensors of arbitrary order to
matrices of moderate dimension, which circumvents all
existing bottlenecks in tensor compression;
Develop MPS decomposition tailored for optimizing the
dimensionality of the core matrices and the compression
quality. Implementation issues of paramount importance
for practical computation are discussed in detail. These
include tensor mode permutation, tensor bond dimension
control, and positioning the core matrix in MPS;
Extensive experiments are performed along with comparisons to existing state-of-the-art tensor classification
methods to show its advantage.

A preliminary result of this work was presented in [18]. In
the present paper, we rigorously introduce the MPS as a new
and systematic approach to tensor compression for classification, with computational complexity and efficiency analysis.
Furthermore, new datasets as well a new experimental design
showcasing computational time and classification success rate
(CSR) benchmarks are included.
The rest of the paper is structured as follows. Section II
provides a rigorous mathematical analysis comparing MPS and
TD in the context of tensor compression. Section III is devoted
to MPS tailored for effective tensor compression, which also
includes a computational complexity analysis comparing MPS
to HOOI, MPCA and uncorrelated multilinear discriminant
analysis with regularization (R-UMLDA) [19]. In Section IV,
experimental results are shown to benchmark all algorithms in
classification performance and training time. Lastly, Section V
concludes the paper.
II. MPS DECOMPOSITION VS TD DECOMPOSITION IN
TENSOR COMPRESSION

We introduce some notations and preliminaries of multilinear algebra [2]. Zero-order tensors are scalars and denoted by
lowercase letters, e.g., x. A first-order tensor is a vector and
denoted by boldface lowercase letters, e.g., x. A matrix is a
second-order tensor and denoted by boldface capital letters,
e.g., X. A higher-order tensor (tensors of order three and
above) are denoted by boldface calligraphic letters, e.g., X .
Therefore, a general Nth-order tensor of size I1 ×I2 ×· · ·×IN
can be defined as X ∈ RI1 ×I2 ×···×IN , where each Ii is the
dimension of its mode i. We also denote xi , xij and xi1 ···iN
as the ith entry x(i), (i, j)th entry X(i, j) and (i1 , · · · , iN )th
entry X (i1 , · · · , iN ) of vector x, matrix X and higher-order
tensor X , respectively.
Mode-n matricization (also known as mode-n unfolding or flattening) of X is the process of Qunfolding or
reshaping X into a matrix X(n) ∈ RIn ×( i6=n Ii ) such
that X(n) (in , j) = X (i1 , · · · , in , · · · , iN ) for j = 1 +

Qk−1
− 1) m=1,m6=n Im . We also define the dimenQN
sion of X as n=1 In . The mode-n product of X with a
matrix A ∈ RJn ×In is denoted as X ×n A, which is a N thorder tensor of size I1 × · · · × In−1 × Jn × In+1 × · · · × IN
such that
(X ×n A)(i1 , · · · , in−1 , jn , in+1 , · · · , iN ) =
In
X
X (i1 , · · · , in , · · · iN )A(jn , in ).
PN

k=1,k6=n (ik

in =1

The
norm
PI of X is defined as ||X ||F =
PI
PI Frobenius
( i11=1 i22=1 · · · iNN =1 x2i1 i2 ···iN )1/2 .
We are concerned with the following problem of tensor
compression for supervised learning :
Based on K training Nth-order tensors X (k) ∈
I1 ×I2 ×···×IN
R
(k = 1, 2, . . . , K), find common factors to
compress both training tensor X (k) and test tensors Y (`)
(` = 1, · · · , L) to a feature space of moderate dimension to
enable classification.
Until now, only TD has been proposed to address this
problem [7]. More specifically, the K training sample tensors
are firstly concatenated along the mode (N + 1) to form an
(N + 1)th-order tensor X as
X = [X (1) X (2) · · · X (K) ] ∈ RI1 ×I2 ×···×IN ×K .

(1)

TD-based compression such as HOOI [10] is then applied to
have the approximation
X

≈ R ×1 U(1) ×2 U(2) · · · ×N U(N ) ,
(j)

(2)

Ij ×∆j

where each matrix U
∈ R
(j = 1, 2, . . . , N ) is
orthogonal, i.e. U(j)T U(j) = I (I ∈ R∆j ×∆j denotes the
identity matrix). It is called a common factor matrix and can
be thought of as the principal components in each mode j.
The parameters ∆j satisfying
∆j ≤ rank(X(j) )

(3)

are referred to as the compression ranks of the TD.
The (N + 1)th-order core tensor R and common factor
matrices U(j) ∈ RIj ×∆j are supposed to be found from the
following nonlinear least squares
min

R∈R∆1 ×···×∆N ×K ,
U(j) ∈RIj ×∆j ,j=1,...,N
(j) T

subject to (U
(1)

Φ(R, U(1) , · · · , U(N ) )

) U

(4)
(j)

= I, j = 1, ..., N,

(N )

where Φ(R, U , · · · , U ) := ||X − R ×1 U(1) ×2
U(2) · · · ×N U(N ) ||2F The optimization problem (4) is computationally intractable, which could be addressed only by
alternating least squares (ALS) in each U(j) (with other U(`) ,
` 6= j held fixed) [10]:
min

R(j) ∈R∆1 ×···×∆N ×K ,
U(j) ∈RIj ×∆j

Φ(j) (R(j) , U(j) )

subject to (U

(5)
(j) T

) U

(j)

= I,

where Φ(j) (R(j) , U(j) ) := ||X − R(j) ×1 U(1) ×2 U(2) · · · ×N
U(N ) ||2F . The computation complexity per one iteration consisting of N ALS (5) is [20, p. 127]
O(K∆I N + N KI∆2(N −1) + N K∆3(N −1) )

(6)

3

for
Ij ≡ I

and

∆j ≡ ∆, j = 1, 2, ..., N.

The optimal (N +1)th-order core tensor R ∈ R

(7)

∆1 ×···×∆N ×K

e (k) ∈
in (4) is seen as the concatenation of compressed X
R∆1 ×···×∆N of the sample tensors X (k) ∈ RI1 ×···×IN , k =
1, · · · , K:
e
R = [X

(1)

e (2) · · · X
e (N ) ] = X ×1 (U(1) )T · · · ×N (U(N ) )T .
X

Accordingly, the test tensors Y (`) are compressed to
e (`) = Y (`) ×1 (U(1) )T · · · ×N (U(N ) )T ∈ R∆1 ×···×∆N . (9)
Y
The number
Nf =

N
Y

∆j

(10)

j=1

I1 ≥ · · · ≥ In−1 and In ≤ · · · ... ≤ IN . The elements of
X can be presented in the following mixed-canonical form
[23] of the matrix product state (MPS) or tensor train (TT)
decomposition [16], [14], [15], [17]:
xi1 ···k···iN

si ,

(11)

where rj := rank(X (j) ) and {srj , · · · , s1 } is the set of nonzero eigenvalues of the positive definite matrix X (j) (X
)T
Q (j)
Ij ×(K `6=j I` )
in decreasing order. Since the matrix X (j) ∈ R
is highly unbalanced as a result of tensor matricization along
one mode versus the rest, it is almost full-row (low) rank
(rj ≈ Ij ) and its squared X (j) (X (j) )T of size Ij × Ij is
well-conditioned in the sense that its eigenvalues do not decay
quickly. As a consequence, (11) cannot be small for small ∆j
so the ALS (5) cannot result in a good approximation. The
information loss with the least square (5) is thus more than
X
i=1

si
P rj

i=1 si

(n−1)

(n)

(n+1)

≈ Bi1 · · · Bin−1 Gk Cin

(N +1)

· · · CiN

,

(15)
(j)

(j)

where matrices Bij and Ci(j−1) (the upper index “(j)”
denotes the position j of the matrix in the chain) of size
∆j−1 × ∆j (∆0 = ∆N +1 = 1), are called “left” and “right”
common factors which satisfy the following orthogonality
conditions:
Ij
X

(j)

(j)

(Bij )T Bij

= I,

(j = 1, . . . , n − 1)

(16)

ij =1

Ij−1

X

i=1

−

(k)

xi1 ···in ···iN

and

rj −∆j −1

rj −∆j −1

=

(1)

thus represents the dimension of the feature space
R∆1 ×···×∆N .
Putting aside the computational intractability of the optimal
factor matrices U(j) in (4), the TD-based tensor compression
by (8) and (9) is a systematic procedure only when the
right hand side of (2) provides a good approximation of X ,
which is impossible for small ∆j satisfying (3) [2]. In other
words, the compression of large dimensional tensors to small
dimensional tensors results in substantial lossy compression
under the TD framework. Furthermore, one can see the value
of (5) is lower bounded by
X

of Y (`) is defined by (Y (`) − X̄ ) ×1 (U (1) )T · · · ×N (U (N ) )T .
Thus, MPCA suffers the similar conceptual drawbacks
inQN
herent by TD. Particularly, restricting Nf =
j=1 ∆j to
a moderate size leads to ignoring many important principle
components.
We now present a novel approach to extract tensor features,
which is based on MPS. Firstly, permute all modes of the
tensor X and position mode K such that such that
(8)
(14)
X ∈ RI1 ×···In−1 ×K×In ···×IN ,

si
log2 Prj

i=1 si

,

(12)

which is really essential in the von Neumann entropy [21] of
X (j) :
rj
X
si
si
P rj
−
log2 Prj
.
(13)
s
i=1 i
i=1 si
i=1

(j)

(j)

Cij−1 (Cij−1 )T

=

I, (j = n + 1, . . . , N + 1)
(17)

ij−1 =1

respectively, where I denotes the identity matrix. Each matrix
(n)
Gk of dimension ∆n−1 × ∆n is the compression of the
training tensor X (k) . The parameters ∆j are called the bond
dimensions or compression ranks of the MPS. Using the com(j)
(j)
mon factors Bij and Ci(j−1) , we can extract the core matrices
for the test tensors Y (`) as follows. We permute all Y (`) ,
` = 1, · · · , L to ensure the compatibility between the training
(n)
and test tensors. The compressed matrix Q` ∈ R∆n−1 ×∆n
of the test tensor Y (`) is then given by
X
(n)
(1)
(n−1)
(`)
Q`
=
(Bi1 )T · · · (Bin−1 )T yi1 ······iN
i1 ,...,iN
(n+1) T

(Cin

(N +1) T

) · · · (CiN

) .

(18)

The dimension
Nf = ∆n−1 ∆n

(19)

Note that each entropy (13) quantifies only local correlation
between mode j and the rest [22]. The MPCA [12] aims at
(4) with

is the number of reduced features.

X = [(X (1) − X̄ ) · · · (X (K) − X̄ )]
PK
PL
1
with X̄ = K+L
( k=1 X (k) + `=1 Y `) ). With such definition of X , (N + 1)th-order core tensor X is the concatenation
of principal components of X (k) , while principal components

The advantage of MPS for tensor compression is that the
order N of a tensor does not affect directly the feature number
Nf in Eq. (19), which is only determined strictly by the
product of the aforementioned bond dimensions ∆n−1 and
∆n . In order to keep ∆n−1 and ∆n to a moderate size, it

III. TAILORED MPS FOR TENSOR COMPRESSION

4

is important to control the bond dimensions ∆j , and also to
optimize the positions of tensor modes as we address in this
section. In what follows, for a matrix X we denote X(i, :)
(X(:, j), resp.) as its ith row (jth column, resp.), while for
a third-order tensor X we denote X (:, `, :) as a matrix such
that its (i1 , i3 )th entry is X (i1 , `, i3 ). For a N th-order tensor
X ∈ RI1 ×···×IN we denote X[j] ∈ R(I1 I2 ···Ij )×(Ij+1 ···K···IN )
as its mode-(1, 2, . . . , j) matricization. It is obvious that
X[1] = X(1) .
A. Adaptive bond dimension control in MPS
To decompose the training tensor X into the MPS according
to Eq. (15), we apply two successive sequences of SVDs to
the tensor which include left-to-right sweep for computing the
(1)
(n−1)
left common factors Bi1 , . . . , Bin−1 , and right-to-left sweep
(n+1)
(N +1)
Cin , . . . , CiN

which satisfy the left-canonical constraint due to (26).
(3)
Applying the same procedure for determining Bi3 by reshap(2)
∆2 ×(I3 ···K···IN )
ing the matrix V ∈ R
to
W(3) ∈ R(∆2 I3 )×(I4 ···K···IN ) ,
performing the SVD, and so on. This procedure is iterated till
obtaining the last QR-approximation
W(n−1)

≈ U(n−1) V(n−1) ∈ R(∆n−2 In−1 )×(KIn ···IN ) ,
U(n−1) ∈ R(∆n−2 In−1 )×∆n−1 ,
V(n−1) ∈ R∆n−1 ×(KIn ···IN ) ,
(28)
with U(n−1) orthogonal:
U(n−1) (U(n−1) )T = I

(29)

and reshaping U(n−1) ∈ R(∆n−2 In−1 )×∆n−1 into a third-order
tensor U ∈ R∆n−2 ×In−1 ×∆n−1 to define the last left common
factors

for computing the right common factors
(n)
and the core matrix Gk in Eq. (15) as follows:
• Left-to-right sweep for left factor computation:
(j)
The left-to-right sweep involves acquiring matrices Bij
(ij = 1, . . . , Ij ; j = 1, . . . , n − 1) fulfilling orthogonality
condition in Eq. (16). Start by performing the mode-1 matricization of X to obtain

Bin−1 = U (:, in−1 , :) ∈ R∆n−2 ×∆n−1 , in−1 = 1, · · · , In−1 ,
(30)
which satisfy the left-canonical constraint due to (29).
In a nutshell, after completing the left-to-right sweep, the
elements of tensor X are approximated by

W(1) := X [1] = X (1) ∈ RI1 ×(I2 ···K···IN ) .

xi1 ···in−1 in ···iN +1 ≈ Bi1 · · · Bin−1 V(n−1) (:, kin · · · iN ). (31)

(k)

For
∆1 ≤ rank(X[1] ),

(20)

apply SVD to W(1) to have the QR-approximation
W(1)

≈

U(1) V(1) ∈ RI1 ×(I2 ···K···IN ) ,

(21)

where U(1) ∈ RI1 ×∆1 is orthogonal:
(U(1) )T U(1) = I,

(22)

and V(1) ∈ R∆1 ×(I2 ···K···IN ) . Define the the most left common
factors by
(1)

Bi1 = U(1) (i1 , :) ∈ R1×∆1 , i1 = 1, · · · , I1

(23)

which satisfy the left-canonical constraint in Eq. (16) due to
(22).
Next, reshape the matrix V(1) ∈ R∆1 ×(I2 ···K···IN ) to W(2) ∈
R(∆1 I2 )×(I3 ···K···IN ) . For
∆2 ≤ rank(W

(2)

) ≤ rank(X[2] ),

(24)

apply SVD to W(2) for the QR-approximation
W(2) ≈ U(2) V(2) ∈ R(∆1 I2 )×(I3 ···K···IN ) ,

(25)

(26)

(2)

(27)

(n−1)

The matrix V(n−1) ∈ R∆n−1 ×(KIn ···IN ) is reshaped to
W(N ) ∈ R(∆n−1 K···IN −1 )×IN for the next right-to-left sweeping process.
• Right-to-left sweep for right factor computation:
Similar to left-to-right sweep, we perform a sequence of
SVDs starting from the right to the left of the MPS to get the
(j)
matrices Cij−1 (ij−1 = 1, . . . , Ij−1 ; j = N + 1, . . . , n + 1)
fulfilling the right-canonical condition in Eq. (17). To start,
we apply the SVD to the matrix W(N ) ∈ R(∆n−1 K···IN −1 )×IN
obtained previously in the left-to-right sweep to have the RQapproximation
W(N ) ≈ U(N ) V(N ) ,
(32)
where U(N ) ∈ R(∆n−1 K···IN −1 )×∆N and V(N ) ∈ R∆N ×IN is
orthogonal:
V(N ) (V(N ) )T = I
(33)
for
∆N ≤ rank(W(N ) ) ≤ rank(X [N −1] ).

(34)

Define the most right common factors
(N +1)

and V(2) ∈ R∆2 ×(I3 ···K···IN ) . Reshape the matrix U(2) ∈
R(∆1 I2 )×∆2 into a third-order tensor U ∈ R∆1 ×I2 ×∆2 to
define the next common factors
Bi2 = U (:, i2 , :) ∈ R∆1 ×∆2 , i2 = 1, · · · , I2 ,

(1)

CiN

where U(2) ∈ R(∆1 I2 )×∆2 is orthogonal such that
(U(2) )T U(2) = I,

(n−1)

= V(N ) (:, iN ) ∈ R∆N ×1 , iN = 1, · · · , IN ,

which satisfy the right-canonical constraint (17) due to (33).
Next, reshape U(N ) ∈ R(∆n−1 K···IN −1 )×∆N into W(N −1) ∈
R(∆n−1 K···IN −2 )×(IN −1 ∆N ) and apply the SVD to have the
RQ-approximation
W(N −1) ≈ U(N −1) V(N −1) ,
(N −1)

(∆n−1 K···IN −2 )×∆N −1

where U
∈ R
R∆N −1 ×(IN −1 ∆N ) is orthogonal:

V(N −1) (V(N −1) )T = I

(35)
and V

(N −1)

∈

(36)

5

Algorithm I: MPS for tensor feature extraction

for
∆N −1 ≤ rank(W

(N −1)

)≤ rank(X [N −2] ).

(37)

Reshape the matrix V(N −1) ∈ R∆N −1 ×(IN −1 ∆N ) into a
third-order tensor V ∈ R∆N −1 ×IN −1 ×∆N to define the next
common factor
(N )

CiN −1 = V(:, iN −1 , :) ∈ R∆N −1 ×∆N

(38)

which satisfy Eq. (17) due to (36).
This procedure is iterated till obtaining the last RQapproximation
W(n)

≈ U(n) V(n) ∈ R(∆n−1 K)×(In ∆n+1 ) ,
U(n) ∈ R(∆n−1 K)×∆n , V(n) ∈ R∆n ×(In ∆n+1 ) ,
(39)
with V(n) orthogonal:
V(n) (V(n) )T = I

(40)

∆n ≤ rank(W(n) )≤ rank(X [n−1] ).

(41)

for

Reshape V(n) ∈ R(∆n )×(In ∆n+1 ) into a third-order tensor V ∈
R∆n ×In ×∆n+1 to define the last right common factors
(n+1)

Ci n

= V(:, in , :) ∈ R∆n−1 ×∆n , in = 1, · · · , In ,

(42)

which satisfy (17) due to (40).
By reshaping U(n) ∈ R(∆n−1 K)×∆n into a third-order
(n)
tensor G ∈ R∆n−1 ×K×∆n to define Gk
= G(:, k, :),
k = 1, · · · .K, we arrive at Eq. (15).
Note that the MPS decomposition described by Eq. (15)
can be performed exactly or approximately depending on the
bond dimensions ∆j (j = 1, . . . , N ). The bond dimension
truncation is of crucial importance to control the final feature
number Nf = ∆n−1 ∆n . To this end, we rely on thresholding
the singular values of W(j) . With a threshold  being defined
in advance, we control ∆j such that ∆j largest singular values
s1 ≥ s2 ≥ ... ≥ s∆j satisfy
P∆j
si
Pri=1
≥ ,
j
i=1 sj

(43)

(j)

for rj = rank(W ). The information loss from the von
Neumann entropy (13) of W(j) by this truncation is given
by (12). The entropy of each W(j) provides the correlation
degree between two sets of modes 1, · · · , j and j + 1, · · · , N
[22]. Therefore, the N entropies W(j) , j = 1, · · · , N provide
the mean of the tensor’s global correlation. Furthermore, rank
rj of each W(j) is upper bounded by
min {I1 · · · Ij , Ij+1 · · · IN }

(44)

making the truncation (43) highly favorable in term of compression loss to matrices of higher rank due to balanced row
and column numbers.
A detailed outline of our MPS approach to tensor feature
extraction is presented in Algorithm 1.

Input:
Output:

X ∈ RI1 ×···×In−1 ×K···×IN ,
: SVD threshold
(n)
Gk ∈ R∆n−1 ×∆n , k = 1, · · · , K
(j)
Bij (ij = 1, . . . , Ij , j = 1, . . . , n − 1)
(j)

Ci

(j−1)

(i(j−1) = 1, . . . , I(j−1) , j = n + 1, . . . , N + 1)

1: Set W(1) = X(1)
% Mode-1 matricization of X
2: for j = 1 to n − 1
% Left-to-right sweep
3:
W(j) = USV
% SVD of W(j)
4:
Wj ≈ U(j) W(j+1)
% Thresholding S for QR-approximation
5:
Reshape U(j) to U
(j)
6:
Bij = U (:, ij , :)
% Set common factors
7: end
8: Reshape V(n−1) to WN ∈ R(∆n−1 K···IN )×IN
9: for j = N down to n % right-to-left sweep
10:
W(j) = USV
% SVD of W(j)
11:
W(j) ≈ W(j−1) V(j)
% Thresholding S for RQ-approximation
13:
Reshape V(j) to V
(j+1)
% Set common factors
14:
Cij−1 = V(:, ij−1 , :)
15: end
16: Reshape U(n) into G ∈ R∆n−1 ×K×∆n
(n)
17: Set Gk = G(:, k, :)
% Training core matrix
Texts after symbol “%” are comments.

B. Tensor mode pre-permutation and pre-positioning mode K
for MPS
One can see from (44) that the efficiency of controlling the
bond dimension ∆j is dependent on its upper bound (44).
Particularly, the efficiency of controlling the bond dimensions
∆n−1 and ∆n that define the feature number (19) is dependent
on
min {I1 · · · In−1 , In · · · IN }
(45)
Therefore, it is important to pre-permute the tensors modes
such that the ratio
Qn−1 QN
min{ i=1 Ii , i=n Ii }
(46)
Qn−1 QN
max{ i=1 Ii , i=n Ii }
is near to 1 as possible, while {I1 , · · · , In−1 } is in decreasing
order
I1 ≥ · · · ≥ In−1
(47)
and {In , · · · , IN } in increasing order
In ≤ · · · ≤ IN

(48)

Qj
QN
min{ i=1 Ij , i=j+1 Ii }
Qj
QN
max{ i=1 Ij , i=j+1 Ii }

(49)

to improve the ratio

for balancing W(j) .
The mode K is then pre-positioned in n-th mode as in (14).
C. Complexity analysis
In the following complexity analysis it is assumed In = I
∀n for simplicity. The dominant computational complexity
of MPS is O(KI (N +1) ) due to the first SVD of the matrix
obtained from the mode-1 matricization of X . On the other
hand, the computational complexity of HOOI requires several

6

iterations of an ALS method to obtain convergence. In addition, it usually employs the HOSVD to initialize the tensors
which involves the cost of order O(N KI N +1 ), and thus very
expensive with large N compared to MPS.
MPCA is computationally upper bounded by O(N KI N +1 ),
however, unlike HOOI, MPCA doesn’t require the formation
of the (N + 1)th order core tensor at every iteration and
convergence can usually happen in one iteration [12].1
The computational
complexity of R-UMLDA is approxiPN
mately O(K n=2 I n + (C + K)I 2 + (p − 1)[IK + 2I 2 +
(p − 1)2 + (2I(p − 1)] + 4I 3 ), where C is the number of
classes, p is the number of projections, which determines the
core vector size [19]. Therefore, R-UMLDA would perform
poorly for many samples and classes.
D. MPS-based tensor object classification
This subsection presents two methods for tensor objection
classification based on Algorithm 1. For each method, an
explanation of how to reduce the dimensionality of tensors
to core matrices, and subsequently to feature vectors for
application to linear classifiers is given.
1) Principal component analysis via tensor-train (TTPCA):
The TTPCA algorithm is an approach where Algorithm 1 is
applied directly on the training set, with no preprocessing such
as data centering. Specifically, given a set of N th-order tensor
samples X (k) ∈ RI1 ×I2 ×···×IN , then the core matrices are
obtained as
(n)

Gk

∈ R∆n−1 ×∆n .

(50)

Vectorizing each k sample results in
(n)

gk

∈ R∆n−1 ∆n .

(51)

Using (43), ∆n−1 ∆nQ features of k is significantly less in
N
comparison to Nf = n=1 In of X (k) , which allows for PCA
to be easily applied, followed by a linear classifier.
2) MPS: The second algorithm is simply called MPS,
where in this case we first perform data centering on the set
of training samples {X (k) }, then apply Algorithm 1 to obtain
the core matrices
(n)

Gk

∈ R∆n−1 ×∆n .

(52)

Vectorizing the K samples results in (51), and subsequent
linear classifiers such as LDA or nearest neighbors can be
utilized. In this method, MPS can be considered a multidimensional analogue to PCA because the tensor samples have
been data centered and are projected to a new orthogonal space
using Algorithm 1, resulting in the core matrices.
IV. E XPERIMENTAL RESULTS
In this section, we conduct experiments on the proposed
TTPCA and MPS algorithms for tensor object classification.
An extensive comparison is conducted based on CSR and
training time with tensor-based methods MPCA, HOOI, and
R-UMLDA.
1 This does not mean that MPCA is computationally efficient but in contrast
this means that alternating iterations of MPCA prematurely terminate, yielding
a solution that is far from the optimal one of a NP-hard problem.

Four datasets are utilized for the experiment. The Columbia
Object Image Libraries (COIL-100) [24], [25], Extended Yale
Face Database B (EYFB) [26], BCI Jiaotong dataset (BCI)
[27], and the University of South Florida HumanID “gait
challenge“ dataset (GAIT) version 1.7 [28] . All simulations
are conducted in a Matlab environment.
A. Parameter selection
TTPCA, MPA and HOOI rely on the threshold  defined in
(43) to reduce the dimensionality of a tensor, while keeping its
most relevant features. To demonstrate how the classification
success rate (CSR) varies, we utilize different  for each
dataset. It is trivial to see that a larger  would result in a longer
training time due to its computational complexity, which was
discussed in subsection III-C. Furthermore, TTPCA utilizes
PCA, and a range of principal components p is used for the
experiments. HOOI is implemented with a maximum of 10
ALS iterations. MPCA relies on fixing an initial quality factor
Q, which is determined through numerical simulations, and a
specified number of elementary multilinear projections (EMP),
we denote as mp , must be initialized prior to using the RUMLDA algorithm. A range of EMP’s is determined through
numerical simulations and the regularization parameter is fixed
to γ = 10−6 .
B. Tensor object classification
1) COIL-100: For this dataset we strictly compare MPS
and the HOSVD-based algorithm HOOI to analyse how adjusting  affects the approximation of the original tensors, as well
as the reliability of the extracted features for classification.
The COIL-100 dataset has 7200 color images of 100 objects
(72 images per object) with different reflectance and complex
geometric characteristics. Each image is initially a 3rd-order
tensor of dimension 128 × 128 × 3 and then is downsampled
to the one of dimension 32 × 32 × 3. The dataset is divided
into training and test sets randomly consisting of K and L
images, respectively according to a certain holdout (H/O) ratio
L
. Hence, the training and test sets are represented
r, i.e. r = K
by four-order tensors of dimensions 32 × 32 × 3 × K and
32 × 32 × 3 × L, respectively. In Fig. 1 we show how a few
objects of the training set (r = 0.5 is chosen) change after
compression by MPS and HOOI with two different values of
threshold,  = 0.9, 0.65. We can see that with  = 0.9, the
images are not modified significantly due to the fact that many
features are preserved. However, in the case that  = 0.65, the
images are blurred. That is because fewer features are kept.
However, we can observe that the shapes of objects are still
preserved. Especially, in most cases MPS seems to preserve
the color of the images better than HOOI. This is because
the bond dimension corresponding to the color mode I3 = 3
has a small value, e.g. ∆3 = 1 for  = 0.65 in HOOI. This
problem arises due to the the unbalanced matricization of the
tensor corresponding to the color mode. Specifically, if we
take a mode-3 matricization of tensor X ∈ R32×32×3×K , the
resulting matrix of size 3 × (1024K) is extremely unbalanced.
Therefore, when taking SVD with some small threshold , the
information corresponding to this color mode may be lost due

7

TABLE I: COIL-100 classification results. The best CSR corresponding to different H/O ratios obtained by MPS and HOOI.
Algorithm
r = 50%
HOOI
MPS
r = 90%
HOOI
MPS

CSR

Nf



98.87 ± 0.19
99.19 ± 0.19

198
120

0.80
0.80

87.22 ± 0.56
89.38 ± 0.40

112
59 ± 5

0.75
0.75

to dimension reduction. On the contrary, we can efficiently
avoid this problem in MPS by permuting the tensor such that
X ∈ R32×K×3×32 before applying the tensor decomposition.
(a) Original samples (size 32 × 32 × 3)

(b) Samples with MPS, ǫ = 0.9 (core size 18 × 24)

(c) Samples with HOOI, ǫ = 0.9 (core size 18 × 16 × 2)

(d) Samples with MPS, ǫ = 0.65 (core size 6 × 3)

(e) Samples with HOOI, ǫ = 0.65 (core size 5 × 4 × 1)

Fig. 1: Modification of ten objects in the training set of COIL100 are shown after applying MPS and HOOI corresponding
to  = 0.9 and 0.65 to compress tensor objects.

100
95

99
CSR(%)

98
r=50%

97

r=80%

90

96
95
MPS
HOOI

94

85

93
85
90
CSR(%)

80
85
75
r=90%

80
75
70

r=95%

70
65

0.6

0.7

0.8
ε

0.9

60

0.6

0.7

0.8

0.9

ε

Fig. 2: Error bar plots of CSR versus thresholding rate  for
different H/O ratios.
K nearest neighbors with K=1 (KNN-1) is used for classification. For each H/O ratio, the CSR is averaged over
10 iterations of randomly splitting the dataset into training
and test sets. Comparison of performance between MPS and

CSR
r = 80%
94.13 ± 0.42
95.37 ± 0.31
r = 95%
77.76 ± 0.90
83.17 ± 1.07

Nf



112
18

0.75
0.65

112
18

0.75
0.65

HOOI is shown in Fig. 2 for four different H/O ratios, i.e.
r = (50%, 80%, 90%, 95%). In each plot, we show the CSR
with respect to threshold . We can see that MPS performs
quite well when compared to HOOI. Especially, with small ,
MPS performs much better than HOOI. Besides, we also show
the best CSR corresponding to each H/O ratio obtained by
different methods in Table. I. It can be seen that MPS always
gives better results than HOOI even in the case of small value
of  and number of features Nf defined by (10) and (19) for
HOOI and MPS, respectively.
2) Extended Yale Face Database B: The EYFB dataset contains 16128 grayscale images with 28 human subjects, under 9
poses, where for each pose there is 64 illumination conditions.
Similar to [29], to improve computational time each image
was cropped to keep only the center area containing the face,
then resized to 73 x 55. The training and test datasets are not
selected randomly but partitioned according to poses. More
precisely, the training and test datasets are selected to contain
poses 0, 2, 4, 6 and 8 and 1, 3, 5, and 7, respectively. For a
single subject the training tensor has size 5 × 73 × 55 × 64 and
4 × 73 × 55 × 64 is the size of the test tensor. Hence for all 28
subjects we have fourth-order tensors of sizes 140×73×55×64
and 112 × 73 × 55 × 64 for the training and test datasets,
respectively.
In this experiment, the core tensors remains very large even
with a small threshold used, e.g., for  = 0.75, the core
size of each sample obtained by TTPCA/MPS and HOOI are
18 × 201 = 3618 and 14 × 15 × 13 = 2730, respectively,
because of slowly decaying singular values, which make them
too large for classification. Therefore, we need to further reduce the sizes of core tensors before feeding them to classifiers
for a better performance. In our experiment, we simply apply
a further truncation to each core tensor by keeping the first
few dimensions of each mode of the tensor. Intuitively, this
can be done as we have already known that the space of each
mode is orthogonal and ordered in such a way that the first
dimension corresponds to the largest singular value, the second
one corresponds to the second largest singular value and so on.
Subsequently, we can independently truncate the dimension
of each mode to a reasonably small value (which can be
determined empirically) without changing significantly the
meaning of the core tensors. It then gives rise to core tensors of
smaller size that can be used directly for classification. More
specifically, suppose that the core tensors obtained by MPS and
HOOI have sizes Q × ∆1 × ∆2 and Q × ∆1 × ∆2 × ∆3 , where
Q is the number K (L) of training (test) samples, respectively.
˜1 × ∆
˜ 2 and
The core tensors are then truncated to be Q × ∆
˜ 1 ×∆
˜ 2 ×∆
˜ 3 , respectively such that ∆
˜ l < ∆l (l = 1, 2, 3).
Q×∆
˜ l is chosen to be the same for both training
Note that each ∆

8

TABLE II: EYFB classification results
Algorithm
KNN-1
HOOI
MPS
TTPCA
MPCA
R-UMLDA
LDA
HOOI
MPS
TTPCA
MPCA
R-UMLDA

CSR ( = 0.9)

CSR ( = 0.85)

CSR ( = 0.80)

CSR ( = 0.75)

90.71 ± 1.49
94.29 ± 0.49
86.05 ± 0.44
90.89 ± 1.32
71.34 ± 2.86

90.89 ± 1.60
94.29 ± 0.49
86.01 ± 0.86

91.61 ± 1.26
94.29 ± 0.49
87.33 ± 0.46

88.57 ± 0.80
94.29 ± 0.49
86.99 ± 0.53

96.07 ± 0.80
97.32 ± 0.89
95.15 ± 0.45
90.00 ± 2.92
73.38 ± 1.78

95.89 ± 0.49
97.32 ± 0.89
95.15 ± 0.45

96.07 ± 0.49
97.32 ± 0.89
95.15 ± 0.45

96.07 ± 0.49
97.32 ± 0.89
94.86 ± 0.74

TABLE III: BCI Jiaotong classification results
Algorithm
Subject 1
HOOI
MPS
TTPCA
MPCA
R-UMLDA
CSP
Subject 2
HOOI
MPS
TTPCA
MPCA
R-UMLDA
CSP
Subject 3
HOOI
MPS
TTPCA
MPCA
R-UMLDA
CSP
Subject 4
HOOI
MPS
TTPCA
MPCA
R-UMLDA
CSP
Subject 5
HOOI
MPS
TTPCA
MPCA
R-UMLDA
CSP

CSR ( = 0.9)

CSR ( = 0.85)

CSR ( = 0.80)

CSR ( = 0.75)

84.39 ± 1.12
87.24 ± 1.20
78.57 ± 3.95
82.14 ± 3.50
63.18 ± 0.37
80.14 ± 3.73

83.37 ± 0.99
87.55 ± 1.48
78.43 ± 3.73

82.04 ± 1.05
87.24 ± 1.39
79.43 ± 4.12

84.80 ± 2.21
87.65 ± 1.58
79.14 ± 2.78

83.16 ± 1.74
90.10 ± 1.12
80.57 ± 0.93
81.29 ± 0.78
70.06 ± 0.39
81.71 ± 8.96

82.35 ± 1.92
90.10 ± 1.12
81.14 ± 1.86

82.55 ± 1.93
90.00 ± 1.09
81.29 ± 1.78

79.39 ± 1.62
91.02 ± 0.70
80 ± 2.20

60.92 ± 1.83
61.12 ± 1.36
67.43 ± 2.56
56.14 ± 2.40
57.86 ± 0.00
77.14 ± 2.26

61.84 ± 1.97
61.22 ± 1.53
68.29 ± 2.56

61.12 ± 1.84
61.12 ± 1.54
67.71 ± 2.28

60.51 ± 1.47
60.71 ± 1.54
66.43 ± 2.02

48.27 ± 1.54
52.35 ± 2.82
50.29 ± 2.97
51.00 ± 3.96
46.36 ± 0.93
59.86 ± 1.98

47.55 ± 1.36
52.55 ± 3.40
49.71 ± 3.77

49.98 ± 1.29
52.55 ± 3.69
49.14 ± 3.48

47.96 ± 1.27
51.84 ± 3.11
52.00 ± 3.48

60.31 ± 1.08
59.39 ± 2.08
53.43 ± 2.79
50.43 ± 1.48
55.00 ± 0.55
59.14 ± 2.11

60.82 ± 0.96
59.18 ± 2.20
54.29 ± 3.19

59.90 ± 2.20
58.57 ± 1.60
53.86 ± 3.83

60.41 ± 1.36
59.29 ± 1.17
54.86 ± 2.49

and test core tensors. In regards to TTPCA, each core matrix
is vectorized to have ∆1 ∆2 features, then PCA is applied.
Classification results for different threshold values  is
shown in Table. II for TTPCA, MPS and HOOI using
two different classifiers, i.e. KNN-1 and LDA. Results from
MPCA and R-UMLDA is also included. The core tensors
obtained by MPS and HOOI are reduced to have sizes of
˜1 × ∆
˜ 2 and Q × ∆
˜1 × ∆
˜2 × ∆
˜ 3 , respectively such
Q×∆
˜1 = ∆
˜ 2 = ∆ ∈ (10, 11, 12, 13, 14) and ∆
˜ 3 = 1.
that ∆
Therefore, the reduced core tensors obtained by both methods
have the same size for classification. With MPS and HOOI,
each value of CSR in Table. II is computed by taking the
average of the ones obtained from classifying different reduced
core tensors due to different ∆. In regards to TTPCA, for each
, a range of principal components p = {50, . . . , 70} is used.

We utilize Q = {70, 75, 80, 85, 90} for MPCA, and the range
mp = {10, . . . , 20} for R-UMLDA. The average CSR’s are
computed with TTPCA, MCPA and R-UMLDA according to
their respective range of parameters in Table. II. We can see
that the MPS gives rise to better results for all threshold values
using different classifiers. More importantly, MPS with the
smallest  can produce the highest CSR. The LDA classifier
gives rise to the best result, i.e. 97.32 ± 0.89.
3) BCI Jiaotong: The BCIJ dataset consists of single trial
recognition for BCI electroencephalogram (EEG) data involving left/right motor imagery (MI) movements. The dataset
includes five subjects and the paradigm required subjects to
control a cursor by imagining the movements of their right
or left hand for 2 seconds with a 4 second break between
trials. Subjects were required to sit and relax on a chair,

9

TABLE IV: GAIT classification results
Algorithm
Probe A
HOOI
MPS
TTPCA
MPCA
R-UMLDA
Probe C
HOOI
MPS
TTPCA
MPCA
R-UMLDA
Probe D
HOOI
MPS
TTPCA
MPCA
R-UMLDA
Probe F
HOOI
MPS
TTPCA
MPCA
R-UMLDA

CSR ( = 0.9)

CSR ( = 0.85)

CSR ( = 0.80)

CSR ( = 0.75)

63.71 ± 3.36
70.03 ± 0.42
75.31 ± 0.29
55.77 ± 1.08
46.62 ± 2.13

63.90 ± 3.40
70.03 ± 0.38
76.03 ± 0.38

64.16 ± 3.39
70.01 ± 0.36
76.38 ± 0.78

64.33 ± 3.20
69.99 ± 0.38
77.75 ± 0.92

36.67 ± 2.84
41.46 ± 0.64
39.17 ± 0.90
29.35 ± 2.29
20.87 ± 0.76

36.73 ± 2.79
41.36 ± 0.64
40.83 ± 0.41

36.70 ± 3.07
41.29 ± 0.63
41.61 ± 1.02

36.87 ± 3.68
41.46 ± 0.59
44.40 ± 1.54

19.73 ± 0.91
23.82 ± 0.42
21.92 ± 0.54
21.11 ± 3.43
7.88 ± 1.00

19.96 ± 1.15
23.84 ± 0.43
22.14 ± 0.20

20.32 ± 0.93
23.84 ± 0.45
22.84 ± 0.42

20.29 ± 1.11
23.84 ± 0.40
21.92 ± 0.59

20.77 ± 0.92
20.50 ± 0.40
14.78 ± 0.60
17.12 ± 2.79
9.67 ± 0.58

20.71 ± 0.72
20.52 ± 0.34
14.74 ± 0.77

20.15 ± 0.65
20.50 ± 0.29
15.29 ± 0.75

19.96 ± 0.67
20.56 ± 0.46
15.40 ± 0.55

TABLE V: Seven experiments in the USF GAIT dataset
Probe set
Size
Differences

A(GAL)
71
View

B(GBR)
41
Shoe

C(GBL)
41
Shoe, view

D(CAR)
70
Surface

looking at a computer monitor approximately 1m from the
subject at eye level. For each subject, data was collected over
two sessions with a 15 minute break in between. The first
session contained 60 trials (30 trials for left, 30 trials for right)
and were used for training. The second session consisted of
140 trials (70 trials for left, 70 trials for right). The EEG
signals were sampled at 500Hz and preprocessed with a filter
at 8-30Hz, hence for each subject the data consisted of a
multidimensional tensor channel × time × Q. The common
spatial patterns (CSP) algorithm [30] is a popular method
for BCI classification that works directly on this tensor, and
provides a baseline for the proposed and existing tensor-based
methods. For the tensor-based methods, we preprocess the data
by transforming the tensor into the time-frequency domain
using complex Mortlet wavelets with bandwidth parameter
fb = 6Hz (CMOR6-1) to make classification easier [31],
[32]. The wavelet center frequency fc = 1Hz is chosen.
Hence, the size of the concatenated tensors are 62 channels×
23 f requency bins × 50 time f rames × Q.
We perform the experiment for all subjects. After applying
the feature extraction methods MPS and HOOI, the core
tensors still have high dimension, so we need to further reduce
their sizes before using them for classification. For instance,
the reduced core sizes of MPS and HOOI are chosen to be
Q × 12 × ∆ and Q × 12 × ∆ × 1, where ∆ ∈ (8, . . . , 14),
respectively. With TTPCA, the principal components p =
{10, 50, 100, 150, 200}, Q = {70, 75, 80, 85, 90} for MPCA
and mp = {10, . . . , 20} for R-UMLDA. With CSP, we average
CSR for a range of spatial components sc = {2, 4, 6, 8, 10}.
The LDA classifier is utilized and the results are shown in

E(CBR)
44
Surface, shoe

F(CAL)
70
Surface, view

G(CBL)
44
Surface, view, shoe

Table. III for different threshold values of TTPCA, MPS and
HOOI. The results of MPCA, R-UMLDA and CSP are also
included. MPS outperforms the other methods for Subjects 1
and 2, and is comparable to HOOI in the results for Subject
5. CSP has the highest CSR for Subjects 3 and 4, followed by
MPS or TTPCA, which demonstrates the proposed methods
being effective at reducing tensors to relevant features, more
precisely than current tensor-based methods.

Fig. 3: The gait silhouette sequence for a third-order tensor.
4) USF GAIT challenge: The USFG database consists of
452 sequences from 74 subjects who walk in elliptical paths in
front of a camera. There are three conditions for each subject:
shoe type (two types), viewpoint (left or right), and the surface
type (grass or concrete). A gallery set (training set) contains
71 subjects and there are seven types of experiments known as

10

probe sets (test sets) that are designed for human identification.
The capturing conditions for the probe sets is summarized in
Table V, where G, C, A, B, L and R stand for grass surface,
cement surface, shoe type A, shoe type B, left view and right
view, respectively. The conditions in which the gallery set was
captured is grass surface, shoe type A and right view (GAR).
The subjects in the probe and gallery sets are unique and there
are no common sequences between the gallery and probe sets.
Each sequence is of size 128 × 88 and the time mode is 20,
hence each gait sample is a third-order tensor of size 128×88×
20, as shown in Fig. 3. The gallery set contains 731 samples,
therefore the training tensor is of size 128×88×20×731. The
test set is of size 128 × 88 × 20 × Ps , where Ps is the sample
size for the probe set that is used for a benchmark, refer to
Table V. The difficulty of the classification task increases with
the amount and and type of variables, e.g. Probe A only has
the viewpoint, whereas Probe F has surface and viewpoint,
which is more difficult. For the experiment we perform tensor
object classification with Probes A, C, D and F (test sets).
The classification results based on using the LDA classifier
is shown in Table IV. The threshold  still retains many
features in the core tensors of MPS and HOOI. Therefore,
further reduction of the core tensors is chosen to be Q×20×∆
and Q × 20 × ∆ × 1, where ∆ ∈ (8, . . . , 14), respectively.
The principal components for TTPCA is the range p =
{150, 200, 250, 300}, Q = {70, 75, 80, 85} for MPCA and
mp = {10, . . . , 20} for R-UMLDA. The proposed algorithms
achieve the highest performance for Probes A, C, and D. MPS
and HOOI are similar for the most difficult test set Probe F.
C. Training time benchmark
An additional experiment on training time for MPS2 , HOOI,
MPCA and R-UMLDA is provided to understand the computational complexity of the algorithms. For the COIL-100
dataset, we measure the elapsed training time for the training
tensor of size 32 × 32 × 3 × K (K = 720, 3600, 6480)
for H/O= {0.9, 0.5, 0.1}, according to 10 random partitions
of train/test data (iterations). MPCA, HOOI and R-UMLDA
reduces the tensor to 32 features, and MPS to 36 (due to
a fixed dimension ∆2 ). In Fig. 4a, we can see that as the
number of training images increases, the MPS algorithms
computational time only slightly increases, while MCPA and
HOOI increases gradually, with UMLDA having the slowest
performance overall.
The EYFB benchmark reduces the training tensor features to
36 (for MPS), 32 (MPCA and HOOI), and 16 (UMLDA, since
the elapsed time for 32 features is too long). For this case, Fig.
4b demonstrates that MPCA provides the fastest computation
time due to its advantage with small sample sizes (SSS). MPS
performs the next best, followed by HOOI, then UMLDA with
the slowest performance.
The BCI experiment involves reducing the training tensor to
36 (MPS) or 32 (MPS, HOOI and UMLDA) features and the
elapsed time is shown for Subject 1 in Fig. 4c. For this case
MPS performs the quickest compared to the other algorithms,
with UMLDA again performing the slowest.
2 TTPCA

would be equivalent in this experiment.

Lastly, the USFG benchmark tests Probe A by reducing
the MPS training tensor to 36 features, MPCA and HOOI
to 32 features, and UMLDA to 16 features. Fig. 4d shows
that MPCA provides the quickest time to extract the features,
followed by MPS, HOOI and lastly UMLDA.
V. C ONCLUSION
In this paper, a rigorous analysis of MPS and Tucker
decomposition proves the efficiency of MPS in terms of
retaining relevant correlations and features, which can be used
directly for tensor object classification. Subsequently, two new
approaches to tensor dimensionality reduction based on compressing tensors to matrices are proposed. One method reduces
a tensor to a matrix, which then utilizes PCA. And the other
is a new multidimensional analogue to PCA known as MPS.
Furthermore, a comprehensive discussion on the practical
implementation of the MPS-based approach is provided, which
emphasizes tensor mode permutation, tensor bond dimension
control, and core matrix positioning. Numerical simulations
demonstrates the efficiency of the MPS-based algorithms
against other popular tensor algorithms for dimensionality
reduction and tensor object classification.
For the future outlook, we plan to explore this approach
to many other problems in multilinear data compression and
tensor super-resolution.
R EFERENCES
[1] H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos, “A survey of
multilinear subspace learning for tensor data,” Pattern Recognition,
vol. 44, no. 7, pp. 1540 – 1551, 2011.
[2] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,”
SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.
[3] L. Tucker, “Some mathematical notes on three-mode factor analysis,”
Psychometrika, vol. 31, no. 3, pp. 279–311, 1966.
[4] L. D. Lathauwer, B. D. Moor, and J. Vandewalle, “A multilinear
singular value decomposition,” SIAM Journal on Matrix Analysis and
Applications, vol. 21, no. 4, 2000.
[5] M. A. O. Vasilescu and D. Terzopoulos, “Multilinear analysis of image
ensembles: Tensorfaces,” Proceedings of the 7th European Conference
on Computer Vision, Lecture Notes in Comput. Sci., vol. 2350, pp. 447–
460, 2002.
[6] B. Savas and L. Eldén, “Handwritten digit classification using higher
order singular value decomposition,” Pattern Recognition, vol. 40, no. 3,
pp. 993 – 1003, 2007.
[7] A. H. Phan and A. Cichocki, “Tensor decompositions for feature extraction and classification of high dimensional datasets,” IEICE Nonlinear
Theory and Its Applications, vol. 1, no. 1, pp. 37–68, 2010.
[8] L. Kuang, F. Hao, L. Yang, M. Lin, C. Luo, and G. Min, “A tensor-based
approach for big data representation and dimensionality reduction,”
IEEE Trans. Emerging Topics in Computing, vol. 2, no. 3, pp. 280–
291, Sept 2014.
[9] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa,
and A. H. 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, March 2015.
[10] L. D. Lathauwer, B. D. Moor, and J. Vandewalle, “On the best rank-1
and rank-(r1,r2,. . .,rn) approximation of higher-order tensors,” SIAM J.
Matrix Anal. Appl., vol. 21, no. 4, pp. 1324–1342, Mar. 2000.
[11] L. D. Lathauwer and J. Vandewalle, “Dimensionality reduction in higherorder signal processing and rank-(R1,R2,,RN) reduction in multilinear
algebra,” Linear Algebra and its Applications, vol. 391, no. 0, pp. 31 –
55, 2004.
[12] H. Lu, K. Plataniotis, and A. Venetsanopoulos, “MPCA: Multilinear
principal component analysis of tensor objects,” IEEE Trans. Neural
Networks, vol. 19, no. 1, pp. 18–39, Jan 2008.
[13] F. Verstraete, D. Porras, and J. I. Cirac, “Density matrix renormalization group and periodic boundary conditions: A quantum information
perspective,” Phys. Rev. Lett., vol. 93, no. 22, p. 227205, Nov 2004.

11

H/O=0.9

H/O=0.5

50

H/O=0.1

50

50

50

45

45

40

40

35

35

35

30

30

30

25

25

25

20

20

20

15

15

15

15

10

10

10

10

5

5

5

5

40

MPCA
UMLDA
MPS
HOOI

45

MPCA
UMLDA
MPS
HOOI

45

40

0

0
0

5

10

30
25
20

0

0
0

Iteration

Time (s)

Time (s)

35

5

10

0

5

Iteration

1

10

2

3

4

5

6

7

8

9

10

Iteration

Iteration

(a) COIL-100.

(b) EYFB.

30

200
MPCA
UMLDA
MPS
HOOI

180
25
160
140

Time (s)

Time (s)

20

15
MPCA
UMLDA
MPS
HOOI

10

120
100
80
60
40

5
20
0

0
1

2

3

4

5

6

7

8

9

10

1

Iteration

2

3

4

5

6

7

8

9

10

Iteration

(c) BCI Subject 1.

(d) GAIT Probe A.

Fig. 4: Training time of datasets for MPS, MPCA, HOOI and R-UMLDA.

[14] G. Vidal, “Efficient classical simulation of slightly entangled quantum
computation,” Phys. Rev. Lett., vol. 91, no. 14, p. 147902, Oct 2003.
[15] ——, “Efficient simulation of one-dimensional quantum many-body
systems,” Phys. Rev. Lett., vol. 93, no. 4, p. 040502, Jul 2004.
[16] D. Pérez-Garcı́a, F. Verstraete, M. Wolf, and J. Cirac, “Matrix product
state representations,” Quantum Information and Computation, vol. 7,
no. 5, pp. 401–430, 2007.
[17] I. V. Oseledets, “Tensor-train decomposition,” SIAM Journal on Scientific Computing, vol. 33, no. 5, pp. 2295–2317, 2011.
[18] J. A. Bengua, H. N. Phien, and H. D. Tuan, “Optimal feature extraction
and classification of tensors via matrix product state decomposition,” in
2015 IEEE International Congress on Big Data, June 2015, pp. 669–
672.
[19] H. Lu, K. N. Plataniotis, and A. N. Venetsanopoulos, “Uncorrelated
multilinear discriminant analysis with regularization and aggregation
for tensor object recognition,” IEEE Transactions on Neural Networks,
vol. 20, no. 1, pp. 103–123, Jan 2009.
[20] M. Ishteva, P.-A. Absil, S. van Huffel, and L. de Lathauwer, “Best low
multilinear rank approximation of higher-order tensors, based on the
Riemannian trust-region scheme.” SIAM J. Matrix Anal. Appl., vol. 32,
no. 1, pp. 115–135, 2011.
[21] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum
information. Cambridge, England: Cambridge University Press, 2000.
[22] J. A. Bengua, H. N. Phien, H. D. Tuan, and M. N. Do,
“Efficient tensor completion for color image and video recovery:
Low-rank tensor train,” arXiv preprint, 2016. [Online]. Available:
http://arxiv.org/abs/1606.01500
[23] U. Schollwöck, “The density-matrix renormalization group in the age

[24]
[25]

[26]

[27]
[28]

[29]

[30]

[31]

[32]

of matrix product states,” Annals of Physics, vol. 326, no. 1, pp. 96 –
192, 2011.
S. A. Nene, S. K. Nayar, and H. Murase, “Columbia object image library
(coil-100),” Technical Report CUCS-005-96, Feb 1996.
M. Pontil and A. Verri, “Support vector machines for 3d object recognition,” IEEE Trans. Patt. Anal. and Mach. Intell., vol. 20, no. 6, pp.
637–646, Jun 1998.
A. Georghiades, P. Belhumeur, and D. Kriegman, “From few to many:
illumination cone models for face recognition under variable lighting
and pose,” IEEE Trans. Pattern Anal. and Mach. Intell., vol. 23, no. 6,
pp. 643–660, Jun 2001.
(2013) Data set for single trial 64-channels eeg classification in bci.
[Online]. Available: http://bcmi.sjtu.edu.cn/resource.html
S. Sarkar, P. J. Phillips, Z. Liu, I. R. Vega, P. Grother, and K. W.
Bowyer, “The humanid gait challenge problem: data sets, performance,
and analysis,” IEEE Transactions on Pattern Analysis and Machine
Intelligence, vol. 27, no. 2, pp. 162–177, Feb 2005.
Q. Li and D. Schonfeld, “Multilinear discriminant analysis for higherorder tensor data classification,” IEEE Trans. Patt. Anal. and Mach.
Intell., vol. 36, no. 12, pp. 2524–2537, Dec 2014.
Y. Wang, S. Gao, and X. Gao, “Common spatial pattern method for
channel selelction in motor imagery based brain-computer interface,”
in 2005 IEEE Engineering in Medicine and Biology 27th Annual
Conference, Jan 2005, pp. 5392–5395.
Q. Zhao and L. Zhang, “Temporal and spatial features of singletrial eeg for brain-computer interface,” Computational Intelligence and
Neuroscience, vol. 2007, pp. 1–14, Jun 2007.
A. H. Phan, “NFEA: Tensor toolbox for feature extraction and applica-

12

tion,” Lab for Advanced Brain Signal Processing, BSI, RIKEN, Tech.
Rep., 2011.

