Multilinear Algebra andTensor Decomposition
Qibin Zhao
Tensor Learning Unit
RIKEN AIP

2018-6-2 @ Waseda UniversitySelf-Introduction• 2009, Ph.D. in Computer Science, Shanghai Jiao
Tong University• 2009 - 2017, RIKEN Brain Science Institute• 2017 - Now, RIKEN AIP• Research Interests:• Brain computer interface, brain signal processing• Tensor decomposition and machine learningSelf-Introduction

Brain computer interface

Self-Introduction
260

Introduction and Motivation

4th-order tensor=

RECRUITMENT. . . INFORMATIONMARCH 15, 2017

RIKEN Center for AdvancedIntelligence Project (AIP)

5th-order tensors
...

...

...

...

...

=

http://www.riken.jp/en/research/labs/aip/=

Tensor Learning Unit
Unit Leader:
6th-order
tensor Dr. Qibin Zhao

Tensors are high-dimensional generalizations ofvectors and matrices, which can provide a naturaland scalable representation for multi-dimensional,= multi-aspect data with inherent
multi-relational or
structure and complex dependence. In our team, weinvestigate the various tensor-based machinelearning models, e.g., tensor decomposition,
multilinear latent variable model, tensor regression
and representationsclassification, tensor
tensor learning,
andtensors.
BayesianFigure 1.6: Graphicalandnetworks,
symbols deepfor higher-orderblock
tensorlearning,
with
aim
to
facilitate
thelearningfromhighEach block represents either a 3rd-order tensor or a 2nd-order tensor. The outerdimensional
structured
large-scaleparameter
space. aIn
circle indicates a global
structureof the data
blockortensor(e.g.latenta vector,a matrix,
addition,
we innerdevelopthe reflects
scalabletheandefficientof tensor
learning
3rd-order block tensor),while thecirclestructureeach element
algorithms
supported inbythe
the top
theoreticalprinciples,the goal
within the block tensor.
For example,
diagrama vectorwith
of 3rdordertoadvanceexisting
machinelearning
approaches.
The
novel
applications
tensors is represented by an outer circle with one edge (a vector) which surrounds
in computer
brain
data analysis
be exploited
an inner circle with threeedges vision(a 3rdandorder
tensor),so thatwill
thealsowholestructureto
providenew insights into tensor learning methods.
designates a 4th-ordertensor.

Opening Positions
be visualized through changes in the architecture of a tensor network
diagram.

1

2

3

1.3

POSTDOCTORAL andTECHNICALRESEARCH ofCurse of Dimensionality
Generalized Separation
Variables forRESEARCHER
Multivariate Functions
STAFF
INTERN

1.3.1

Curse of Dimensionality

Doctoral degree

Technical support
Ph.D students areresearchers
preferable
The term curse of dimensionality wasforcoinedby Bellman(1961) toindicate that the number of samples needed to estimate an arbitrary
function with a given
levelof accuracy
with theWe areseekingtalentedgrows
and exponentiallycreative researchers
who

are willing to solve the challenging problems in machinelearning. For research topics, please refer to the bottomright side. If you are interested in joining our team,please contact us (see the top-right side).

Contact Information
Mitsui Building, 15th floor,
1-4-1 Nihonbashi, Chuoku, Tokyo103-0027, Japan
Email: qibin.zhao@riken.jp

Research Field
Computer Science
Related Fields
Machine Learning
Computer Vision
Neuroscience
Research Subjects
TensorDecomposition
Tensor Networks
Tensor Regression
and Classification
Deep Tensor
Learning
Bayesian Tensor
Learning

Outline• Vector and linear algebra• Matrix and its decomposition• What is tensor?
• Basic operations in tensor algebra• Classical tensor decomposition✦

CP Decomposition✦

Tucker Decomposition

Vectorsn
Warm-up: Vectors in R

We can think of vectors in two ways:
Points in a multidimensional space with respect to somecoordinate system translation of a point in a multidimensional spaceex., translation of the origin (0,0)
Dot Product or Scalar ProductDot product or scalar product

Dot product is the product of two vectors
Example:

x⋅y =



x1
x2

⋅

y1
y2

= x1 y1 + x2 y2 = s

It is the projection of one vector onto anotherx ⋅ y = x y cos θθx.yDot Product or Scalar Product

Scalar ProductCommutative:x⋅y = y ⋅xDistributive:(x + y )⋅ z = x ⋅ z + y ⋅ zLinearity(cx ) ⋅ y = c(x ⋅ y )x ⋅ (cy ) = c(x ⋅ y )(c1x) ⋅ (c2 y ) = (c1c2 )(x ⋅ y )

Non-negativity:

NormsNorms in RnEuclidean norm (sometimes called 2-norm):x = x 2 = x ⋅ x = x + x +L+ x =
2
1

2
2

2
n

n

∑x
i =1

2
i



The length of a vector is defined to be its (Euclidean) norm.

A unit vector is of length 1.

Non-negativity properties also hold for the norm:

Singular
Value
Decomposition
SingularValue
Decomposition:
T
D=USV

=
S

U

D
I1x I 2

A matrix D∈ R

SVD orthogonalizes these spaces and decomposes D

has a column space and a row space( U contains the left singular vectors/eigenvectors )( V contains the right singular vectors/eigenvectors )D = USV T

V

Rewrite as a sum of a minimum number of rank-1 matrices
RD=∑ σ u ov
r =1

r

r

r

Matrix SVD Properties

D=USV

Matrix SVD Properties:

Rank Decomposition:
 sum of min. number of rank-1 matrices

D=

σ

+σ

1

v1

2

T

u1



D = ∑ σ u ov
R

r =1

…..

v

T
2

u2

Multilinear Rank Decomposition:

D

U

+σ

R

r

vRT

uR
R1

R2

D = ∑ ∑ σ u ov
r1 =1

=

r

S

r2 =1

r1 r 2

V

r1

r2

r

T

Matrix in Machine
Learning
Matrix factorisation
modelsData often available in matrix form.

featurescoefficient

samplesMatrix inmodels
Machine Learning
Matrix factorisation
Data often available in matrix form.

moviesmovie
rating4users

Matrix in Machine
Learning
Matrix factorisation
modelsData often available in matrix form.

words

word
count
57

text documents

Decomposition
MatrixMatrix
factorisation
models in Machine Learning⇡

dictionary learning
low-rank approximation
factor analysis
latent semantic analysis
data X

dictionary W≈

activations H

MatrixMatrix
factorisation
models
Decompositionin Machine Learning⇡

dictionary learning
low-rank approximation
factor analysis
latent semantic analysis
data X

dictionary W

activations H≈

8

Matrix
Decomposition
Matrixfactorisation
models in Machine Learningfor dimensionality reduction (coding, low-dimensional embedding)≈Decomposition
MatrixMatrix
factorisation
models in Machine Learningfor interpolation (collaborative filtering, image inpainting)

≈

pa

tte
r

ns

Basic Model of Matrix Decomposition
nnegative matrix factorisation

F features

VK

N samples

H

≈

W

I

data V and factors W, H have nonnegative entries.

I

nonnegativity of W ensures interpretability of the dictionary, because
patterns wk and samples vn belong to the same space.

Introduction

Diﬃculties in NMF

Matrix
Decomposition
withConstraints
Constrained NMF methods
Diﬀerent types of constraints have been considered in previous works:
Sparsity constraints: either on W or H (e.g., Hoyer, 2004; Eggert and Korner,
2004);
Shape constraints on wk , e.g.:convex NMF: wk are convex combinations of inputs (Ding et al., 2010);
I harmonic NMF: wk are mixtures of harmonic spectra (Vincent et al., 2008).
I

Spatial coherence or temporal constraints on hk : activations are smooth(Virtanen, 2007; Jia and Qian, 2009; Essid and Fevotte, 2013);
Cross-modal correspondence constraints: factorisations of related
modalities are related, e.g., temporal activations are correlated (Seichepine
et al., 2013; Liu et al., 2013; Yilmaz et al., 2011);
Geometric constraints: e.g., select particular cones Cw (Klingenberg et al.,2009; Essid, 2012).
Essid & Ozerov (TPT/Technicolor)

A tutorial on NMF

ICME 201439 / 170

Matrix and Matrix Decomposition

BLIND SOURCE SEPARATION AND LINEAR GENERALIZED COM

NMF(non-negativity
constraints)• ICA (Independent Component
Analysis)• SCA (Sparse Component Analysis)• MCA (Morphological Component

ICA(independency
constraints)

SCA(sparsity constraints)

Analysis)• NMF (Non-negative Factorization)

MCA(morphologicalfeatures)Fig. 1.2 Four basic component analysis methods: Independent Compone
negative Matrix Factorization (NMF), Sparse Component Analysis (SComponent Analysis (MCA).

49 images
among2429from
fromMIT’s
MIT’s CBCL
49 images
among
2429
CBCLface
facedataset
datasetIntroduction

Motivation

Principal
Analysis
Explaining face
images Component
by PCA
MethodFacialfeatures

Vectorised images

Importance of featuresin each image

⇡
...

...

...

...
V

W

Essid & Ozerov (TPT/Technicolor)

A tutorial on NMF

HICME 2014

9 / 170

Principal Components
5 2nd Principal

1st Principal
Component, y1

Component, y2

Objective Function:<latexit sha1_base64="E74h6onRQXyZ8feFloPZszoa3Nk=">AAACKnicbVDLTsJAFJ3iC/FVdemmkZjghrTGRJeoG5eY8EookukwhQnTR2ZuRdL0e9z4K25YaIhbP8Qp1IjgSWZy5px7M/ceJ+RMgmlOtdza+sbmVn67sLO7t3+gHx41ZBAJQusk4IFoOVhSznxaBwactkJBsedw2nSGd6nffKJCssCvwTikHQ/3feYygkFJXf3G9vBzN1Y3DBw3HiWJzakLpV/hsfbDW62FxyixBesP4LyrF82yOYOxSqyMFFGGalef2L2ARB71gXAsZdsyQ+jEWAAjnCYFO5I0xGSI+7StqI89KjvxbNXEOFNKz3ADoY4Pxkxd7IixJ+XYc1RlOqhc9lLxP68dgXvdiZkfRkB9Mv/IjbgBgZHmZvSYoAT4WBFMBFOzGmSABSag0i2oEKzllVdJ46JsmWXr4bJYuc3iyKMTdIpKyEJXqILuURXVEUEv6A29ow/tVZtoU+1zXprTsp5j9Afa1zc/Q6mR</latexit>

4• PCA is to look for a low dimensional
projection in which the majority of signal
energy is kept.

3

2
4.0

4.5

5.0

5.5

6.0

• Here “Principal” represents “Major” thatthe projected signal has the largest
energy along the first principal direction(red line in the figure).

Introduction

Motivation

Principal Principal
Component Component
Analysis (PCA) Analysis
Recalling the technique1

Assuming the data is real-valued (vn 2 RF ) and centered (E[v] = 0),• PCA returns a dictionary WPCA 2 RF ⇥K such that the least squares

error is minimized:

1 X
WPCA = min
kvn
W N
n
• A solution is given by:v̂n k22 =1kVN

WWT Vk2F

WPCA = E1:K

where E1:K denotes the K dominant eigenvectors of Cv :X
1T
Cv = E[vv ] ⇡vn vn .
N n
1slide adapted from (Févotte, 2012).

Essid & Ozerov (TPT/Technicolor)

A tutorial on NMF

ICME 2014

7 / 170

PCA dictionary with K=25

PCA dictionary with K = 25

red pixels indicate negative values
14

pa

tte
r

ns

Nonnegative
Matrix Decomposition
Nonnegative
matrix factorisation

F features

VK

N samples

H

≈

W

I

data V and factors W, H have nonnegative entries.

I

nonnegativity of W ensures interpretability of the dictionary, because
patterns wk and samples vn belong to the same space.

Inonnegativity of H tends to produce part-based representations, because
subtractive combinations are forbidden.

Early work by Paatero and Tapper (1994), landmark Nature paper by Lee and Seung (1999)
12

NMF as a constrained minimization problem

NMF as a constrained minimisation problem

Minimise a measure of fit between V and WH, subject to nonnegativity:X
min D(V|WH) =d([V]fn |[WH]fn ),
W,H 0

fn

where d(x|y ) is a scalar cost function, e.g.,Isquared Euclidean distance (Paatero and Tapper, 1994; Lee and Seung, 2001)
I Kullback-Leibler divergence (Lee and Seung, 1999; Finesso and Spreij, 2006)
I Itakura-Saito divergence (Févotte, Bertin, and Durrieu, 2009)I ↵-divergence (Cichocki et al., 2008)
I-divergence (Cichocki et al., 2006; Févotte and Idier, 2011)
I Bregman divergences (Dhillon and Sra, 2005)
I and more in (Yang and Oja, 2011)
Regularisation terms often added to D(V|WH) for sparsity, smoothness,
dynamics, etc.
20

Common
NMFdesignalgorithm design
Common
NMF algorithm
I

Block-coordinate update of H given W(i 1) and W given H(i) .

I

Updates of W and H equivalent by transposition:
V ⇡ WH , VT ⇡ HT WTI

Objective function separable in the columns of H or the rows of W:
X
D(V|WH) =D(vn |Whn )nI

Essentially left with nonnegative linear regression:def

min C (h) = D(v|Wh)h 0

Numerous references in the image restoration literature. e.g., (Richardson, 1972;
Lucy, 1974; Daube-Witherspoon and Muehllehner, 1986; De Pierro, 1993)

21

NMF
with K=25
NMF dictionary
with dictionary
K = 25experiment reproduced from (Lee and Seung, 1999)
15

Multi-stream and cross-modal NMF schemes

NTF models

What is tensor?
Multi-way data representations
Some data can have more meaningful representation using multi-way
arrays rather than matrices (two-way arrays).
Electroencephalography (EEG) data (Lee et al., 2007)
channelVfrequency

frequency

channeltime

Essid & Ozerov (TPT/Technicolor)time

A tutorial on NMF

ICME 2014

107 / 170

What is tensor?
TENSOR PROPERTIES AND BASIS OF TENSOR ALGEBRA

29

Y
Mode-3Mode-1y 651

Mode-2
Fig. 1.13

A three-way array (third-order tensor) Y ∈ R7×5×8 with elements yitq .

nerally, tensors are denoted by an underlined capital boldface letters, e.g., Y ∈ RI1 ×I2 ×···×IN .
ntrast, matrices are denoted by boldface capital letters, e.g., Y; vectors are denoted by

Figure 1.3: A block matrix and its representation as a 4th-order tensor, createdby reshaping (or a projection) of blocks in the rows into lateral slices of 3rd-ordertensors.

What is tensor?

VectorMatrix

3rd-order Tensor

4th-order Tensor
...

...

...

...

...

...

...

...

...

...

...

A sample set

One sample

Scalar...

One-way
Univariate

2-way

3-way

4-way

5-way

Multiway Analysis (High-order tensors)
Multivariate

Figure 1.4: Graphical representation of multiway array (tensor) data of increasing
structural complexity and “Volume” (see (Olivieri, 2008) for more detail).

Tensor Fibers
TENSOR PROPERTIES AND BASIS OF TENSOR ALGEBRA

31

y

y

y13:

Fig. 1.15 Fibers: for a third-order tensor Y = [yitq ] ∈ RI×T ×Q (all fibers are treated as column
vectors).

Tensor slices

y

y13:

Fig. 1.15 Fibers: for a third-order tensor Y = [yitq ] ∈ RI×T ×Q (all fibers are treated as column
vectors).

Fig. 1.16 Slices for a third-order tensor Y = [yitq ] ∈ RI×T ×Q .

Q

Q

y: t q

q

...

...

q

Q

...

y1 t :

Y

...

...

...

q

Y

Y
Y: t :

Tensor Unfolding
TENSOR PROPERTIES AND BASIS OF TENSOR ALGEBRAA

I $ I2 I3

I3

I2

I1

33

Mode-1 unfolding: A(1) !"# 1
I2
I2
I2

I1
I2

I3
I $ I1 I3

I3
I1

Mode-2 unfolding: A(2) !"# 2
I1
I1
I1

I2

I1
I2

I3
I $ I1 I2

I3

I1

I1

Mode-3 unfolding: A(3) !"# 3
I1
I1
I1

I3
I2

I2

Fig. 1.19 Unfolding (matricizing) of a third-order tensor. The tensor can be unfolded in three
ways to obtain matrices comprising its mode-1, mode-2 and mode-3 vectors.

An Example of Tensor Unfolding
34

A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI

A::2

A

A::1
a111

a121

a131

a141

a211

a221

a231

a241

a311

a321

a331

a341

=

A(1) =

a111 a121 a131 a141
a211 a221 a231 a241
a311 a321 a331 a341

A(2) =

a111
a121
a131
a141

A(3) =

a111 a211 a311
a112 a212 a312

Fig. 1.20

a211
a221
a231
a241

a311
a321
a331
a341

a112

a122

a132

a142

a212

a222

a232

a242

a312

a322

a332

a342

a112 a122 a132 a142
a212 a222 a232 a242
a312 a322 a332 a342
a112
a122
a132
a142

a212
a222
a232
a242

a121 a221 a321
a122 a222 a322

a312
a322
a332
a342
a131 a231 a331a132 a232 a332

a141 a241 a341
a142 a242 a342

Example of unfolding the third-order tensor in mode-1, mode-2 and mode-3.

1.4.5.1 Outer Product The outer product of the tensors Y ∈ RI1 ×I2 ×···×IN and X ∈ R J1 ×J2 ×···×JM

The outer product of the tensors Y ∈ RI1 ×I2 ×···×IN and X ∈ RJ1 ×J2 ×···×JM is given by

1.4.5 Outer, Kronecker, Khatri-Rao and Hadamard Products
Z = Y ◦ X ∈ RI1 ×I2 ×···×IN ×J1 ×J2 ×···×JM ,Matrix Products(1.75)

zi1 ,i2 ,...,iN ,j1 ,j2 ,...,jM = yi1 ,i2 ,...,iN xj1 ,j2 ,...,jM .

(1.76)

Several special matrix products are important for representation of tensor factorizations and decompositions.
whereMatrix1.4.5.1
Outer
OuterProduct:
Product

Observe that, the tensor Z contains all the possible combinations of pair-wise products between the elements
The ofouter
productof the tensors Y ∈ RI1 ×I2 ×···×IN and X ∈ RJ1 ×J2 ×···×JM is given by
Y andX.
As special cases, the outer product of two vectors a ∈ RI and b ∈ RJ yields a rank-one matrix

Z = Y ◦ X ∈ RI1 ×I2 ×···×IN ×J1 ×J2 ×···×JM ,A = a ◦ b = ab ∈ R
T

where

I×J

(1.75)(1.77)and the outer product of three vectors: a ∈ RI , b ∈ RJ and c ∈ RQ yields a third-order rank-one tensor:

yI×J×Qzi1 ,i2 ,...,iZN =,j1 ,ji1 ,i2 ,...,ia2◦,...,jb ◦Mc =
∈R
, N xj1 ,j2 ,...,jM .

(1.78)(1.76)

wherethat, the tensor Z contains all the possible combinations of pair-wise products between the elements
Observeof Y and X.
(1.79)
z = ai bj cq .
As special cases, the outer product of twoijqvectorsa ∈ RI and b ∈ RJ yields a rank-one matrix

Matrix Kronecker Product:

A = a ◦ b = abT ∈ RI×J

1.4.5.2 Kronecker Product(1.77)Q
The outer
Kronecker
productof twovectors:
matrices aA∈∈ R
RII×J
B J∈and
RT ×Rc isdenotedas A ⊗ B ∈RIT ×JR tensor:, band∈R∈ aRmatrix
yieldsa third-orderrank-oneand theproductof threeand defined as (see the MATLAB function kron):

where

⎡

I×J×Q
,
Z=a◦b◦c ∈⎤R⎢ a21 B a22 B · · · a2J B ⎥⎢⎥⎢A⊗B = ⎢ .
.. ⎥
..
..
⎥
..
.
. ⎦
.
⎣
z =a b c .
ijq

=

'

1.4.5.2 Kronecker Product(1.78)

a11 B a12 B · · · a1J Bi

aI1 B aI2 B · · · aIJ B

j(1.80)(1.79)q(

a1 ⊗ b1 a1 ⊗ b2 a1 ⊗ b3 · · · aJ ⊗ bR−1 aJ ⊗ bR .

(1.81)The Kronecker product of two matrices A ∈ RI×J and B ∈ RT ×R is a matrix denoted as A ⊗ B ∈ RIT ×JR
and defined as (see the MATLAB function kron):⎡⎤It should be mentioned that, in general, the outer product of vectors yields a tensor whereas the Kronecker
J
T
Q
, b∈R
, c ∈a R
theirthree-way
product
gives
a vector. For
example,for thea ∈ofRvectors
It should
be mentionedthat,
in general,
thethreeoutervectors
productyieldstensorwhereasthe outer
Kronecker
I×T ×QQ
product= a ◦a bvector.
◦c ∈R
is a third-ordertensorwith thewhile
the three-way
three-way outer∈=RTa,j bct c∈q , R
theirproductYgives
For example,
for the threevectorsa ∈entries
RJ , byitq
×Q
Kronecker
a vector vec(Y)c ⊗ bthe⊗ aentries∈ RITQy. = a b c , while the three-way
product Yproduct= a ◦ bof◦ thec ∈sameRI×Tvectors
is aisthird-ordertensor=with

Matrix Products

itq
ITQ

Kronecker product of the same vectors is a vector vec(Y) = c ⊗ b ⊗ a ∈ R

j t q

.

Matrix
Hadamard
Product:
1.4.5.3
Hadamard Product1.4.5.3
Hadamard
The
Hadamard
product Productof two equal-size matrices is the element-wise product denoted by ! (or .∗ for
MATLAB notation) and defined as

The Hadamard product of two equal-size matrices is the element-wise product denoted by ! (or .∗ for⎡⎤
MATLAB notation) and defined as
a11 b11 a12 b12 · · · a1J b1J⎢ a⎡21 b21 a22 b22 · · · a2J b2J ⎥ ⎤⎢ a11 b11 a12 b12 · · · a1J ⎥b1J⎥.
A! B = ⎢..
..
..
⎢ ⎢ a.. b⎥⎣ ⎢ . 21 21 .a22 b22 . · · · . a2J ⎦b2J ⎥⎥A ! B = a⎢
⎢I1 bI1.. aI2 bI2.. · · · . .a.IJ bIJ ..
.
.
⎣ .
aI1 bI1

aI2 bI2

1.4.5.4 Khatri-Rao Product···aIJ bIJ⎥.
⎥⎦(1.88)(1.88)For two matrices A = [a1 , a2 , . . . , aJ ] ∈ RI×J and B = [b1 , b2 , . . . , bJ ] ∈ RT ×J with the same number of
1.4.5.4 J,Khatri-Rao
Product
columns
their Khatri-Rao
product, denoted by ⊙, performs the following operation:

Matrix Khatri-Rao Product:

For two matrices A = [a1 , a2 , . . . , aJ ] ∈ RI×J and B = [b1 , b2 , . . . , bJ ] ∈ RT ×J with the same number of(1.89)A ⊙ B = [a1 ⊗ b1 a2 ⊗ b2 · · · aJ ⊗ bJ ]columns J, their Khatri-Rao product,
denoted by ⊙, performs the following operation:'= vec(b1 aT1 )A ⊙ B = [a ⊗ b

The Khatri-Rao product is:• associative

The Khatri-Rao product is:associative•• distributive•• non-commutative
distributive

=

1

'

1

(

vec(b2 aT2 ) · · · vec(bJ aTJ )a ⊗ b ··· a ⊗ b ]vec(b1 aT1 )

2

2

J

J

∈ RIT ×J .

(vec(b2 aT2 ) · · · vec(bJ aTJ )A ⊙ (B ⊙ C) = (A ⊙ B) ⊙ C,

(A +AB)
A=
⊙C
⊙⊙
C,C,
⊙⊙
(BC⊙=C)(A+⊙BB)∈ RIT ×J .

(1.90)(1.89)(1.90)(1.91)(1.92)(1.91)

Tensor Matrix Product1.4.6 Mode-n Multiplication of Tensor by Matrix and Tensor by Vector, Contracted
Tensor ProductTo multiply a tensor by a matrix, we need to specify which mode of the tensor is multiplied by the columns(or rows) of a matrix (see Figure 1.21 and Table 1.1).
Definition 1.5 (mode-n tensor matrix product) The mode-n product Y = G ×n A of a tensor G ∈RJ1 ×J2 ×···×JN and a matrix A ∈ RIn ×Jn is a tensor Y ∈ RJ1 ×···×Jn−1 ×In ×Jn+1 ×···×JN , with elements
yj1 ,j2 ,...,jn−1 ,in ,jn+1 ,...,jN =Jn!

(1.97)

gj1 ,j2 ,...,JN ain ,jn .

jn =1
TENSOR PROPERTIESAND BASIS OF TENSOR ALGEBRA

37

The tensor-matrix product can be applied successively along several modes, and it is commutative, that is(G ×n A) ×m B = (G ×m B) ×n A = G ×n A ×m B,

(a)

(m =
/ n).

(1.98)The repeated (iterated) mode-n tensor-matrix product for matrices A and B of appropriate dimensions canGY1Abe simplified as(G ×n A) ×n B = G ×n (BA).

(1.99)For G ∈ RJ1 ×J2 ×···×JN and a set of matrices A(n) ∈ RIn ×Jn , their multiplication in all possible modes (n =1, 2, . . . , N) is denoted as×1

G × {A} = G ×1 A

(1)

×2 A

(2)

=

· · · ×N A(N) ,(1.100)and the resulting tensor has dimension I1 × I2 × · · · × IN . Multiplication of a tensor with all but one modeis denoted as
(1)
(n−1)
(n+1)
(N)
G
×
·
·
·
×
×
·
·
·
×
−n {A} = G ×1 A
n−1 A
n+1 A
N A
(4 ! 7)(7 ! 5 ! 8)

(1.101)(4 ! 5 ! 8)

giving a tensor of dimension I1 × · · · × In−1 × Jn × In+1 × · · · × IN . The above notation is adopted from[85].
It is not difficult to verify that these operations satisfy the following properties

(b)

"

#

(n)

" (N)

G × {A} (n) = A G(n) A

G

⊗A

B

(N−1)

(n+1)

··· ⊗ A

⊗A

(n−1)

#

(1) T

··· ⊗ A.

Y

(1.102)(4 ! 7)

Tensor Matrix Product(7 ! 5 ! 8)

(4 ! 5 ! 8)

(b)
G

Y2

B

×2

=

(9 ! 5)(7 ! 5 ! 8)

(7 ! 9! 8)

(c)
C

×3

(6 ! 8)Y3=(7 ! 5 ! 8)

(7 ! 5 ! 6)

Fig. 1.21 Illustration of the mode-n multiplications of a third-order tensor by matrices. (a) mode-Tensor Vector Contracted Product38

A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARIGca×¯ 2 b

×¯ 1

×¯ 3
=

(7 ! 5 ! 8)

(7 ! 1! 1)(1 ! 5 ! 8)

(5 ! 1 ! 1)

(1 ! 1 ! 8)

(8 ! 1 ! 1)

(1 ! 1 ! 1)

Fig. 1.22 Illustration of mode-n multiplication of a third-order tensor G by vectors, yielding scalary = G ×¯ 1 a ×¯ 2 b ×¯ 3 c. Note that the dimension of the result is reduced by one. For example,
multiplying a three-way (a third-order) tensor by a vector in mode-1 results in a 2-way tensor (amatrix).

2

Special Form of Tensors
42

A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI

Rank-one
tensor:A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI

c

c

b

Y
a

Y

a Y = a ◦ b ◦ c ∈ RI×T ×Q , a ∈ RI , b ∈ RT , c ∈ RQ .
Rank-one third-order tensor:

Fig. 1.23(I × T × Q)

Fig. 1.23b(I × T × Q)

I

Rank-one third-order tensor: Y = a ◦ b ◦ c ∈ RI×T ×Q , a ∈ RI , b ∈ RT , c ∈ RQ .

Examples of- -tensors
withspecial
forms
1 1

I

-1- 1-

1- 1- 1- 11- 1- 1- 1(a)- 1- 1-- 1- 1-

1

(b)

(c)

1

Fig. 1.24 Special forms of third-order tensors: (a) Super-Identity cube tensor I, (b) sparse tensor
with diagonal frontal slices, which can be mathematically expressed as I ×3 C, and (c) block
diagonal(b)(c)(a) tensor.

ig. 1.24 SpecialDECOMPposition
tensorforms of third-ordertensors: (a) Super-Identity
tensorFACtor
I, (b) sparse(CANDECOMP)
or PARAFAC cube
(PARAllel
decomposition)whichwith diagonal frontal
slices,can under
be mathematically
expressed
preserves
the whichuniquenesssome mild conditions[90].as I ×3 C, and (c) blockiagonal tensor.

The CP Format
The3.2.
CP Format

319

CP Approximation 317

xtended to multilinear models using tensor decompositions, such
...
b
b
1
R
e Canonical Polyadic (CP) and the Tucker
decompositions,
as il-Λ
...
= A
+ +
X
X
ated in Figures 3.1, 3.2 anda 3.3.
aR

1

R( I R)(I J )The CP Format

ar

B

br

( R R)

cR

c1

T

C

CP decomposition (also called the CANDECOMP, PARAFAC, or
b1 + ... +=R
nical Polyadic decomposition) decomposes
an Nbth-ordertensor, ...
X I
1 A
2 Λ
I
I
1
2
N
R
, into a linearbr
a1 combination ofa R terms, brarN
br , (which
1927;
I J K )are rank-1 tensors, and is given by (Hitchcock,( R R R)
( I R)
hman, 1970; Carroll and Chang, 1970)(R J )

cr( K R)RT

B

br(R J )

Figure 3.2: Analogy
between a low-rank matrix factorization, X
A BT
R
R
1
2
abrsimple
low-rank
Xr br (top),⁄r brandbrN tensor factorization (CP decomposir 1 ⁄r aRtion), X
br cr (bottom).
2 B
3 C
r1 A
1
r 1 ⁄r ar<latexit sha1_base64="xFykfbaMagIisjPc8WLIcSUeUlE=">AAACTXicbZFbS8MwFMfTedmct6mPvhSHMBFGK4K+CHMi+OCDgrvAOkeapltY2pTkVBilX9AXwTe/hS8+KCJmswVvBwL//E7Oyck/bsSZAst6Mgpz8wuLxdJSeXlldW29srHZViKWhLaI4EJ2XawoZyFtAQNOu5GkOHA57bjjs2m+c0elYiK8gUlE+wEehsxnBINGg4rnBBhGrp9000FSs/fSkxycprlyLnU/D6cOpz7UcnqWOsITkG+bqSPZcAR7tzf7OTvPeg4qVatuzcL8K+xMVFEWV4PKo+MJEgc0BMKxUj3biqCfYAmMcJqWnVjRCJMxHtKeliEOqOonMzdSc1cTz/SF1CsEc0a/VyQ4UGoSuPrkdE71OzeF/+V6MfjH/YSFUQw0JF8X+TE3QZhTa02PSUqAT7TARDI9q0lGWGIC+gPK2gT795P/ivZB3bbq9vVhtdHM7CihbbSDashGR6iBLtAVaiGC7tEzekVvxoPxYrwbH19HC0ZWs4V+RKH4CbIAteg=</latexit>

1B

1

2B

2

N B

(3.4)

N
<latexit sha1_base64="ZKZuGe5+IlpkdskkjM29VOeJ8XI=">AAACTXicbZHPS8MwFMfT+WNz/pp69FIcwkQY7RD0IqhD8OBBYZuDdY40TbewtCnJqzBK/0Evgjf/Cy8eFBGz2YK/HgS++by8l5dv3IgzBZb1ZBTm5hcWi6Wl8vLK6tp6ZWOzo0QsCW0TwYXsulhRzkLaBgacdiNJceByeuOOm9P8zR2ViomwBZOI9gM8DJnPCAaNBhXPCTCMXD/ppoOk1thLj3NwlubKudT9PJw6nPpQy2kzdYQnIN+epo5kwxHs3bb2c3ae9RxUqlbdmoX5V9iZqKIsrgaVR8cTJA5oCIRjpXq2FUE/wRIY4TQtO7GiESZjPKQ9LUMcUNVPZm6k5q4mnukLqVcI5ox+r0hwoNQkcPXJ6Zzqd24K/8v1YvCP+gkLoxhoSL4u8mNugjCn1poek5QAn2iBiWR6VpOMsMQE9AeUtQn27yf/FZ1G3bbq9vVB9eQss6OEttEOqiEbHaITdIGuUBsRdI+e0St6Mx6MF+Pd+Pg6WjCymi30IwrFT7WPteo=</latexit>

1
2
N
J
;
B
,
B
,
.
.
.
,
B
decomposition is rarely exact andK,has to be estimated by minimizing<latexit sha1_base64="fZieMeZZz1Yrk+n2aAQD91Uf968=">AAACTXicbZHPS8MwFMfT+Xv+mnr0UhzCRBitCnoRdEPw4GGCc4N1jjRNt2DalORVGKX/oBfBm/+FFw+KiNnWgjofBL75vLyXl2/ciDMFlvViFGZm5+YXFpeKyyura+uljc1bJWJJaJMILmTbxYpyFtImMOC0HUmKA5fTlntfH+VbD1QqJsIbGEa0G+B+yHxGMGjUK3lOgGHg+kk77SWVw730NAf1NFfOle7n4dTh1IdKTmupIzwB+fY8dSTrD2Dv7mY/ZxdZz16pbFWtcZjTws5EGWXR6JWeHU+QOKAhEI6V6thWBN0ES2CE07ToxIpGmNzjPu1oGeKAqm4ydiM1dzXxTF9IvUIwx/RnRYIDpYaBq0+O5lR/cyP4X64Tg3/STVgYxUBDMrnIj7kJwhxZa3pMUgJ8qAUmkulZTTLAEhPQH1DUJth/nzwtbg+qtlW1r4/KZ7XMjkW0jXZQBdnoGJ2hS9RATUTQI3pF7+jDeDLejE/ja3K0YGQ1W+hXFBa+AbkPtew=</latexit>

cost function.
functions
are typically of the Leaste ⁄ra suitable
are non-zero
entries of Suchthe cost
diagonalcore tensornnnRR and B(LS)nR are factor
Squares
type,norm
b1 ,inb2the, . .form. , bR of the
RInFrobenius
ma-The PARAFAC16 can be formulated as follows (see Figures 1.26 andGiven a data tensor Y ∈ RI×T ×Q and the positive index J, find three-com
matrices or factors, A = [a1 , a2 , . . . , aJ ] ∈ RI×J , B = [b1 , b2 , . . . , bJ ]TENSOR DECOMPOSITIONS
47 factorization:
RQ×JAND
whichFACTORIZATIONS
perform the following approximate

Alternative Representations of CP Decomposition(a)Y=(Q × J)

C

=

X=B

A
(I × J)

yitq =

E+(J × T)yitq(I × T × Q)

j=1

aij btj cqj + eitq .

(I × T × Q)

16 Also called the CANDECOMP (Canonical Decomposition) or simply CP dec

cJ

c1
+

+
bJ

b1

!
a1

aJ

(I × T ×Q)

b1 c1TbJ cJT
+

+

=
a1

J
"

eitq

(b)

Y

j=1

aj ◦ bj ◦ cj + E = !A, B, C" +or equivalently in the element-wise form (see Table 1.2 for various rep

T

Y

J
"

aJ

(c)
T

b1 1

J J

+

+

Alternative Representationsof CP Decomposition
=
a1

aJ

(c)

Yr = Y(1)
Y1

Y2

...

T

A

Xr = Br
X1

!

YQ

X2

...

XQ

(J × TQ)

Xq ! Dq X(I × J)(I × TQ)(q = 1, 2, . . . , Q)
(d)

DQ

!

A

Y1
(I × T × Q)

D1

X=B

(J × J)(J × T)T(I × J)

Yq = ADqX,(q = 1, 2, . . . , Q)

Fig. 1.27 The alternative representations of the third-order PARAFAC model: (a) as a set of

Alternative Representations of CP Decomposition

48

A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI(Q × J)

C c

j

Y

!

A

X=B

G

bj

aj
(I × T × Q)

(I × J)T(J × T )(J × J × J)

c1

cJ

λ1

+

λJ

+

b1

=
a1

bJ
aJ

Fig. 1.28 Harshman’s PARAFAC model with a superdiagonal core tensor G = Λ =!
diag(λ1 , λ2 , . . . , λ J ) ∈ RJ×J×J for the third-order tensor Y ! Λ ×1 A ×2 B ×3 C = Jj=1 λ j a j ◦ b j ◦ c j .
(In this model all vectors are normalized to unit length).

320

CP
Approximation
Multiway Component Analysis and Tensor Decompositions

Algorithm 1: Basic ALS for the CP decomposition of a
3rd-order tensor

Input: Data tensor X RI J K and rank R
Output: Factor matrices A RI R , B RJ R , C RK R , and scaling
vector ⁄ RR
1: Initialize A, B, C
2: while not converged or iteration limit is not reached do
3:A
X 1 C B CT C BT B
4:
Normalize column vectors of A to unit length (by computing thenorm of each column vector and dividing each element of avector by its norm)
5:B
X 2 C A CT C AT A
6:
Normalize column vectors of B to unit length
7:C
X 3 B A BT B CT C
8:Normalize column vectors of C to unit length,
store the norms in vector ⁄9: end while
10: return A, B, C and ⁄.

multiplying a matricized tensor and Khatri–Rao product (of factor

e original Tucker model makes the assumption of orthogonality of the factor matrice
y to SVD), [79], [17], [84], [83], [117], [106]. We will, however, ignore these constr
TuckerElementwise,
the Tucker
decompositionin (4.1)
isposingnonnegativity
constraints
theApproximation
problemof estimating
the component matrices a60
A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI
nsor is Pconverted
intoageneralized
NMF
problemcalledtheNonnegative
TuckerDe
Q
R!
!
!
(Q × Pimplementations
)
on
(NTD)
(see
Chapter7
fordetails).
The
first
of
Tuckerdecompos
gpqr aip bjq ckr for i = 1, . . . , I, j = 1, . . . , J, k = 1, . . . , K.
xijk ≈
C
onnegativity
constraints together with a numberof other constraints were given by K
p=1 q=1 r=1e and Bro in [79], [17]. The NTD imposes nonnegativity constraints for all compoTe Pand, Q,a andare thenumberof components(i.e., tocolumns)in theimposes
factor matrices
X=B
escoreRtensor,
while
a=semi-NTD(inanalogysemi-NMF)
nonnega+G
E
AY
B, andrespectively.
If P, Q,matrices
R are smallerI, elements
J, K, the ofcore
G can be
aints
to C,
only
some components
and/orthan
somethetensor
core tensor.
ught
as a compressedversionof X. Indescriptions
some cases, forthethestoragefor model
the decomposed
ere
areofseveral
equivalentmathematical
Tucker(see
Table(I × T × Q)
(I × T × Q)
(R × T )(I × J) (J × R × P)
sion
of
the
tensorcanbesignificantly
smaller
than
for
the
original
tensor; see Baderbe expressed in a compact matrix form using mode-n multiplications

1.34 Tucker
Tucker3 model
is a weighted sum of the
product of three(factors)
Kolda [17]. Fig.
The
decompositionis outerillustrated
in vectors
Figure4.1.stored
as columns of component matrices A ∈ R , B = X ∈ R
and C ∈ R . The core tensor
Most fittingGalgorithms
inCsection
assumefactormatrices∈ RY =
definesa (discussed
linking
structure
betweenthesetof4.2)
components
andC⟧
J, thatR,+andPthedenote
the
G
×
A
×
B
×
+
E
=
⟦G;
A,
B,
E,
(1
1 In order
2 to achieve3
numberofcomponents.
uniquenessfortheTuckermodelsitisnecessary
tocolumnwise impose
orthonormal,
but
this
is
not
required.
In
fact,
CP
can
be
viewedasaadditional constraints such as sparsity and nonnegativity.
cial case of Tucker where the core tensor is superdiagonal and P = Q = R.
"
= ⟦G;A,FormB, C⟧forms
is theshorthandnotationfor theY
The
matricized(one
perDecomposition:mode)of (4.1)
areTucker3 tensor decomposition.
Matrixof
Tucker
I×J

T

T ×R

Q×P

J×R×P

ng the unfolding
approach
we can form
obtain matrix forms expressed compactly by theor equivalently
in the element-wise
T
≈
AG
(C
⊗
B)
,
X
products:(1) J R P(1)!!!
yitq =
g jrp ai j btr cqp +T
eitq ,≈j=1BG
(C
⊗
A)
,
X(2)
(2)
r=1 p=1

I

T
T
Y
#
A
G
(C
⊗
B)
,
(1)
(1)
(B
⊗
A)
.
TX(3) ≈ CG
Q
(3)

(1.158)(1

where a j ∈ R , b j ∈ R , and c j ∈ R , (that is, the vectors within the associated component(factor) matrices A, B and C,), and g jrp are scaling factors whichT are the entries of a core tensor(2)(2)= [g jrp ] ∈ fromR J×R×P . the formulas
ese equationsGfollow
in sections 2.4 and 2.6; see [134] for further

Y

# BG

(C ⊗ A) ,(1

C
Different Types of Tucker Decomposition(Q x P)

62

A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI

Tucker 3
yitq
Tucker 3=
J !
R !
P
!
3
g jr p ai j btrTucker
cqp
j=1 r=1 p=1

C (Q x P) T
X=B
G C (Q×P) T
X=B
GT(J x R x P )X(R=x TB)GYY(I x T x Q)
Y

! A
! A
x J)
! (IA

(I x T x Q)

(I x J)(J x R x P )(R x T )(I ×T ×Q)

(I ×J)(J × R × P )(R × T )GG(JG
x R x Q)

X = BT
X = BTT) T
X (R=x B(a) Tucker3

Tucker 2
btr
yitq =
g jrq ai jTucker2
j=1 r=1
Tucker 2
J !
R
!

Y
Y
(I x T x Q)
Y

! A
! A
x J)
! (I A(I ×J)(J × R × Q )(R × T )(I x T x Q)

(I x J)(J x R x Q )(R x T )Y
Y(I ×T ×Q)
Y

! A
! A
! (I×AJ)(I ×T ×Q)(b) Tucker2

yitq =

J!
j=1

g jtq ai j

Tucker 1
Tucker 1
Tucker 1(I (c)x T x Tucker1
Q)(I x J)

G
G(J × T ×
GQ)(J x T x Q)

y3
(I J)
(I threeT Q) related Tucker(J T Q)
Fig. 1.35 Summary of thedecompositions.
y3
x

x

x

x

x

334
MultiwaySVD
Component
and Tensor Decompositions
From
MatrixtoAnalysis
Higher-order
Case
(a)

Eigenvector of XX
ur
RT

Rank of XX

sr

U
R

S

...

X

T

Eigenvector of X X

R

...

I

~
=

T

0

vr

T

Singular
value

V

0

0

R

J
(I J )(I J )(I I )(J J )( I3 R3)

(b)

I3

U

(3)

R3
I3
I1

( I1

(c)

I1
X

I2
I 2 I 3)

I1 U

( I1

(1)

R3
R1 St
R2

R1
R1)

I2
S

( I1 I 2 I3)

I4

(4)

I2 U(2)

( I2

R2
R2 )

336

HOSVD

Multiway Component Analysis and Tensor Decompositions

Algorithm 2: Sequentially Truncated HOSVD (Vannieuwenhoven et al., 2012)
IN
Input: N th-order tensor X RI1 I2
and approximation
accuracy Á
Output: HOSVD in the Tucker format X̂ JS; U 1 , . . . , U N K,
such that X X̂ F Á
1: S
X
2: for n 1 to N do
3:
U n , S, V
truncated_svd S n , ÁN4:SVS
5: end for
6: S
reshape S, R1 , . . . , RN
7: return Core tensor S and orthogonal factor matrices
U n RIn Rn .

Algorithm 3: Randomized SVD (rSVD) for large-scaleand low-rank matrices with single sketch (Halko et al.,2011)
Input: A matrix X

RI J , desired or estimated rank R, and6: S
reshape S, R1 , . . . , RN
7: return Core tensor S and orthogonal factor matrices

U

n

R

In Rn

.

rSVD

Algorithm 3: Randomized SVD (rSVD) for large-scaleand low-rank matrices with single sketch (Halko et al.,2011)
Input: A matrix X RI J , desired or estimated rank R, andoversampling parameter P or overestimated rank R R P ,exponent of the power method q (q 0 or q 1)
Output: An approximate rank-R SVD, X USVT , i.e., orthogonal
matrices U RI R , V RJ R and diagonal matrix S RR R withsingular values
1: Draw a random Gaussian matrix
RJ R ,
2: Form the sample matrix Y
XXT q X
RI R
3: Compute a QR decomposition Y QR
4: Form the matrix A QT X RR J
5: Compute the SVD of the small matrix A as A USVT
6: Form the matrix U QU.

Low multilinear rank approximation is always well-posed, however,HOOI
3.4. Higher Order SVD (HOSVD) for Large-Scale Problems

337

Algorithm 4: Higher Order Orthogonal Iteration (HOOI)(De Lathauwer et al., 2000b; Austin et al., 2015)
IN
Input: N th-order tensor X RI1 I2(usually in Tucker/HOSVD
format)Output: Improved Tucker approximation using ALS approach, withorthogonal factor matrices U n
1: Initialization via the standard HOSVD (see Algorithm 2)
2: repeat3:
for n 1 to N do
4:
Z
X p n Up T
5:C
Z n ZTn RR R
6:
Un
leading Rn eigenvectors of C
7:
end for
8:G
Z NUN T
9: until the cost function X 2F
G 2F ceases to decrease
10: return JG; U 1 , U 2 , . . . , U N K

where XBest is the best multilinear rank approximation of X, for a
specific tensor norm.

For the formulation of the NTF problem and the algorithms, the following symbols areused:outer product
Khatri-Rao product~Hadamard product46
A. CICHOCKI, R. ZDUNEK, A.H. PHAN, S. AMARI↵
element-wise division⇥nmode-n product of tensor
c1 and matrix
cJ
(N
)
(n+1)
(n
1)
(1)
A n =A
... A
A
... A

Nonnegative Tensor Factorization

b1

bJ

+

=
3.2 Problem
Y definition

+

+

E

The Nonnegative Tensor Factorization
problem can beaformulated
as nonnegative canona1
J(I × T × Q)
× Q)
ical decomposition/ parallel factor decomposition (CANDECOMP(I /× TPARAFAC)
asfollows (after [6]):Fig. 1.26 A graphical representation of the third-order PARAFAC as a sum of rank-one tensors.
All the(NTF).
vectors {aGiven
are N
treatedas column
vectorsof
2 ⇥...⇥IN and are linked for eachj , b j , c j } an
Definition-th order
tensor2 factor
RI1 ⇥Imatrices
and a positive integer!YJ
index j via the outer product operator, that is, Y = j=1 a j ◦ b j ◦ c j + E or equivalently in a compact
J, factorize
a2 setcomponent= Iinto×1 A ×B ×3of
C +NE. nonnegative(In this model notform YY
all vectors matrices
are normalized to unit length).
(n)

(n)

(n)

A(n) = [a1 , a2 , . . . , aJ ] 2 RIn ⇥J , (n = 1, 2, . . . , N ) representing the common (loading)factors, that is,factorization (PMF) method to the three-way
PARAFAC model, referring to the result as PMF3
J
X
(1)
(2)
(N )
(three-way positive matrix factorization). Although
based model
Y = Ŷ + E =aj
aj such. . constrained. aj + nonnegativityE=
may not match perfectly the input data (i.e., it may have larger residual errors E than the standard
j=1
PARAFAC without any constraints) such decompositions are often very meaningful and have(1)(2) [115], [116].
physicalinterpretation
[30],
I⇥
⇥2 A
· · · ⇥N A(N ) + E = JA(1) , A(2) , . . . , A(N ) K + E
1A
It is often convenient to assume that all vectors have unit length so that we can use the(n)Harshman’s PARAFAC modelwith ||amodified
1 givenand jby=[62],1, 2,[63]. . . , J.
j ||2 = 1 for n = 1, 2, . . . N

J"The tensor E is the approximation
error. Figure 2 illustrates the decomposition for aY=
λ j a j ◦ b j ◦ c j + E ! ⟦λ, A, B, C⟧,(1.124)third-order tensor.
j=1

Matrix Nonnegative Least-Squares (MNLS)
3.4. Nesterov-type algorithm for strongly convex MNLS

21

Algorithm 2: Nesterov-type algorithm for MNLS
Input: X 2 Rm⇥k , B 2 Rk⇥n , A0 2 Rm⇥n , tol > 0.
1 Compute W =XB, Z = BT B.
2 Compute L = max(eig(Z)) µ = min(eig(Z)).
3

Set Y0 = A0 ,ppLµ= pL+pµ , k = 0.

while (1) do
5
rf (Yk ) = W + Ak Z;
6
if (max(|rf (Yk ) ~ Yk |) < tol) then7
break;
8
else
⇥
⇤
1
9
Ak+1 = Yk L rf (Yk ) + ;10
Yk+1 = Ak+1 + (Ak+1 Ak );11
k = k + 1;
4

12

return Ak .

Terminating condition If ⇤ denotes the Lagrange multiplier matrix associated withthe matrix element-wise nonnegativity constraints in (3.8), then the Karush-Kuhn-Tucker(KKT) conditions for problem (3.8) are4.2. Nesterov-type algorithm for MNLS with proximal term

2

The Algorithm

Thus, fX can be expressed as

1
The objective function: fX (A, B, C) = 2 XA
1
=
XB
2
1
=
XC
2 NTF
4.3. Nesterov-based AO

A (C

B)

T

2
F

B (C

A)

T

2
F

C (B

A)T F .

(4.4

2

26

These expressions form the basis for the AO NTF in the sense that, if we fix two matri
Algorithm
4: Nesterov-based
AO NTFfactors,
then we can update
the third by solving an MNLS problem. For reasons relate
Input: with
X , Athe0, B0 0, Cof
0, MNLS, tol. problem, we propose to add a proximal term. Mor00 the
conditioning
1 Set k = 0kkkspecifically,
if
A
,
B
,
and
C
are the estimates of A, B, and C, respectively, after th
2 while (terminating condition is FALSE) do
k+1
iteration,3
WAk-th= AO
XA (CBk ) thenAA(Ck byBk )T (Ck Bk ) + I
kk , ZAis=given
4
Ak+1 = Nesterov MNLS(WA , ZA , Ak , , tol)
2
T
1
T
5
WB = XB (Ck AA
)
B
,
Z
=
(C
A
)
(C
+ I A k2 ,k+1 := argminkB
XAk A(Ck+1
Bk )k A+k+1 )kA
(4.5
k+1
k
k F
2
F
AB0
6
Bk+1 = Nesterov MNLS(WB , Z, B2k , , tol)
7
WC = XC (Ak+1 Bk+1 )Ck , ZC = (Ak+1 Bk+1 )T (Ak+1 Bk+1 ) + I
0 determinesto the proximal term. If (Ck Bk ) is a wel
8
Ck+1where= Nesterov
MNLS(WCthe, ZCweight, Ck , ,assigned
tol)
9
(Ak+1
, Bk+1 , Ck+1
) = Normalize(A, Bk+1 , Cto) small weight on the proximal term an
k+1putconditioned
matrix,then it is k+1
reasonable10(Ak+1 , Bk+1 , Ck+1 ) = Accelerate(Ak+1 , Ak , Bk+1 , Bk , Ck+1 , Ck , k)compute Ak+1 that leads to a large decrease of the function fX (A, Bk , Ck ). If, on the othe11k =k+1
hand, (Ck Bk ) is an ill-conditioned matrix, then it is reasonable to put large weight o
12 return Ak , Bk , Ck .
the proximal term, leading to a better conditioned problem and easy computation of Ak+
that improves the fit in fX (A, Bk , Ck ) but is not very far from Ak . This is the strateg