PROCEEDINGS OF THE IEEE

1

Linked Component Analysis
from Matrices to High Order Tensors:
Applications to Biomedical Data
Guoxu Zhou, Qibin Zhao, Yu Zhang, Tülay Adalı, Shengli Xie, and Andrzej Cichocki

arXiv:1508.07416v1 [cs.CE] 29 Aug 2015

F
Abstract—With the increasing availability of various sensor technologies, we now have access to large amounts of multi-block (also called
multi-set, multi-relational, or multi-view) data that need to be jointly analyzed to explore their latent connections. Various component analysis
methods have played an increasingly important role for the analysis
of such coupled data. In this paper, we first provide a brief review of
existing matrix-based (two-way) component analysis methods for the
joint analysis of such data with a focus on biomedical applications.
Then, we discuss their important extensions and generalization to
multi-block multiway (tensor) data. We show how constrained multiblock tensor decomposition methods are able to extract similar or statistically dependent common features that are shared by all blocks, by
incorporating the multiway nature of data. Special emphasis is given to
the flexible common and individual feature analysis of multi-block data
with the aim to simultaneously extract common and individual latent
components with desired properties and types of diversity. Illustrative
examples are given to demonstrate their effectiveness for biomedical
data analysis.

Index Terms—(Multiway) blind source separation (BSS), (multilinear)
independent component analysis, nonnegative/sparse matrix/tensor
factorizations, group and joint independent component analysis, independent vector analysis (IVA), CPD (PARAFAC) decompositions,
constrained Tucker decompositions for multi-block data, data fusion,
analysis of multi-relational data.

• Guoxu Zhou is with the Laboratory for Advanced Brain Signal Processing, RIKEN, Brain Science Institute, Wako-shi, Saitama,
Japan and the School of Automation at Guangdong University
of Technology, Guangzhou, China (Corresponding author, e-mail:
zhouguoxu@brain.riken.jp).
• Qibin Zhao is with the Laboratory for Advanced Brain Signal Processing,
RIKEN, Brain Science Institute, Wako-shi, Saitama 3510198, Japan (email: qbzhao@brain.riken.jp).
• Yu Zhang is with the Key Laboratory for Advanced Control and Optimization for Chemical Processes, Ministry of Education, East China
University of Science and Technology, Shanghai 200237, China (e-mail:
yuzhang@ecust.edu.cn)).
• Tülay Adalı is with the Department of Computer Science and Electrical
Engineering, University of Maryland, Baltimore County, MD 21250
USA (e-mail: adali@umbc.edu).
• Shengli Xie is with the School of Automation at Guangdong University
of Technology, Guangzhou, China (e-mail: shlxie@gdut.edu.cn).
• Andrzej Cichocki is with the Laboratory for Advanced Brain Signal
Processing, RIKEN, Brain Science Institute, Wako-shi, Saitama 3510198,
Japan, and with Systems Research Institute, Polish Academy of Science,
Warsaw, Poland (Corresponding author, e-mail: cia@brain.riken.jp).

1

I NTRODUCTION

Development and wide deployment of various sensor
technologies have resulted in the generation of large
amounts of multi-block data that need to be jointly
analyzed in order to extract physiologically meaningful and useful hidden (latent) components. Medical
recording and diagnostic devices play an important
role among those. In order to explore spatial, temporal, and spectral differences and similarities in multiarray, multi-block bio-signals and images, researchers
developed a variety of strategies based on component
analysis and multiway factor analysis. In this paper,
we focus on analysis of multi-block data that are
naturally connected using tensor/matrix factorizations.
We critically review emerging multiway component
analysis or multilinear blind source separation approaches that include procedures for extracting hidden
components (or source signals) with specific features
or constraints (e.g., sparsity, nonnegativity, statistical
independence). The presented approaches can be applied not only to unsupervised multilinear blind source
separation but also for missing data imputation (tensor
completion), feature extraction, classification, clustering and anomaly detection, especially in application
to electroencephalography (EEG), electrocorticography
(ECoG), and magnetic resonance imaging (MRI) data.
Multi-block tensor decompositions and multi-way analysis allow us to discover meaningful hidden spatiotemporal structures of complex data and perform generalizations by capturing multi-linear and multi-aspect
relationships. The challenge is how to analyze complex
inhomogeneous multi-set large-scale multi-array data.
We describe the key procedure to apply tensor/matrix
factorization/decompositions methods for these problems in practice and demonstrate a number of successful examples.
Tensors have been adopted in diverse branches of
data analysis, such as in signal and image processing, psychometrics, chemometrics, biometrics, quantum physics, quantum chemistry and neuroscience [1]–
[10]. Tensor decompositions provide extensions of blind
source separation (BSS) and 2-way (matrix) component
analysis (2-way CA) to multi-way component analysis
(MWCA) [3], [4], and have proven to be suitable for

2

...

X(i, :, :)

K
I

X(1) 2 RI⇥JK

X

...

A higher-order tensor can be interpreted as a multiway
array of numbers. Tensors are denoted by calligraphic
capital letters, e.g., X ∈ RI1 ×I2 ×···×IN (we shall assume
that all entries of a tensor are real-valued). The order of
a tensor is the number of its “modes”, “ways” or “dimensions”, which include e.g., space, time, frequency,
trials, subjects, classes, and dictionaries [1], [2]. A matrix
has two modes: rows and columns, while an N th-order
tensor has N modes.
Matrices (2nd-order tensors) are denoted by boldface
capital letters, e.g., X, and vectors (1st-order tensors)
by boldface lowercase letters; for instance the columns
of the matrix A = [a1 , a2 , . . . , aR ] ∈ RI×R are denoted
by ar and elements of a matrix (scalars) are denoted by
lowercase letters, e.g., ai,r (or [A]i,r ).
Subtensors are formed when a subset of indices is
fixed. Of particular interest are mode-n fibers (vectors),
defined by fixing every index but n, and slices which
are two-dimensional sections (matrices) of a tensor,
obtained by fixing all the indices but two, see Fig.1.
The process of unfolding flattens a tensor into a
matrix [2]. In the simplest scenario, mode-n unfolding
(matricization, flattening) of tensor X ∈ RI1 ×I2 ×···×IN

I

...

c Guoxu Zhou

...

1.1 Preliminaries – Basic Tensor Notations and Operations

mode-1 fibers X(:, j, k)
slices X(:, :, k)
...

dimensionality reduction of multi-way data, even if the
data are incomplete and corrupted by noise. Tensor
decompositions are particularly attractive for biomedical data analysis that exhibits not only large volume
but also very high variety. Thus, they are suited to
problems in bio- and neuro-informatics or computational neuroscience where data are collected in various
forms of large, sparse tabular, graphs or networks with
multiple aspects and high dimensionality. For historical
overview of tensor factorizations/decompositions and
important scientific contributions in these areas in some
chronological orders please check [2], [3] and references
therein.
Moreover, multi-block tensors, which arise in numerous important applications that require the analysis of
diverse and partially related data, can be decomposed
to common (or correlated) and individual (or uncorrelated) components. The effective analysis of coupled
tensors requires the development of new models and
associated algorithms that can identify the core relations that may exist among the different tensors, and
scale to large multi-sets data.
This paper extends beyond the standard tensor decomposition models such as the unconstrained single
block Tucker and canonical polyadic (CP) models, and
aims to demonstrate flexibilities of matrix/tensor decompositions for multi-block and multi-modal data,
together with their role as a mathematical backbone for
the discovery of hidden structures in large-scale multirelational data [1], [3].

PROCEEDINGS OF THE IEEE

J

...

...

J

X(2) 2 RJ⇥IK

X(:, j, :)

...

K

...
unfolding

X(3) 2 RK⇥IJ

Fig. 1: Fibers, slices, and unfolding of a 3rd-order tensor.
yields a matrix X(n) ∈ RIn ×(I1 ···In−1 In+1 ···IN ) whose
columns are all the mode-n fibers arranged in a specific
order (in this paper rows and columns are ordered
colexicographically). See Fig.1 for an illustration.
The mode-n product of the tensor X ∈ RI1 ×···×IN
and a matrix B ∈ RJ×In is the tensor C =
×J×In+1 ×···×IN
X ×n B ∈ RI1 ×···×In−1
with entries
PIn
ci1 ,i2 ,...,in−1 ,j,in+1 ,...,iN = in =1 xi1 ,i2 ,...,iN bj,in . This can
be also expressed in a matrix form as C(n) = BX(n) .
The Kronecker product ⊗ and the KhatriRao product
of matrices areN defined in
(p)
TABLE 1. Moreover, we define
=
p6=n A
(N )
(n+1)
(n−1)
(1)
A
⊗J
···⊗A
⊗A
· · · ⊗ A , and similarly
define p6=n A(p) .
See Table 1 for the other useful notations. We refer
readers to [2], [3] for more details about various tensor
operations.

2

T WO - WAY B LIND S IGNAL S EPARATION

Blind source separation (BSS) is an unsupervised learning method with the aim of estimating R source signals
B ∈ RT ×R from the measurement matrix X ∈ RI×T ,
such that [11], [12]
X = AB> =

R
X
r=1

ar bTr =

R
X
r=1

ar ◦ br ,

(1)

where the outer product ◦ is defined in TABLE 1,
A = [a1 , a2 , . . . , aR ] ∈ RI×R is the unknown mixing
matrix (also known as the basis matrix or dictionary,
depending on application) and B = [b1 , b2 , . . . , bR ] ∈
RT ×R is the matrix of sources br (components or latent variables). Although we explicitly assume that the
columns of B are the latent variables, please keep in
mind that the roles of A and B can be exchanged if we
consider X> in (1).
Without any constraint, the matrix factorization problem (1) is in general highly undetermined as it actually

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

3

TABLE 1: Notations And Definitions
B, br , bm,r
B> , trace(B)
krB
X, X(n)
kBkF
~,
C = A⊗B
C=A

B

C = X ×n B
X = b(1) ◦ b(2) ◦ · · · ◦ b(N )

Matrix, the rth-column, and the (m, r)th-entry of B. [B]m,r = bm,r .
Transpose, trace of matrix B
Kruskal rank of B is the maximum value ensuring that any subset of krB columns is linearly independent.
Tensor, the mode-n matricization of tensor X.
P
1
Frobenius norm of B ∈ RM ×R , i.e., kBkF = ( m,r b2mr ) 2 . For tensor X, kXkF = X(n) F , ∀n.
Element-wise (Hadamard) product, division of matrices or tensors.
Kronecker product of A ∈ RI1 ×J1 and B ∈ RI2 ×J2 yields C = [ai1 j1 B] ∈ RI1 I2 ×J1 J2 with entries
c(i1 −1)I2 +i2 ,(j1 −1)J2 +j2 = ai1 j1 bi2 j2 .



Khatri-Rao product of A = a1 a2 · · · aJ ∈ RI1 ×J and B = b1 b2 · · · bJ ∈ RI2 ×J
I
I
×J
1
2
yields a matrix C = a1 ⊗ b1 a2 ⊗ b2 · · · aJ ⊗ bJ ∈ R
.
Mode-n product of tensor X ∈ RI1 ×···×IN and matrix B ∈ RJ×In yields tensor C = X ×n B ∈
P In
RI1 ×···×In−1 ×J×In+1 ×···×IN with entries ci1 ,i2 ,...,in−1 ,j,in+1 ,...,iN =
in =1 xi1 ,i2 ,...,iN bj,in
and matrix representation C(n) = BX(n) .
Outer product of vectors b(n) ∈ RIn yields a rank-1 tensor X = b(1) ◦ b(2) ◦ · · · ◦ b(N ) with entries
(1) (2)
(N )
xi1 ,i2 ,...,iN = bi1 bi2 · · · bi .
N

has an infinite number of solutions. However, with very
little a priori information, BSS enables the recovery of
sources in the sense that
B̂ = Ψ(X) = BΛP,

(2)

where B̂ is an estimate of B, Ψ denotes a suitable BSS
algorithm, Λ is a diagonal scaling matrix and P is a permutation matrix, denoting the unavoidable scaling and
permutation ambiguities, respectively [12]. This feature
is called essential uniqueness of BSS. From a more
general sense, (1) essentially provides a representation
or transformation of the data matrix X that preserves
as much of the information as possible. Different types
of information to be preserved lead to various BSS
criteria and methodologies. Below we briefly introduce
the most popular ones.
Principal component analysis (PCA). PCA transforms
data matrix X into a set of uncorrelated variables, called
the principal components using either eigenvalue or
singular value decomposition (EVD/SVD). A common
use of PCA is for dimensionality reduction where only
a set of principal components is retained to preserve the
maximum variance in the data. In PCA, uniqueness is
achieved by imposing orthogonality on the transform
matrix.
Independent component analysis (ICA). If we rewrite the
matrix product in (1) using the random vector notation
x(t) = As(t),

t = 1, 2, . . . , T

(3)

thus interpreting each row of X as T -dimensional realization of the entries of the random vector x(t) and
the rows of B> in (1) as the T -dimensional realization
of the entries of the random vector s(t), then we can
pose the problem as one that estimates a demixing
matrix W such that the rows of Wx are as independent
as possible. In other words, the demixing matrix W
is expected to be the pseudo-inverse of the mixing
matrix A, neglecting the scaling and permutation ambiguities. With this principle, we need a criterion to
measure the statistical independence, and the most

natural criterion is the maximization of independence
through minimization of mutual information which can
be shown to include others such as maximum likelihood and maximization of non-Gaussianity [12]–[14].
Another approach is to explicitly estimate higher-order
statistics such as cumulants to establish independence
using joint diagonalizations as in joint approximate
diagonalization of eigenmatrices (JADE) [15]. ICA can
be viewed as a generalization of PCA in the sense that
it requires to venture beyond second-order statistics
to incorporate higher-order statistical information [15].
ICA is attractive because by only using the natural
assumption of statistical independence, it enables an
essentially unique decomposition without the need for
any additional constraint. For example, when both nonGaussianity and sample dependence are taken into
account in the formulation, any source, including those
that are Gaussians can be identified subject to permutation and scaling ambiguities as long as the Gaussian
sources do not have proportional covariances [13].
ICA has proven to be very useful for functional MRI
(fMRI) data analysis, and can be performed in two different ways [16]–[18]: namely spatial ICA that extracts
unique independent spatial maps, and temporal ICA
that extracts independent time courses by considering
X> in (1). Although they are essentially just two different ways for organizing the data, spatial ICA is more
useful as the spatial independence assumption is well
suited for the systematically non-overlapping nature
of the spatial patterns for most cognitive activation
paradigms [13], [19].
Sparse component analysis (SCA). Sparsity is a natural
property of many real-word signals, for example, activities of our brain, eye movements, and heart beat
among others. SCA is a technique that allows us to
recover latent source signals by exploiting the sparsity
of source signals and/or the mixing matrix; and it
is even able to separate dense signals, because some
smooth dense signals often can be very sparse in a
transform domain [20]. For example, EEG signals are

4

Trial 2

Trial 3
Trial 2
Trial 1

Trial 3

......

Channel

Channel

Trial 1

Channel

Channel

CWT

Trial

0

Time

Time(s)

0

TF
R

C

Frequency

Frequency

much more sparse in the frequency domain than in the
time domain [21]. In this case, BSS can be performed in
the transform domain using SCA and then the outputs
can be inversely transformed to construct the original
source signals. SCA is particularly important because,
by incorporating compressed sensing [22], it allows the
solution of underdetermined BSS problem, i.e., when
there are fewer observations than sources (I < R) in (1).
Before we can do this, of course, we need to know the
mixing matrix A, which is often estimated by clustering
the observations using the sparsity assumption on the
sources [23], or by blind identification methods utilizing
the statistical independence of sources [24], [25].
Nonnegative matrix factorization (NMF). In power spectral analysis for EEG data, the latent frequency components are expected to be nonnegative and generally
very sparse. In this scenario, NMF proves to be particularly useful as it imposes nonnegativity on both
the mixing matrix and the source signals [26]. Due
to the nonnegativity constraints, the data matrix is
represented as purely additive combinations of basic
building blocks while without any subtractions, endowing NMF with the special ability of learning-parts
[26]. While NMF achieved great success in machine
learning, lack of uniqueness guarantee in the general
case poses difficulty when it is applied to BSS. Fortunately, by imposing additional sparsity—even very
weak sparsity—on the sources, NMF is able to provide unique factorization [27], [28]. Typically, under the
pure-source dominance condition [29], [30], that is, for
each signal there exists at least one instant at which
only this signal is active or strongly dominant, NMF can
be a powerful tool to separate statistically dependent
nonnegative sources [29], [30]. NMF based on the puresource dominance assumption is called separable NMF,
which is highly scalable and can be applied to discover
the most representative members in a large database
[31].
If we further require each row of B to be pure-source
dominant, i.e., the whole matrix B to be a sub-matrix of
a permutation matrix, we obtain the orthogonal NMF
model. Orthogonal NMF is quite useful as it proves to
be equivalent to K-means clustering [32], [33], where
the columns of A are the cluster centroids and B is the
cluster indicator matrix. This representation is certainly
unique as it meets the strongest form of the puresource dominance condition. In summary, NMF and its
various variants have been powerful and versatile tools
in machine learning and signal processing.
Smooth Component Analysis (SmCA). Smoothness is
another important intrinsic property that provides advantages, especially for biomedical image and video
processing. In SmCA, we impose smoothness constraints not on the raw data X, but rather on the hidden
components, i.e., vectors of the factor matrix B, and/or
on basis vectors of the mixing matrix A. SmCA can
be implemented by solving the following optimization

PROCEEDINGS OF THE IEEE

Channel

0

Time

Trial 1
Trial

Trial 2
Trial 3

Fig. 2: Illustration of how to generate high-order tensor
data for EEG data analysis, by incorporating the timefrequency representation (TFR) of each channel.

problem [34], [35]:
min X − AB>
A,B

2
F

+ γ1 kL1 Ak2,1 + γ2 kL2 Bk2,1 ,

(4)

where k·kF is the Frobenius norm defined in TABLE 1,
PT
L2,1 matrix norm is defined as kYk2,1 = t=1 kyt k2 for
Y ∈ RI×T , γ1 , γ2 are nonnegative penalty parameters
that control level of smoothness, and matrices L1 , L2
represent difference operators, typically, the matrix with
all its entries being zero, except 1 on its diagonal and
−1 on its superdiagonal [1], [34], [35].

After we build a proper objective function, BSS is
finally achieved by solving an optimization problem.
While many standard optimization approaches can be
applied, some typical methods that are especially suitable for BSS are briefly introduced here. Taking into
account that the mixing matrix in ICA is full column
rank and hence has a certain underlying structure, the
concept of natural gradient was proposed and proved
to be the steepest direction [36], [37]. On the other hand,
for NMF and SmCA, alternating descent methods are
quite useful since more than one parameter matrix need
to be optimized. In these types of methods, each time
the optimization is performed over a subset of parameters (typically, one factor matrix or even only one column of it) with the others fixed [1], [38]. The alternating
direction method of multipliers (ADMM) [39], which
has been enjoying growing popularity recently and is
well suited to large-scale optimization problems, also
falls into this category. Basically, alternating descent
is rather a principle than a concrete method and has
been widely applied to solve large-scale optimization
problems that arise in machine learning and related
areas [39].

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

M ULTIWAY DATA A NALYSIS
B(2) ∈ RI2 ×R2
I3
R2
I1

=

X

R1 G

I2

R3

B(3) ∈ RI3 ×R3

B(1) ∈ RI1 ×R1

(a) Tucker model

cy

X
Time

Sp

l
tra (3)
ec
b1
(2)

b1

1

=

Sp

l
tra (3)
ec
bR
(2)

bR

R

+···+

Temporal

(1)

b1

Spatial

n
ue
eq
Fr

Spatial

So far the data we considered were represented as
a single matrix, which has only two dimensions, or
modes. In many applications the data can be naturally
represented as high dimensional arrays, i.e., tensors.
For example, in EEG, each trial can be formed by a
matrix with two modes of channels (electrodes) and
time, and multiple trials form a 3rd-order tensor of
channel×time×trial. If we consider the time-frequency
representation of each channel, we obtain a 4th-order
tensor of channel×time×frequency×trial, see Fig.2. To
analyze such type of data, we may simply unfold them
to form a large matrix and then apply the variety aforementioned BSS methods. However, this often causes
loss of internal structure of high-order data. It is more
desirable to take the tensor structure into account when
we analyze such data.

Channel

3

5

Temporal

(1)

bR

(b) CP model

3.1

Basic Tensor Decomposition Models

The Tucker Decomposition. Given an N th-order tensor
X ∈ RI1 ×I2 ···×IN , by Tucker decomposition, it can be
represented as [1]–[3]
(1)

(2)

(N )

X = G ×1 B ×2 B · · · ×N B
(5)
.
= JG; B(1) , B(2) , · · · , B(N ) K,
h
i
(n)
(n)
where B(n) = b(n)
∈ RIn ×Rn is the
b
·
·
·
b
1
2
Rn
mode-n factor component matrix consisting of latent
(n)
components brn as its columns, G ∈ RR1 ×R2 ···×RN
is the core tensor reflecting the connections (or links)
between the components, and (R1 , R2 , . . . , RN ) with
Rn = rank(X(n) ) is called the multilinear rank of X.
Fig.3(a) illustrates the Tucker decomposition of a 3rdorder tensor.
Unconstrained Tucker decomposition is generally not
unique, because an equivalent representation can be
obtained by replacing B(n) with B(n) Q and G with
G ×n Q−1 in (5), where Q is any Rn -by-Rn nonsingular
matrix. In other words, B(n) estimates the range of
X(n) . Hence, we often apply truncated SVD to moden matricization X(n) for n = 1, 2, . . . , N , in order to
obtain orthogonal factor matrices B(n) . This method is
well known as the high-order SVD (HOSVD) and has
been widely applied to multidimensional data analysis
[40]–[43].
The canonical polyadic (CP) decomposition (also known
as PARAllel FACtor analysis (PARAFAC) or CANonical
DECOMPosition (CANDECOMP)). If we restrict the
core tensor G in (5) to be diagonal with R1 = R2 = · · · =
RN = R, we obtain the CP model that decomposes
the target tensor into the sum of rank-one terms (see
Fig.3(b)):
XR
(2)
(N )
X=
λr b(1)
r ◦ br ◦ · · · ◦ br
r
(6)
.
= JΛ; B(1) , B(2) , . . . , B(N ) K,

Fig. 3: Illustration of the Tucker decomposition and the
CP decomposition for 3rd-order tensors.

where ◦ denotes the outer product of vectors, and Λ =
diag(λr ). The minimum number of R is called the rank
(or explicitly, CP-rank) of X [2].
Using the mode-n matricization operator, (5) and (6)
can be rewritten in matrix factorization forms as
(n) (n)>
X>
B
, n = 1, 2, . . . , N,
(7)
(n) = A
N

(p)
where A(n) =
G> for the Tucker dep6=n B
J (n)

(p)
composition, or A(n) =
Λ for the CP
p6=n B
decomposition. Both models are multilinear models,
and we often update the factor matrices in an alternating fashion by solving N least squares (LS) problems
2
minA(n) X(n) − B(n) A(n)> F , see [2], [3], [38], [44] for
details.
Tensor decompositions have been widely applied
to brain signal processing [21], [45]–[47]. In [48], the
CP decomposition was employed to extract timefrequency-space atoms to identify the different signal
components of the EEG data, which renders tensor
decomposition an efficient approach for processing EEG
data. In [49], [50], tensor-based EEG preprocessing was
introduced to perform localization of brain sources
based on EEG measurements, which can provide improved performance in the case of several simultaneously active brain regions and low signal-to-noise ratios. In [51], a tensor toolbox for event-related potentials
(ERP) analysis was developed, which includes multiway decomposition methods and visualization of timefrequency transformed ERP data. The group analysis
can be performed naturally by tensor representation of
brain activities recorded from varying subjects, conditions or repeats. Based on several strategies to represent
EEG in a tensor format, tensor decompositions, espe-

6

PROCEEDINGS OF THE IEEE

cially, nonnegative tensor decompositions, are applied
to extract the dominant sources of brain activity in
an unsupervised manner, which can also be used for
analysis of group differences between subjects.
As a natural extension of BSS, multiway BSS is an
attractive and promising approach because it allows
components to be simultaneously extracted from different domains (i.e., modes) of tensor data. For example,
in EEG data analysis, via multiway BSS we can extract
the spectral components, the temporal components, and
the spatial components simultaneously, together with
links between them, reflected by the core tensor. As
essential uniqueness is one key feature of BSS, next
we discuss unique tensor decomposition models and
associated algorithms.
3.2

Multiway Blind Source Separation

n=1

κn ≥ 2R + (N − 1),

r=1

where cr is the cumulant of sr (t), and ar is the rth
column of A. As a result, the mixing matrix A can
be estimated via the CP decomposition of Cx , and the
identifiability is guaranteed by the uniqueness of the
CP decomposition of Cx . This approach is attractive
because it is able to estimate the mixing matrix of even
underdetermined mixing systems [12], [15], [24], [25].
If the sources have certain temporal structures, e.g.,
each source signal sr is either a deterministic ergodic
sequence or a stationary process [56], [57], the covariance matrices of x(t) satisfy
Cx (τ ) = A diag(cs (τ )) A> ,

As the Tucker decomposition is not unique in general
sense, we often need to impose additional constraints
on the factor matrices and the core tensor to obtain a
unique and physically meaningful representation. Comparing (7) with (1), the columns of X(n) are just linear
mixtures of the columns of B(n) . This motivates us to
extract B(n) by applying BSS algorithm to X>
(n) directly,
which has been discussed in [52], [53] in terms of
multilinear ICA and multiway BSS. Multiway BSS can
be considered as an extension of the high-order SVD
where orthogonality constraints on factors are replaced
with the other ones such as statistical independence,
nonnegativity, sparsity, smoothness, depending on a
priori knowledge about the components [3], [40]. As
constrained Tucker decompositions, multiway BSS is of
particular interest in practice because it provides essentially unique components that have specific statistical or
deterministic properties (e.g., statistical independence,
smoothness, sparseness and nonnegativity [?]). Particularly, different from 2D BSS, multiway BSS is able to
accommodate all these diverse constraints or properties
of factors in one Tucker model, allowing desirable
interpretation of data and flexible feature extraction.
The CP decomposition can be interpreted as the
Tucker model with additional structure constraints, and
hence is essentially unique under the well-known sufficient condition [54]:
N
X

statistics) of x just forms a higher-order symmetric
tensor
R
X
Cx =
cr ar ◦ ar ◦ · · · ◦ ar ,
(9)

(8)

where κn is the Kruskal rank of B(n) (see TABLE 1), and
R is the rank of X = JΛ; B(1) , B(2) , . . . , B(N ) K. More relaxed conditions can be found in [55]. This uniqueness
analysis is purely based on multilinear algebra while
neglecting the physical attributes of factors.
One of the most important applications of the CP
decomposition in BSS is blind identification of the
mixing matrix. If we consider the ICA problem defined
in (3), the cumulant (can also be any other high-order

τ = 0, 1, 2, . . . ,

(10)

where τ is a time-lag, and diag(cs (τ )) is a diagonal matrix whose diagonal entries are csr (τ ) = E[s(t + τ )s(t)],
see [56] for details. It can be seen that (10) is just
the joint approximate diagonalization (JAD) of a set of
symmetric matrices Cx (τ ), τ = 1, 2, . . ., and has been
extensively studied within the last two decades, see
e.g., [56], [58]–[60]. If we stack Cx (τ ), we can obtain
a 3rd-order tensor C satisfying
C=

T
X
r=1

ar ◦ ar ◦ cs (τ ),

(11)

where s(τ ) is a vector of s(τ ), ∀ τ . (11) allows us to
estimate A by performing the CP decomposition of
tensor C. Again, the identifiability of A is guaranteed
by the uniqueness the CP decomposition of tensor C.
By comparing (7) with (1), we can also incorporate ICA to CP decompositions. A tensorial extension of probabilistic ICA (PICA), named tensor-PICA,
was proposed in [61] for multisubject fMRI analysis.
In tensor-PICA, assuming one component matrix, say
B(n) , has mutually statistically independent columns
(components), a two-dimensional ICA criterion is used
to update B(n) while the LS criterion is applied to
update the other factor matrices, and all factor matrices
are updated in an alternating manner. However, this
method suffers from a divergence issue caused by the
different objectives of the ICA step and the LS step
[62]. In [62] a new method combining ICA and the
CP decomposition was proposed, which overcomes the
disadvantages of tensor-PICA but still suffers from
high computational complexity. In [63] a more efficient
method was proposed, where one factor matrix B(n)
was extracted by applying BSS to X(n) (see (7)) at
first, then, by virtue of J
the special Khatri-Rao prod(p)
uct structure of A(n) =
, all the remaining
p6=n B
factor matrices can be recovered by using a series of
SVDs of rank-1 matrices. Particularly, once the matrix
extracted by using BSS is of full column rank (which is

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

almost always true), all the other factor matrices can
be uniquely identified, even if the general algebraic
uniqueness condition (8) for the CP decomposition is
not satisfied, which somehow extends the uniqueness
analysis of the CP decomposition. In [64] the robust
CP decomposition was also proposed to separate underdetermined mixtures of event-related sources that
include quasi-periodic signals (e.g., electrocardiogram
(ECG)), sources with synchronized trials (e.g., ERP),
and amplitude-variant sources.
Note that under the nonnegativity constraints on
the factors, J
the special Khatri-Rao product structure
(p)
of A(n) =
in (7) often makes A(n) very
p6=n B
(p)
sparse [38], even if B
(p 6= n) have weak sparsity.
The sparsity of B(n) allows us to uniquely estimate
B(n) and A(n) first by applying NMF to X(n) , and
then recover all factor matrices B(p) from A(n) , p 6= n,
which is an efficient way to perform nonnegative CP
decomposition (also known as nonnegative tensor factorization). Nonnegativity is important for the CP decomposition because the optimal low-rank approximation under nonnegativity constraint always exists [65]
and is almost always unique [66], which is not true for
unconstrained CP decomposition. Hence, nonnegative
CP decomposition is a very promising way to perform
nonnegative multiway data analysis [1], [38], [67].

4

L INKED A NALYSIS

OF

M ULTI - BLOCK DATA

Multi-block data typically include multiple measurements of the same phenomenon under various experimental conditions. For example, EEG signals or
fMRI data of the human brain in response to a certain
stimulus, but from different subjects and trials, can be
grouped together and naturally linked as multi-block
data. Multi-block data may also refer to multi-modal
data where different types of pulse sequences can be
used in a single scanning session to form multimodal
MRI datasets for each subject, or modalities such as
fMRI and EEG might need to be considered jointly [45].
All these possibilities pose the challenge to develop systematic approaches for fusing these various data types
together to automatically find common patterns—if
they exist—of related changes across multi-block data
[8], [13], [68]. Thus, linked component analysis is an
emerging research topic with great potential in biomedical data analysis as it allows the extraction of not
only shared common components across multi-block
data but also individual components specific for each
dataset.
A number of models have been introduced to enable
analysis of multi-block data by either defining hard
links among matrices/tensors through shared factors,
or by soft statistical links as in the case of independent
vector analysis (IVA) and canonical correlation analysis
(CCA). To define these models, first consider a set of
data matrices Xk ∈ RIk ×T , k = 1, 2, . . . , K. Examples

7

include brain imaging data where each Xk is represented as voxel×time points or voxel×subject. These
data blocks are naturally linked as they are collected
across a single shared mode (time T ), and can be
represented as in (1) for all K blocks (see Fig.4)
Xk = Ak B>
k , k = 1, 2, . . . , K.

(12)

In what follows, we first present models where the
main assumption is the existence of common components across all the datasets and is imposed through
the model directly as hard constraints, then in Section 4.2, we consider the use of soft links that take
advantage of statistical relationship, typically in the
form of statistical dependence as in IVA, CCA, and
multiset CCA (MCCA), and finally introduce models
that explicitly consider and model common as well as
distinct components within each dataset.
4.1

Models with common factors

To investigate shared common information across
multi-block data, e.g., common spatial maps, a number
of models have been proposed [13], [19], [69]:
Concatenation-based ICA models, as the name suggests, concatenate the data blocks either in the vertical dimension as in Group ICA [16] or horizontally
as in joint ICA [70]. In a simple concatenation, the
assumption is that all observations are spanned by a
set of exactly the same components. For example, as
illustrated in Fig.4, in
Xk = Ak B> , k = 1, 2, . . . , K,

(13)

B is shared by all data blocks whereas the mixing
matrix Ak ∈ RIk ×R is associated with Xk only. By
concatenating the observations vertically to form a tall
i>
h
P
>
>
∈ R( k Ik )×T , we
matrix X = X>
X2 · · · XK
1
have
e >
X = AB
(14)
h
i>
e = A> A> · · · A> . Note that if horiwhere A
1
2
K
zontal concatenation is applied (assuming that Ik = I,
∀k in Fig.4), we simply transpose X, which leads to
a model that assumes a common mixing matrix—the
matrix A in (3)—for the ICA decomposition. This is
the model adopted in joint ICA where we have Xk =
AB>
k , k = 1, 2, . . . , K, and the common (mixing) matrix
A is used to fuse multi-modal data such as fMRI and
EEG represented as subjects by voxels and subjects by
time points respectively [70]. In this application, the
columns of the mixing matrix include subject covariations and when there are multiple groups, such as patients and healthy controls, biomarkers can be identified
through a t-test on these columns. Or, as explained in
[69], multiple blocks Xk can be multi-task and multisubject data of the same type and dimension and can
be stacked either vertically or horizontally depending

8

PROCEEDINGS OF THE IEEE

Linked Component Analysis
(Common components dominant in one mode)
=

X1
I1 ⇥ T

A2R

G2

A2

PK

k=1

Ik ⇥J

A k Gk
.
G 2 RJ⇥R⇥K

Tucker decomposition

B>

Ak = A 2 RI⇥J ,
Ik = I, Jk = J, 8k.

I 2 ⇥ J2 J2 ⇥ R R ⇥ T

B>

B>

A

...

...

I2 ⇥ T
=

XK

Group BSS

B>

I 1 ⇥ J1 J1 ⇥ R R ⇥ T
=

X2

G1

A1

A

IK ⇥ T

AK

GK

CP decomposition
Gk is a square
diagonal matrix,
Ak = A 2 RI⇥R ,
Ik = I, Jk = R, 8k.

B>

IK ⇥ JK JK ⇥ R R ⇥ T

G 2 RR⇥R⇥K
A

B>

Fig. 4: Decomposition models interpreted as linked component analysis of multi-set data: Xk = Ak B>
k =
Ak Gk B> , where the latent components, i.e., the columns of B are assumed to be exactly the same for all data
blocks.
TABLE 2: Decomposition models for linked component analysis of multi-block data using common components
Decomposition

Linked component analysis
(k = 1, 2, . . . , K)

models
Concatenation-based models
joint ICA, group ICA, and variants

X>
k

=

BA>
k

or Xk = Ak

B> ,

where B contains the shared common components, see (13).
Xk = A Σk B> , where A and B contain the shared common components, see (15).

CP decomposition

Σk is a diagonal (scaling) matrix, A can be an arbitrary matrix, or has a Khatri-Rao product structure,
e.g., A = M1

Tucker decomposition

Xk = A Gk

M2 , which is equivalent to the CP decomposition of a higher-order tensor.

B> ,

A and B contain the shared common components, see (17).

A is an arbitrary matrix, or has the Kronecker product structure, e.g., A = M1 ⊗ M2 , which is
equivalent to the Tucker decomposition of a higher-order tensor.

on the interpretation. An important step included in
group ICA, which has been implemented in the GIFT
toolbox (available at http://icatb.sourceforge.net), is
double dimension reduction stages using PCA. If each
Xk represents a subject data such as fMRI (time points
by voxels), then we first perform a subject level PCA,
and after vertical concatenation of dimension-reduced
subject data, a second level PCA is applied at the group
level to estimate a common group subspace [16], [71].
Individual subject maps are then reconstructed using
the group and subject-level PCA matrices thus preserving most of the variability for individual subjects.
Other implementations and applications of the group
ICA model are also possible and discussed in [69].
Compared with (1), the concatenated ICA model described in (13) is merely the partitioned version of standard BSS (1), which is a popular technique to solve the
computational issue in the case where X is extremely

large in one dimension. In addition, this model allows
for explicit modeling of common modes across e.g.,
specific trials or subjects. One major limitation of the
joint ICA model is the strong assumption that all data
blocks share exactly the same common independent
components. The approach proposed in [72] uses CCA
and MCCA [73] such that the assumption of identical
mixing matrix (hence subject covariations of [70]) is
relaxed and instead maximally correlated subject covariations are estimated.

Tensor ICA stacks all the data blocks to create a
3rd-order tensor X, rather than 2D matrices obtained
by concatenation. This technique has been used for
multi-subject fMRI analysis based on the CP model,
for tensors with dimensionality of voxel×time×subject
[61], [74], [75]. In this method, data for each subject can

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

be expressed as
X:,:,k = Xk = A diag(wk ) B>
= [A diag(wk )] B> .

(15)

Note that the mixing matrix associated with Xk is
Adiag(wk ), implying all blocks actually share the same
mixing matrix but differently scaled by their individual
weighting vectors wk . This of course will lead to a more
restrictive model in comparison with (14), which can be
clearly observed from the relationship
X>
(2)

= (W

>

A)B ,
(16)
i
where W = w1 w2 · · · wR ∈ RK×R , X(2) is the
mode-2 unfolding of the data tensor X, (W A) serves
as the mixing matrix and has a Khatri-Rao product
structure, and A captures the latent spatial maps with
additional independence constraints. Note that X>
(2) is
just the vertical concatenation of Xk . In other words,
the major difference between the concatenation based
ICA (14) and the tensor ICA (16) is that: In the concatenation based ICA, the mixing matrix is not constrained
whereas in the tensor ICA it has a Khatri-Rao product
structure. The authors of [74] claimed that it can be a
beneficial feature of the tensor ICA because it avoids
unnecessary duplication of the spatial maps and can
allow them to be inferred more accurately when the
assumption holds.
It is worth noticing that (15) can be viewed as
asymmetric joint diagonalization of matrices Xk , an
extension of the symmetric JAD problem described in
(10), which shows that the CP decomposition in general
is merely a technique that simultaneously diagonalizes a
set of matrices (or more generally, tensors). Remembering the discussion in Section 2, SVD diagonalizes
a single matrix with orthogonality constraints whereas
the CP decomposition jointly diagonalizes multiple matrices but without additional constraints. Hence the CP
decomposition is essentially a useful tool for linked
component analysis of multi-block data that are of the
same size and share common factors in each dimension.
For example, in [76], the CP decomposition was used
for shared processing of perception and imagery of
music, in order to determine if both the spatial location
and temporal evolution are shared between tasks, identifying common networks and their timing signatures.
A natural question that may arise is whether it is
meaningful in practice to relax model (15) to
h

X:,:,k = Xk = AGk B> ,

(17)

where Gk , associated with Xk , is not necessarily diagonal any more. This interesting extension has been
studied in [77] in terms of population value decomposition (PVD) and has been applied to discover common features from thousands of similar images. This
simple model allows to reduce a very large number of
images to a manageable set of parameters and to make

9

statistical inferences on samples of large-scale datasets.
The authors demonstrated the potential of this method
by applying it to analyze the sleep heart health study
dataset. Moreover, the entries of matrix Gk can be used
as predictors in a regression context, e.g., to predict the
risk of Alzheimer’s disease using fMRI [77]. If we stack
Xk to form a 3rd-order tensor X, and similarly for Gk ,
by using the Tucker model, (17) can be rewritten using
tensor representation as (see Fig.4)
X = G ×1 A ×2 B,

(18)

implying the multi-block data Xk = X:,:,k share common factors A and B, but with subject-specific information kept in Gk = G:,:,k .
In summary, (15) and (17) are special cases of the CP
model and the Tucker model, respectively, which can be
extended to higher-order tensors [78] and applied to the
scenarios where multi-block data are of the same size
and share common components in all modes. Moreover,
the concatenated ICA models (group ICA and joint
ICA), the CP decomposition, and the Tucker decomposition all perform linked analysis of multi-block data
with the assumption that all samples are dominated
by the common components, and with the exception of
joint ICA that uses horizontal concatenation, all models
assume the blocks to be of the same size. We summarize
these models in TABLE 2.
4.2 Models with statistical links: IVA, CCA, MCCA,
PLS, and their tensorial extensions
The methods discussed in Section 4 assume that the
data blocks share exactly the same latent (independent) components, which could be too restrictive for
some applications and lead to suboptimal solutions.
In contrast, IVA [13], [79], CCA [80] and its multiset
extension MCCA [73], as well as partial least squares
(PLS) are more flexible methods, and provide advantages especially when not all modes are common across
the blocks. It has been shown that MCCA can be used
to perform BSS [81], fusion of EEG, fMRI and structural
MRI data [72], and for brain computer interface (BCI)
[82], and IVA for multi-subject fMRI analyis [13], artifact rejection in concurrent EEG-fMRI data [83] among
others. In what follows, we first briefly review matrix
decomposition using IVA, discuss how it generalizes
CCA and MCCA—like the generalization of PCA using
ICA—and then discuss the multilinear extensions for
CCA and PLS [84].
Independent Vector Analysis (IVA) can be formulated by
considering the model in (3) for K datasets such that
xk (t) = Ak sk (t),

k = 1, 2, . . . , K.

(19)

Hence, in this model, as opposed to those introduced
in Section 4.1, all datasets are assumed to have their
own unique mixing and source/component matrices.
The link across the datasets is achieved in a “soft

10

PROCEEDINGS OF THE IEEE

...
...
...

· · · ...
...

= AK

S2

.

···

A2

..

XK

=

S1

t
.
. . den
n
pe

···

X2

A1

de
dependent in

X1

...

SCV
=

SK

Fig. 5: Diagram of the IVA model.

manner” by defining an optimization criterion that
fully takes advantage of the statistical dependence that
exists among the multiple datasets. The key definition
in the formulation of IVA is the source component
vector (SCV) that is formed by using the corresponding
elements of the source random vectors sk (t) ∈ RI
such that an SCV is given by [s1,i (t), s2,i (t), . . . , sK,i (t)]>
for i = 1, 2, . . . , I, where the second subscript refers
to the source index for each dataset in (19). The IVA
decomposition is achieved by minimizing the mutual information among the I, each K-dimensional
SCVs, as opposed to sources in ICA [13], [85]. Hence,
IVA achieves a decomposition of independent sources
within each dataset while estimating SCVs that are
dependent within its entries (see Fig.5). It is this dependence that helps with the mitigation of permutation
ambiguity for sources (modes) across the datasets, and
hence enables IVA to identify sources that are independent and identically distributed Gaussians as long as
they are correlated across the datasets. Complete identification conditions and bounds are given in [79]. By
properly aligning corresponding independent components from each dataset—whenever they are related—
IVA enables a convenient and useful way to perform
multi-subject analysis as in fMRI [13] and allows one
to exploit dependence for other uses such as in artifact
rejection [83]. It has been shown that IVA is able to
preserve subject variability compared to group ICA
both by simulations [86], [87] and also using real data
where there is significant variations among the subject
fMRI scans, such as in patients with stroke [88] and
schizophrenia [89].
When the SCV distribution is assumed to come
from a multivariate Gaussian distribution, it can be
shown that IVA cost function becomes equivalent to
one of the five cost functions introduced for MCCA
[13], [73]. In [73], maximization of correlation across
elements of an SCV, and hence of multiple datasets,
is achieved by making the covariance matrix of an
SCV as ill-conditioned as possible through five costfunctions introduced by heuristic arguments, and each
SCV is thus estimated in a deflationary manner. All five
cost functions reduce to CCA when there are only two
datasets. MCCA has been applied to both multi-modal
data fusion and multi-subject fMRI analysis [72], [81]

and hence IVA can be considered as a generalization
of MCCA to the case where higher-than second-order
statistical relationship among the datasets are taken
into account by selecting an appropriate multivariate
density model for the SCV [13].
Multilinear extensions of CCA and PLS. Suppose we
have K pairs of sample instances {(Xk , Yk ), k =
1, 2, . . . , K}, with Xk ∈ RI1 ×I2 ×···×IN and Yk ∈
RJ1 ×J2 ×···×JM , which can be concatenated to form two
large tensors with additional sample dimension: X ∈
RI1 ×···×IN ×K and Y ∈ RJ1 ×···×JM ×K . We seek projection
vectors such that the correlation between the multilinear projections of X and Y on them is maximized:
max{wX ,wY }

√

u> v
u> uv> v

,

(20)

where
u = X(N +1) wX ,

v = Y(M +1) wY ,

(21)

wX and wY are just the vectorizations of rank-1 tensors
WX and WY , respectively, and
(1)

(2)

(N )

(1)

(2)

(M )

WX = wX ◦ wX ◦ · · · ◦ wX ,

WY = wY ◦ wY ◦ · · · ◦ wY .

(22)

The vectors u and v are called the first pair of canonical
variables, reflecting the common information shared
by the two groups X and Y. It turns out that (20)
can be easily solved via HOSVD, by applying the Lagrange multiplier method to (20). Then after deflation,
by solving (20) we seek the second pair of canonical
variables, and so on as in MCCA. Note that the canonical variables can cross multiple modes of X and Y,
as long as they are of the same size in these modes.
Multilinear CCA has been extended to the kernel space
and applied, for example, to behavior data classification
[84].
Basically, CCA enables extraction of components that
are highly correlated and is widely applied for feature
extraction. In contrast, PLS is a similar well-established
framework with the aim to predict a set of dependent variables (responses) from a set of independent
variables (predictors) through the extraction of a small
number of latent variables. It solves the following optimization problem:
max{wX ,wY }

u> v,

(23)

where u and v are defined as in (21). Compared with
(20), it can be seen that besides correlation, the energy
(variance) of u and v plays a key role in PLS. Again, (23)
can be solved by HOSVD. Once u> v is maximized, we
replace v in (21) by u that can be expressed by using
X and wX . As such, a connection between X and Y
is built, which allows us to predict Ynew from a given
new sample Xnew . In practice, we often extract multiple
pairs of vectors to improve the prediction accuracy,
following a similar routine as for the multilinear CCA.

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

Two standard tensor decomposition frameworks (i.e.,
CP and Tucker) are both considered as the counterpart
of SVD for tensorial data. The multilinear CCA and PLS
described in (20)–(23) are based on the CPD framework.
Therefore, an alternative strategy is to apply Tucker
decomposition in these methods, which results in a
more flexible model in terms of fitting ability. For
instance, multiway PLS based on Tucker model, called
HOPLS [90], can be described as
max{GX ,wX ,GY ,wY }

u> v,

(24)

where u and v are the same as in (21), but wX and wY
are the vectorizations of WX and WY , respectively, and
(1)

(2)

(N )

(1)

(2)

(M )

WX = GX ×1 WX ×2 WX × · · · ×N WX ,

WY = GY ×1 WY ×2 WY × · · · ×M WY .

(25)

Hence, (22) can be considered as a special case of (25);
and HOPLS is a generalized framework for multiway
PLS, which has shown the promising performance for
decoding of movement trajectories of a monkey’s hands
based on ECoG signals [3], [90]. Similar strategy can be
also extended to multilinear CCA.
4.3 Common and Individual Feature Analysis
(CIFA)
The models discussed in Section 4.1 assume that a
certain set of components are exactly the same across all
datasets, on the other hand, those discussed in Section
4.2 rely on statistical correlation (second and higherorder) to link the multiple datasets. Another approach
is to explicitly model the common and distinct components, which is expected to be a better match to
most physical systems, and hence lead to substantially
improved performance. For example, if we analyze the
human EEG signals in response to a certain stimulus,
but from different subjects and for many trials, it is
expected that these data would share some common
components due to the application of the same stimulus while also possessing components specific to each
individual. To address such scenarios, a more general
linked matrix factorization model was recently studied
in terms of a joint and individual variation explained
(JIVE) [91] and common and individual feature analysis
(CIFA) models [91]–[93]. In these methods, a set of
matrices X = {Xk ∈ RIk ×T : k = 1, 2, · · · , K}, share
at least one common dimension T . They can be, e.g., a
set of multichannel EEG signals (channel×time point)
recorded for the same visual stimulus but from different
subjects. These data matrices can be factorized in a
linked way as follows (see Fig.6(a))
"
#
B̄>
>
Xk = Ak Bk = [Āk Ăk ]
B̆>
k
(26)
=Āk B̄> + Ăk B̆>
k
.
=X̄k + X̆k , ∀k,

11

where B̄ ∈ RT ×C , B̆k ∈ RT ×R̃k , Rk = R̃k + C is the
number of latent components with Rk < Ik and C is
the number of common components, Āk and Ăk are the
partitions of the mixing coefficients Ak corresponding
to B̄ and B̆k . In this way, each data matrix Xk is
represented through a combination of components from
the common (shared) subspace X̄k = Āk B̄> and the
individual (intrinsic) subspace X̆k = Ăk B̆>
k . Basically,
the common information carried by the common subspaces X̄k can be used to identify the whole group
of coupled data, whereas the individual subspaces X̆k
can be used to identify the data member Xk as it
is presented only in Xk . The objective is to estimate
the matrix B̄ representing common components shared
by all data blocks and the matrices B̆k representing
individual components associated with only the data
matrix Xk .
The objective of both JIVE [91] and CIFA [92], [93]
is to solve this challenging problem we describe above,
with the subtle yet important difference: JIVE quantifies the amount of joint variation between data types
whereas CIFA captures highly correlated components
while ignoring their energy (variances), which allows
relatively weak common components to be detected
and extracted. For this reason, JIVE can be viewed as
the extension of PLS while CIFA is more closely related
to CCA. CIFA advances this topic by introducing BSS
methods to the analysis of each common/individual
subspace, motivated by the link between (1), (13), i.e.,
by defining
X̄k = Āk B̄> ,

X̆k = Ăk B̆>
k.

(27)

This formulation allows us to separate common and
individual components that are physically easy to interpret. Hence CIFA enjoys significantly improved flexibility and has demonstrated numerous potential applications in clustering analysis and pattern recognition
[93]. This paradigm is illustrated in Fig.7.
If each data block is a high-order tensor rather than
a matrix, more sophisticated models are required to
perform comprehensive linked analysis. For simplicity,
at first we focus on the common and individual feature
analysis with respect to only one mode, and then we
briefly discuss its extension to multiple modes. To avoid
cumbersome notations, sub- and super- scripts, we also
assume that the tensors are of order 3, but please keep
in mind that they can be of any order and may also
have different sizes—as long as they share at least one
common mode, from which common and individual
components will be extracted. Consider a set of tensors
X={X1 , X2 , . . . , XK } that share a common mode, say,
mode-1. Suppose Xk ∈ RT ×Ik,2 ×Ik,3 , where Xk can also
be a matrix with Ik,3 = 1, ∀k. Using the Tucker model,
each tensor can be decomposed as
(1)

(2)

(3)

Xk = Gk ×1 Bk ×2 Bk ×3 Bk , ∀k,

(28)

12

PROCEEDINGS OF THE IEEE

COBE

CIFA: beyond CCA

Common spaces

COBE

X̄k = Xk B̄B̄> ,
k = 1, 2, . . . , K.

B̄

Xk ,
k = 1, 2, . . . , K.

CFE

F̄
BSS, ICA,
NMF, SCA,
SmCA, ...

X̆k ⇡ Ăk B̆>
k,
k = 1, 2, . . . , K.

X̆k = Xk X̄k ,
k = 1, 2, . . . , K.

IFE

F̆k ,
k = 1, 2, . . . , K.

Individual subspaces

Dimensionality
reduction

Fig. 7: The paradigm of common and individual feature analysis (CIFA) for multi-block matrix data analysis. In
CIFA, the common orthogonal basis extraction (COBE) approach [93] can be employed to separate the common
and individual subspaces. Then, BSS can be respectively applied to the common and individual subspaces, to
achieve common feature (i.e., F̄) extraction and individual feature (i.e., F̆k , k = 1, 2, . . . , K) extraction.
where Gk ∈ RRk,1 ×Rk,2 ×Rk,3 is the core tensor of Xk ,
(n)
and Bk , n = 1, 2, 3, are the corresponding component
matrices. Assuming that the mode-1 factor matrices
(1)
(1)
Bk share some common components such that Bk =
(1)
[B̄(1) B̆k ], we obtain the following linked Tucker
decomposition model:
(2)

(3)

Xk =Ḡk ×1 B̄(1) ×2 Bk ×3 Bk
(1)

(2)

(3)

(29)

+ Ğk ×1 B̆k ×2 Bk ×3 Bk

=X̄k + X̆k ,

∀k,

(1)

where B̄(1) is shared by all data blocks whereas B̆k
is associated with Xk only, Ḡk ∈ RC×Rk,2 ×Rk,3 , Ğk ∈
R(Rk,1 −C)×Rk,2 ×Rk,3 , and C is the number of common
components in mode-1.
For notational simplicity, we temporally denote the
mode-1 matricization of X by X(k,1) , and similarly for
G(k,1) , Ḡ(k,1) and Ğ(k,1) ; then (29) can be rewritten in
parallel with (26) as
(1)

(1) >

X>
(k,1) = Ak Bk

(1)

(1)

(1) >

= Āk B̄(1)> + Ăk B̆k

where k = 1, 2, . . . , K, and


(1)
(3)
(2)
Āk = Bk ⊗ Bk Ḡ>
(k,1) ,


(1)
(3)
(2)
Ăk = Bk ⊗ Bk Ğ>
(k,1) .

,

model allows us to simultaneously rotate the columns
(1)
of Bk and conversely the rows of G(k,1) to achieve the
separation of common and individual subspaces. Afterwards BSS can be applied to respective subspaces to
extract unique and physically meaningful common and
individual components with specific properties, such as
sparsity, nonnegativity, statistical independence.
More generally, multi-block high-order tensor data
may share more than one mode, each of which contains
some common components. To model such a complex
dataset, we can combine these modes to form a big
mode that may contain some common components, by
reshaping each data tensor Xk simultaneously using
the technique introduced in [78]. As such, the matrix
B̄(1) in (30) is replaced by the Kronecker product of
common components coming from different modes.
In other words, we can impose additional structure
(1)
constraints upon B̄(1) and B̆k . This problem can be
solved using the ALS [78]; implementation details are
however omitted due to space limitation.

(30)

5

D IMENSIONALITY R EDUCTION FOR I NCOM PLETE DATA C ORRUPTED BY O UTLIERS
(31)

The key difference between (30) and (26) lies in the
(1)
(1)
multilinear structures of Āk and Ăk shown in (31).
We call (29) or (30) CIFA-Tucker as the constrained
Tucker decomposition model is a key model of this
linked component analysis. To solve this problem, we
may apply the ALS method used in JIVE, but with
additional multilinear structure constraints imposed on
(1)
(1)
Āk and Ăk . Alternatively, we can apply the two-step
method proposed in [93]: in the first step dimensionality reduction can be performed to obtain (29). Then,
the rotational ambiguity associated with the Tucker

When we are ready to apply any specific method we
have discussed to real-world problems, several practical
issues arise. These include questions such as to how
many latent components one should extract, and how
to remove the noise and outliers from noisy and/or
incomplete observation data. Particularly, multi-block
high-order data may cause data deluge and pose new
challenges for effective and efficient data analysis. A
common approach to deal with these issues is dimensionality reduction, which has been a ubiquitous
technique broadly applied for noise reduction, data
compression, and feature extraction. Indeed, most models and algorithms we discussed above rely on efficient
and robust dimensionality reduction.

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

13

TABLE 3: Comparison of popular dimensionality reduction models (c.f. [97]). Y ∈ RI×T is the data matrix of
T samples in I-dimensional space, X = AB> is the extracted low-rank data, the sparse matrix E denotes the
outliers. Laplacian matrix Φ characterizes a simple graph or a hypergraph between the samples of Y [97], L2,1
PT
norm is defined as kYk2,1 = t=1 kyt k2 . Due to space limitation, most of their nonnegative variations are omitted
here as they are straightforward, and their multiway extension can be found in [98]–[100] and references therein.
Model

Objective

PCA

minA, B

Constraints
Y − AB>

minX, E kXk∗ + λ kEk1

RPCA [95]
RPCA on graphs [97]
Graph Laplacian PCA (GLPCA) [101]
Robust GLPCA [101]
NMF [26]
Orthogonal NMF [32], [33]
Sparse component analysis
Smooth component analysis (SmCA) [34]

X,E

X

Y = X + E, X = AB>
trace(XΦX> )

2
minA, B Y − AB> F + γ trace(B> ΦB)
minA, B Y − AB> 2,1 + γ trace(B> ΦB)
2
minA, B Y − AB> F
2
minA, B Y − AB> F
2
minA, B Y − AB> F + λ kBk1
2
minA,B X − AB> F + γ1 kL1 Ak2,1 +
γ2 kL2 Bk2,1

s.t. Y = X + E.

(32)

However, (32) is a highly nonconvex optimization problem that is difficult to solve; a tractable surrogate of
(32) is obtained by replacing rank(X) by the nuclear
norm kXk∗ = σi (X) and kEk0 by kEk1 . By solving this
surrogate problem, both low-rank approximation and
the outliers can be estimated. Detailed discussions on
this topic can be found in [95] and references therein.
A closely related problem, known as matrix completion (MC), also utilizes the low-rank feature of data and
can be modeled as
min rank(X),

A> A = I

minX, E kXk∗ + λ kEk1 + γ

The key to achieve dimensionality reduction is the
low-rank nature of data that is essentially related to
the complexity of the generative model. Estimating
this rank, or equivalently, the number of latent components, is often a challenging task and referred to
as the model selection problem. As an example, in
brain data imaging we never know how many latent
sources that contribute essentially to the measured data;
different model order selections may lead to quite
different latent components and, consequently, different
interpretations.
A successful model order selection largely relies on
appropriate assumption on the distribution of noises.
Gaussian noise is usually assumed in many imaging
systems [94], which often results in a least squares
problem and the use of PCA introduced in Section 2. To
deal with grossly corrupted observations, robust PCA
(RPCA) was recently proposed and attracted considerable attention. In RPCA, the noisy data Y is modeled
as the sum of a low rank term X = AB> and a sparse
term E of outliers, leading to the following optimization
problem [95]:
min rank(X) + γ kEk0 ,

2
F

s.t. XΩ = YΩ ,

(33)

Y = X + E, X = AB>
B> B = I
B> B = I
A and B are nonnegative
A and B are nonnegative, and B> B = I
Various constraints, e.g. nonnegativity
Various constraints, e.g., nonnegativity.
L1 , L2 represent difference operators.

where Ω indicates location of available entries of Y, XΩ
is a matrix with [XΩ ]i,j = xi,j if (i, j) ∈ Ω; otherwise
0. Apparently, if we know the location of outliers in
(32), problem (32) is equivalent to (33) because we can
simply mark the items corrupted by outliers as missing
values. We refer readers to [96] for a comprehensive
overview, and here, focus on tensor completion. Tensor
completion has many important applications: 1) it allows us to perform effective data analysis even if the
data is incomplete or corrupted by outliers; however,
these incomplete data are often simply discarded using traditional techniques; 2) it allows us to predict
the missing items based on partially observed values,
which is known as collaborative filtering and has been
widely applied in recommender systems; 3) it promises
exact recovery of signals from very limited samples.
This feature makes the completion technique somewhat related to compressed sensing [22]. However,
completion mainly makes use of the low-rank property
whereas the compressed sensing employs usually the
sparsity of latent signals under a specific dictionary.
Extension of (32) and (33) to the tensor domain is
rather straightforward if the Tucker model is used to
model the data. For example, rank(X) can be simply
replaced by the weighted sum of multilinear ranks of
tensor X, and then outlier detection and completion
can be achieved by minimizing the multilinear rank.
Efficient variations and algorithms can be found in
[102]–[104]; a recent trend is to handle both outliers
and missing values using a unified framework [103].
Among these methods, the infTucker method [104]
assumes that the multilinear rank is given and the other
methods can learn the multilinear rank but introduce
new tuning parameters whose selection is often as
difficult as the model order selection itself.
The rank minimization framework faces big chal-

14

PROCEEDINGS OF THE IEEE

X1

Ā1

...
=

B̄>

ĂK

C⇥T

+

IK ⇥ C

B̆>
K

R̃K ⇥ T

IK ⇥ R̃K

(a) CIFA for multi-block matrices
(1)

X1

G1

B̆1

B̄(1)

XK

T

...

...

...

=

T

data tensor Y (can be incomplete) is modeled as
Y = X + N + E,

(34)

I1 ⇥ R̃1

...

...

ĀK

B̆>
1
R̃1 ⇥ T

+

I1 ⇥ C

XK

IK ⇥ T

Ă1

C⇥T

=
I1 ⇥ T

B̄>

(1)

B̆K

GK

=

(b) CIFA for multi-block tensors

Fig. 6: Illustration of common and individual feature
extraction (CIFA). (a) The matrix case: the matrix B̄ denotes the common (or highly correlated) features shared
by all data blocks while B̆k are the individual features
specific for each individual block. (b) The tensor case:
CIFA is performed with respect to mode-1, where the
(1)
matrix B̄(1) denotes the common features while B̆k
are the individual features. The multilinear structure of
Gk makes the tensor model more general and flexible.

lenges if the CP model is used to model the data,
because even determining the CP-rank of a given
tensor is NP-hard. If the CP-rank is given, polyadic
factorization with missing values has been developed
by employing CP weighted optimization (CPWOPT)
[105] and nonlinear least squares (CPNLS) [106]. The
Bayesian CP developed in [107] can infer the latent CPrank, but it suffers from the weak point that the missing
data are handled in a heuristic way. Recently, a unified
framework was proposed in [98], [99], where the given

X = JΛ; A(1) , A(2) , . . . , A(N ) K, tensor N is the Gaussian
noise, E denotes the outliers. To infer the latent CPrank R, at the very beginning, R will be assigned a
sufficiently large value and then it shrinks gradually
during the iteration process by enforcing column-wise
sparsity of factor matrices. This method is attractive due
to its full Bayesian treatment which is able to simultaneously handle outliers, Gaussian noise, and incomplete
data uniformly with automatic model order selection.
Particularly, it is free of any tuning parameters.
As a side comment, we would like to emphasize that
outliers are not always treated as noise. For instance, in
video surveillance analysis, outliers can be used to identify the activities that stand out from the background
[99]. Therefore, an alternative interpretation of robust
tensor factorization is to decompose an incomplete tensor as a low-rank tensor, capturing global information
and predicting missing values, and a sparse tensor,
capturing local information. Hence, the power of the
methods we introduced in this section are not limited
to optional pre-processing procedures. Instead, many of
them can be applied independently for different types
of applications. Dimensionality reduction itself is a very
active research topic, many models and approaches
have been proposed for different application purposes.
While readers may find comprehensive reviews in
[108], we summarized some of them for quick reference
in TABLE 3.
It can be seen that the low-rank nature of data is the
foundation of the above models and approaches. This
low-rank nature allows a very compact representation
of large-scale structured raw data, upon which highly
scalable and efficient algorithms can be developed [38],
[109]. To meet the ongoing challenge created by big
data analytics, new data models and optimization algorithms are required. In this regard, tensor networks
were recently proposed to achieve the so-called supercompression (see [110], [111] and references therein).
Randomization and sparsification are emerging techniques that provide favorable scalability for dimensionality reduction. In the first category, the given huge data
tensor can be compressed via either random projections
[112], [113] or random sampling [114]. In the second
category, the dense data tensor is converted to a sparse
one via an element-wise sparsification algorithm [115].
In [116], the authors argued that first-order methods,
randomization, and parallel and distributed computation are three key pillars and offer surprising scalability
benefits for big data analysis. Readers are referred to
[110], [116] and references therein for this ongoing
research topic.

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

6
6.1

I LLUSTRATIVE E XAMPLES

Training	
  data-­‐sets	
  for	
  K	
  mental	
  tasks	
  

Applications of CIFA in SSVEP-based BCI

Steady-state visual evoked potentials (SSVEP) are periodic neural responses elicited at the same frequency as
a visual stimulus when a subject is excited at a specific
frequency by the visual stimulus. An SSVEP-based BCI
is designed to detect the target stimulus by recognizing
the dominant frequency from SSVEP data, and then
to translate it into a pre-assigned computer command.
Although real SSVEPs present relatively stable spectra as long as the recording time is sufficiently long,
they are often contaminated by background sources of
noise. Thus, efficient and accurate recognition of SSVEP
frequencies is essential to achieve high performance in
SSVEP-based BCIs. As one of the most widely adopted
methods, CCA identifies each SSVEP target frequency
by locating the maximum correlation between EEG
recordings and a set of artificial references, each of
which corresponds to a candidate frequency [117].
While CCA gained popularity and provides reasonably good performance in SSVEP frequency recognition, it suffers from two major limitations: it ignores the
multiway structure of EEG data and particularly, uses
artificial references typically constructed by simple sine
and cosine waves rather than real SSVEP features. Recently, by effectively exploiting the multiway feature of
EEG data, a multiway extension of CCA (MWCCA) has
been developed to improve the performance of SSVEP
recognition [118]. The MWCCA approach collaboratively maximizes correlation between multiple dimensions of EEG tensor data and the pre-constructed reference sine-cosine waves at specific stimulus frequencies,
thereby leading to improved performance against the
standard CCA. On the other hand, a group of EEG
data trials recorded at the same stimulus frequency
should share some common features reflecting this frequency information. Such common features extracted
from EEG training data are supposed to bear real
SSVEP characteristics and hence are more qualified as
references for SSVEP recognition (see Fig.8). With this
intuition, MCCA was proposed to optimize the SSVEP
references [82], where the common features extracted
from multiple training sets of EEG data recorded at
the same stimulus frequency were chosen as the references. A solution that does not use any reference
signal but makes use of the statistical dependence of
SSVEP among multiple electrodes has been very recently proposed using IVA and has noted to provide
quite promising performance [119] provided that the
time window is sufficiently long.
Also, we investigated the application of CIFA to learn
real SSVEP features in order to further improve the
accuracy of SSVEP recognition. In this experiment, we
validated effectiveness of CIFA in SSVEP recognition
using the EEG dataset recorded from ten subjects.
During the EEG recordings, four red squares were pre-

15

XK

X1

trial
/me
channel
EEG	
  tensor

Mul/-­‐way	
  CCA,	
  
mul/-­‐set	
  CCA,	
  
CIFA,	
  …

Common
feature
extraction

Common
feature
extraction

CCA

CCA

Extracted	
  
reference	
  signals	
  

Test	
  data
Maximize	
  correla/on

ρK

ρ1

fk= max ρk
Online	
  selected	
  BCI	
  command

Fig. 8: Illustration of SSVEP recognition based on multiway or multiset learning models. Consider multiple
tensors (channel × time × trial) X1 , . . . , XK constructed
by the EEG data corresponding to K distinct target
frequencies of f1 , . . . , fK , respectively. The common features extracted by using MWCCA, MCCA or CIFA are
used as references for SSVEP recognition. The SSVEP
target frequency of new test data for a single trial is then
recognized according to the maximum of correlation
coefficients ρ1 , . . . , ρK between the test data and the
optimized references.

sented on the screen as stimuli flickered with different
frequencies. Each subject completed 20 runs. In each
run, each stimulus was cued as the target followed
by flickering of the four stimuli at four different frequencies: 6, 8, 9 and 10 Hz. Subjects were asked to
focus their attention on the target for 1–4s after each
cue. A total of 80 trials data were acquired from each
subject. More details about the experiment settings
can be found in [118]. Three state-of-the-art methods,
i.e., CCA, MWCCA and MCCA were compared with
CIFA. Fig.9 depicts the averaged recognition accuracies
evaluated by leave-one-run-out cross-validation, from
which we see that CIFA achieved much higher average
accuracy than the others for all time windows from 0.5s

16

PROCEEDINGS OF THE IEEE

SSVEP recognition accuracy (%)

100

Noisy MRI (σ = 10%, SNR = 7dB)

90
80
70
60
CCA
MWCCA
MCCA
CIFA

50
40
30

0.5

1

1.5 2 2.5 3
Time window (s)

3.5

Denoised MRI (PSNR = 36dB)

4

Fig. 9: SSVEP recognition accuracies averaged on subjects, derived by CCA, MWCCA, MCCA and CIFA
respectively, over a time window from 0.5 s to 4 s (0.5
s interval).

to 4s.
6.2

MRI denoising

MRI is a non-invasive high resolution multi-parameter
imaging modality that is widely used in current clinical
diagnosis and scientific research because of its ability to
reveal the 3D internal structures of objects. However,
MRI may suffer from serious effects such as random
noise, specifically when high-resolution and high-speed
is required. The denoising of MRI is thus important
to improve the inspection quality and reliability of
quantitative image analysis. Nonlocal filters by exploiting similarity and sparseness among cubes have
been demonstrated to be effective in denoising MRI.
Recently, HOSVD has been applied to image denoising
by exploiting redundancy in the 3D stack of similar patches [41]. In addition, the augmented HOSVD
method has been investigated and applied to MRI
volume data [120]. One common limitation of existing
denoising methods is that the noise level is assumed
to be known in advance. Therefore, it is appealing to
perform MRI denoising automatically without needing
tuning of the parameters. To achieve this, we apply a recently proposed Bayesian CP factorization (BCPF) [98]
method for MRI denoising. The procedure is described
briefly as follows. Given an I × I × I reference cube
in the noisy MRI data, K similar cubes (including the
reference cube) are found and stacked into a 4D array of
size I ×I ×I ×K, which is thus filtered by BCPF. Finally,
the denoised MRI is obtained by aggregating multiple
estimations at each location. Since BCPF can estimate
CP-rank efficiently and the noise variance automatically
from given data, it can be easily applied to filtering 4D
array instead of HOSVD. For illustration, we use the

Fig. 10: MRI tensor is of size 181 × 217 × 165. The cube
size is 4 × 4 × 4 and the number of similar cubes is 30.
public MRI data1 and consider Gaussian noise whose
standard derivation is 10% of brightest tissue. Fig. 10
shows the results of MRI denoising with an impressive
performance of PSNR= 36dB.
6.3

MRI completion

Missing data are a common problem in MRI studies
due in large part to susceptibility artifact and image
acquisition parameters that result in incomplete brain
coverage and spatial variation in acquired images,
across subjects. In addition, MRI completion can be
also used to improve spatial and temporal resolution
by reconstructing the whole image from very few
measurements. In this experiment, we illustrate MRI
completion by applying the Bayesian tensor completion (BTC) method [100], particularly to improve the
generalization performance when original data do not
posses a global low-rank structure. Several recently
proposed tensor completion methods including the
high accuracy low rank tensor completion (HaLRTC)
method [102], the generalized higher-order orthogonal
iteration (gHOOI) method [121], and the one based on
weighted Tucker (WTucker) decomposition [122] were
also used for comparisons. For HaLRTC and WTucker,
the tuning parameters were selected to yield the best
possible performance. BTC methods as well as gHOOI
can automatically adapt the tensor rank to data. As
shown in Table 4, BTC achieves the best performance
under the conditions studied here. When missing ratio
1. Dataset is available from http://brainweb.bic.mni.mcgill.ca/
brainweb/

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

Estimation

Missing

Noisy

SNR=20dB, MR = 80%, PSNR= 22dB

Fig. 11: Tensor completion for noisy MRI data. The standard derivation of Gaussian noise was 3% of brightest
tissue. MRI tensor is of size 181 × 217 × 165.
TABLE 4: The performance of MRI completion evaluated by
PSNR and RRSE under missing ratios of 50% and 80%.
Methods
BTC
WTucker
gHOOI
HaLRTC

50% missing
PSNR
RRSE
26.40
0.12
23.18
0.18
21.57
0.21
23.12
0.18

80% missing
PSNR
RRSE
22.33
0.19
21.53
0.21
19.77
0.26
17.06
0.36

is low (e.g., 50%), HaLRTC outperforms gHOOI, while
gHOOI outperforms HaLRTC when the missing ratio
is high (e.g., 80%). The visual quality of reconstructed
MRI obtained by BTC is shown in Fig. 11.

7 D ISCUSSION , C ONCLUSIONS
P ERSPECTIVES

AND

F UTURE

Multi-set or multi-block coupled and inhomogeneous
data are becoming ubiquitous across the sciences and
engineering, especially in the biomedical field. Large
scale medical/biological data using modalities such
as EEG, fMRI, ECoG critically require scalable algorithms to efficiently process massive datasets within
reasonable amount of time. To handle these issues,
tensors, multi-dimensional generalizations of matrices,
are attractive and promising tools for the representation
and processing of such data.
In this paper, we have studied analysis of such
multi-relational or multi-block data using simultaneous
(coupled) constrained matrix and tensor factorizations
that allow one to extract common and individual components and/or to establish links between them. We
critically reviewed several useful matrix and tensor
decomposition models for linked multiway analysis (or
multilinear blind source separation).
For large-scale multi-relational data analysis, tensor
decompositions is one emerging technology. Current
applications in many areas of science and technology, such as computational neuroscience, bioinformat-

17

ics, and pattern/image recognition, generate massive
amounts of data with multiple aspects and high dimensionality, and a number of successful examples
have demonstrated that tensor decompositions with
low-rank approximations can provide often useful and
meaningful approximate representations for such large
scale applications [110], [123]. Tensor and matrix decompositions, in particular using the linked models
introduced here, allow us to discover meaningful hidden structures of complex brain data and perform
generalizations by capturing multi-linear and multiaspect relationships. To serve this purpose, flexible
CIFA models were introduced to simultaneously model
common features shared by all data sources and individual features that are specific to each data source. In
BSS, subspace definitions for a joint BSS framework as
recently introduced [124], [125] can be used to define
such common and distinct subspaces within IVA and
the concentrated (hence linked) ICA models like group
and joint ICA. For all those methods, a number of
challenges exist for very large scale data, feature extraction/selection, classification, clustering and anomaly
detection. A recent paper addresses the selection of
model order when the sample size is relatively small by
specifically considering common and distinct subspaces
as in CIFA [93], [126].
Another closely related approach for big data not
discussed in this paper is tensor networks. Complex
interactions and operations between tensors can be
visualized by tensor network diagrams in which tensors
are represented graphically by nodes or shapes [110],
[111]. In such representations, each outgoing edge (line)
emerging from a node represents a mode (a way, a
dimension, indices). Tensor network diagrams could
be potentially very helpful not only in visualizing
tensor decompositions but also for expressing complex
mathematical (multilinear) operations of contractions
of tensors [110], [111]. Tensor networks are connected
to quantum physics, quantum chemistry and quantum
information, which studies the ways to possibly build
a quantum computer and to program it [127]. In addition, tensor networks provide graphically illustrative
large distributed networks and enable one to perform
complex tensor operations (i.e., tensor contractions and
reshaping) in an intuitive way and without the need to
use explicit mathematical representation.
To summarize, the benefits of linked multiway (tensor) analysis methods along with matrix-based approaches for massive biomedical data include:
• A framework to incorporate different type of diversity or various constraints in different modes
or factors, and naturally extending the standard
component analysis and BSS methods to largescale joint BSS as well as multiway data analyses;
• Possibility to handle noisy and missing data naturally by using powerful low-rank tensor approximations and by exploiting robustness and stability

18

PROCEEDINGS OF THE IEEE

of tensor decomposition algorithms;
“Super” compression of large-scale multidimensional data via tensorization and decomposition
of a high-order tensor into factor matrices and/or
core tensors with low-rank and low-order [3], [110];
• Ability to perform all mathematical operations in
feasible tensor network formats, with the ability to
improve the scalability of algorithms [110].
Hence, linked component analysis is a very promising area of research with substantial potentials in
biomedical applications.

[18]

R EFERENCES

[23]

•

[1]

[2]
[3]

[4]
[5]
[6]
[7]
[8]

[9]

[10]

[11]
[12]
[13]

[14]
[15]
[16]

[17]

A. Cichocki, R. Zdunek, A. H. Phan, and S.-I. Amari, Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory
Multi-way Data Analysis and Blind Source Separationpplications to
Exploratory Multi-Way Data Analysis and Blind Source Separation.
Chichester, UK: John Wiley & Sons, Ltd, 2009.
T. G. Kolda and B. W. Bader, “Tensor decompositions and
applications,” SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.
A. Cichocki, D. Mandic, A.-H. Phan, C. Caiafa, G. Zhou,
Q. Zhao, and L. De Lathauwer, “Tensor decompositions for
signal processing application,” IEEE Signal Processing Magazine,
vol. 32, no. 2, pp. 145–163, March 2015.
A. Cichocki, “Tensor decompositions: New concepts for brain
data analysis?” SICE Journal of Control Measurement, and System
Integration, vol. 7, pp. 507–517, 2011.
W. Hackbusch, Tensor Spaces and Numerical Tensor Calculus, ser.
Springer Series in Computational Mathematics. Springer, 2012,
vol. 42.
A. Smilde, R. Bro, and P. Geladi, Multi-Way Analysis with
Applications in the Chemical Sciences. John Wiley & Sons, Ltd,
2004.
P. M. Kroonenberg, Applied Multiway Data Analysis, 1, Ed.
Wiley-Interscience, 2008.
E. Acar, M. Rasmussen, F. Savorani, T. Næs, and R. Bro,
“Understanding data fusion within the framework of coupled
matrix and tensor factorizations,” Chemometrics and Intelligent
Laboratory Systems, vol. 129, pp. 53–63, 2013.
E. Acar, G. Gurdeniz, M. Rasmussen, D. Rago, L. Dragsted,
and R. Bro, “Coupled matrix factorization with sparse factors
to identify potential biomarkers in metabolomics,” in IEEE 12th
International Conference on Data Mining Workshops (ICDMW).
IEEE, 2012, pp. 1–8.
M. Jafari and J. Chambers, “Fetal electrocardiogram extraction
by sequential source separation in the wavelet domain,” IEEE
Transactions on Biomedical Engineering, vol. 52, no. 3, pp. 390–
400, March 2005.
A. Cichocki and S. Amari, Adaptive Blind Signal and Image
Processing: Learning Algorithms and Applications. John Wiley
& Sons, Ltd, 2002.
P. Comon and C. Jutten, Handbook of Blind Source Separation:
Independent Component Analysis and Applications, 1st ed. Academic Press, March 2010.
T. Adali, M. Anderson, and G.-S. Fu, “Diversity in independent
component and vector analyses: Identifiability, algorithms, and
applications in medical imaging,” IEEE Signal Processing Magazine, vol. 31, no. 3, pp. 18–33, May 2014.
T. Adali, H. Li, M. Novey, and J. Cardoso, “Complex ICA using
nonlinear functions,” IEEE Transactions on Signal Processing,
vol. 56, no. 9, pp. 4536–4544, Sept 2008.
J. F. Cardoso, “Blind signal separation: Statistical principles,”
Proceedings of the IEEE, vol. 86, no. 10, pp. 2009–2025, Oct 1998.
V. Calhoun and T. Adali, “Multisubject independent component analysis of fMRI: A decade of intrinsic networks, default
mode, and neurodiagnostic discovery,” IEEE Reviews in Biomedical Engineering, vol. 5, pp. 60–73, 2012.
V. Calhoun and T. Adali, “Unmixing fMRI with independent
component analysis,” IEEE Engineering in Medicine and Biology
Magazine, vol. 25, no. 2, pp. 79–90, March 2006.

[19]
[20]

[21]
[22]

[24]

[25]

[26]
[27]
[28]

[29]

[30]

[31]
[32]

[33]

[34]
[35]

[36]
[37]

[38]

V. Calhoun, T. Adali, G. Pearlson, and J. Pekar, “Spatial and
temporal independent component analysis of functional MRI
data containing a pair of task-related waveforms,” Human Brain
Mapping, vol. 13, no. 1, pp. 43–53, 2001.
Y. Guo and G. Pagnoni, “A unified framework for group
independent component analysis for multi-subject fMRI data,”
NeuroImage, vol. 42, no. 3, pp. 1078 – 1093, 2008.
F. Marvasti, A. Amini, F. Haddadi, M. Soltanolkotabi, B. Khalaj, A. Aldroubi, and J. Sanei, S.and Chambers, “A unified
approach to sparse signal processing,” EURASIP Journal on
Advances in Signal Processing, vol. 2012, no. 44, pp. 1–45, 2012.
S. Sanei and J. A. Chambers, EEG Signal Processing. John Wiley
& Sons, Inc., July 2007.
E. Candes and M. Wakin, “An introduction to compressive
sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp.
21–30, March 2008.
G. Zhou, Z. Yang, S. Xie, and J.-M. Yang, “Mixing matrix
estimation from sparse mixtures with unknown number of
sources,” IEEE Transactions on Neural Networks, vol. 22, no. 2,
pp. 211 –221, Feb. 2011.
A. Ferreol, L. Albera, and P. Chevalier, “Fourth-order blind
identification of underdetermined mixtures of sources FOBIUM,” IEEE Transactions on Signal Processing, vol. 53, no. 5,
pp. 1640–1653, 2005.
L. De Lathauwer, J. Castaing, and J. Cardoso, “Fourth-order
cumulant-based blind identification of underdetermined mixtures,” IEEE Transactions on Signal Processing, vol. 55, no. 6, pp.
2965–2973, 2007.
D. D. Lee and H. S. Seung, “Learning the parts of objects by
non-negative matrix factorization,” Nature, vol. 401, no. 6755,
pp. 788–791, Oct 21 1999.
N. Gillis, “Sparse and unique nonnegative matrix factorization
through data preprocessing,” Journal of Machine Learning Research, vol. 13, no. 11, pp. 3349–3386, Nov 2012.
K. Huang, N. Sidiropoulos, and A. Swami, “Non-negative
matrix factorization revisited: Uniqueness and algorithm for
symmetric decomposition,” IEEE Transactions on Signal Processing, vol. PP(Preprint), no. 99, pp. 1–1, May 2013.
F. Y. Wang, C. Y. Chi, T. H. Chan, and Y. Wang, “Nonnegative
least-correlated component analysis for separation of dependent sources by volume maximization,” IEEE Transactions on
Pattern Analysis and Machine Intelligence, vol. 32, no. 5, pp. 875–
888, May 2010.
G. Zhou, S. Xie, Z. Yang, J.-M. Yang, and Z. He, “Minimumvolume-constrained nonnegative matrix factorization: Enhanced ability of learning parts,” IEEE Transactions on Neural
Networks, vol. 22, no. 10, pp. 1626 –1637, Oct. 2011.
N. Gillis and S. A. Vavasis, “Fast and robust recursive algorithms for separable nonnegative matrix factorization,” ArXiv
e-prints: http://arxiv.org/abs/1208.1237, 2012.
C. Ding, T. Li, W. Peng, and H. Park, “Orthogonal nonnegative
matrix tri-factorizations for clustering,” in 12th ACM SIGKDD
international conference on Knowledge discovery and data mining,
ser. KDD ’06. New York, NY, USA: ACM, 2006, pp. 126–135.
B. Li, G. Zhou, and A. Cichocki, “Two efficient algorithms for
approximately orthogonal nonnegative matrix factorization,”
IEEE Signal Processing Letters, vol. 22, no. 7, pp. 843–846, July
2015.
T. Yokota, Q. Zhao, C. Li, and A. Cichocki, “Smooth PARAFAC
decomposition for tensor completion,” ArXiv e-prints, 2015.
[Online]. Available: http://arxiv.org/abs/1505.06611
T. Yokota, R. Zdunek, A.Cichocki, and Y. Yamashita, “Smooth
nonnegative matrix and tensor factorizations for robust multiway data analysis,” Signal Processing, vol. 113, pp. 234–249,
2015.
S.-I. Amari, A. Cichocki, and H.-H. Yang, “A new learning
algorithm for blind signal separation,” Advances in Neural
Information Processing Systems, vol. 8, pp. 757–763, 1996.
L.-Q. Zhang, A. Cichocki, and S. Amari, “Natural gradient
algorithm for blind separation of overdetermined mixture with
additive noise,” IEEE Signal Processing Letters, vol. 6, no. 11, pp.
293–295, 1999.
G. Zhou, A. Cichocki, Q. Zhao, and S. Xie, “Nonnegative matrix
and tensor factorizations: An algorithmic perspective,” IEEE
Signal Processing Magazine, vol. 31, no. 3, pp. 54–65, May 2014.

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

[39]

[40]

[41]

[42]

[43]

[44]

[45]

[46]
[47]

[48]

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating
direction method of multipliers,” Found. Trends Mach. Learn.,
vol. 3, no. 1, pp. 1–122, 2011.
L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM Journal on Matrix
Analysis and Applications, vol. 21, pp. 1253–1278, 2000.
A. Rajwade, A. Rangarajan, and A. Banerjee, “Image denoising
using the higher order singular value decomposition,” IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 35,
no. 4, pp. 849–862, 2013.
H. A. Phan and A. Cichocki, “Tensor decompositions for feature extraction and classification of high dimensional datasets,”
Nonlinear Theory and Its Applications, IEICE, vol. 1, no. 1, pp. 37–
68, 2010.
L. Omberg, G. Golub, and O. Alter, “A tensor higher-order
singular value decomposition for integrative analysis of DNA
microarray data from different studies,” Proceedings of the National Academy of Sciences, vol. 104, no. 47, pp. 18 371–18 376,
2007.
B. Hunyadi, D. Camps, L. Sorber, W. Van Paesschen, M. De Vos,
S. Van Huffel, and L. De Lathauwer, “Block term decomposition
for modelling epileptic seizures,” EURASIP Journal on Advances
in Signal Processing, vol. 2014, no. 1, pp. 1–19, 2014.
S.Ferdowsi, V.Abolghasemi, and S. Sanei, “EEG-FMRI integration using a partially constrained tensor factorization,” in
IEEE 38th International Conference on Acoustics, Speech and Signal
Processing (ICASSP), Vancouver, BC, Canada, May 2013, pp.
6191–6195.
L. Shoker, S. Sanei, and J. Chambers, “Artifact removal from
electroencephalograms using a hybrid BSS-SVM algorithm,”
IEEE Signal Processing Letters, vol. 12, no. 10, pp. 721–724, 2005.
F. Cong, Q.-H. Lin, L.-D. Kuang, X.-F. Gong, P. Astikainen, and
T. Ristaniemi, “Tensor decomposition of EEG signals: A brief
review,” Journal of Neuroscience Methods, vol. 248, no. 0, pp. 59
– 69, 2015.
M. Weis, F. Römer, M. Haardt, D. Jannek, and P. Husar, “Multidimensional space-time-frequency component analysis of event
related EEG data using closed-form PARAFAC,” in IEEE 34th
International Conference on Acoustics, Speech and Signal Processing
(ICASSP), April 2009, pp. 349–352.
H. Becker, L. Albera, P. Comon, M. Haardt, G. Birot,
F. Wendling, M. Gavaret, C. Bénar, and I. Merlet, “EEG
extended source localization: Tensor-based vs. conventional
methods,” NeuroImage, vol. 96, no. 0, pp. 143 – 157, 2014.
M. De Vos, L. De Lathauwer, B. Vanrumste, S. Van Huffel,
and W. Van Paesschen, “Canonical decomposition of ictal scalp
EEG and accurate source localisation: Principles and simulation
study,” Computational Intelligence and Neuroscience, vol. 2007,
2007.
M. Mørup, L. K. Hansen, and S. M. Arnfred, “ERPWAVELAB: A toolbox for multi-channel analysis of time–frequency
transformed event related potentials,” Journal of Neuroscience
Methods, vol. 161, no. 2, pp. 361 – 368, 2007.
M. Vasilescu and D. Terzopoulos, “Multilinear independent
components analysis,” in IEEE Computer Society Conference on
Computer Vision and Pattern Recognition, vol. 1, June 2005, pp.
547 – 553 vol. 1.
G. Zhou and A. Cichocki, “Fast and unique Tucker decompositions via multiway blind source separation,” Bulletin of the
Polish Academy of Sciences-Technical Sciences, vol. 60, no. 3, pp.
389–405, 12 2012.
N. D. Sidiropoulos and R. Bro, “On the uniqueness of multilinear decomposition of N -way arrays,” Journal of Chemometrics,
vol. 14, no. 3, pp. 229 – 239, 2000.
L. Chiantini, G. Ottaviani, and N. Vannieuwenhoven, “An
algorithm for generic and low-rank specific identifiability of
complex tensors,” SIAM Journal on Matrix Analysis and Applications, vol. 35, no. 4, pp. 1265–1287, 2014.
A. Belouchrani, K. AbedMeraim, J. F. Cardoso, and E. Moulines,
“A blind source separation technique using second-order statistics,” IEEE Transactions on Signal Processing, vol. 45, no. 2, pp.
434–444, Feb 1997.
S. Choi, A. Cichocki, and S. Amari, “Equivariant nonstationary

[58]

[59]

[60]

[61]
[62]

[63]
[64]

[65]
[66]
[67]

[68]

[69]
[70]

[71]

[72]
[73]
[74]
[75]

[76]
[77]

[78]

19

source separation,” Neural Networks, vol. 15, no. 1, pp. 121 –
130, 2002.
A. Ziehe, P. Laskov, G. Nolte, and K. R. Muller, “A fast algorithm for joint diagonalization with non-orthogonal transformations and its application to blind source separation,” Journal
of Machine Learning Research, vol. 5, pp. 777–800, Jul 2004.
L. De Lathauwer, “A link between the canonical decomposition
in multilinear algebra and simultaneous matrix diagonalization,” SIAM Journal on Matrix Analysis and Applications, vol. 28,
no. 3, pp. 642–666, 2006.
G. Zhou, S. Xie, Z. Yang, and J. Zhang, “Nonorthogonal approximate joint diagonalization with well-conditioned diagonalizers,” IEEE Transactions on Neural Networks, vol. 20, no. 11,
pp. 1810–1819, Nov 2009.
C. F. Beckmann and S. M. Smith, “Tensorial extensions of independent component analysis for multisubject fMRI analysis,”
NeuroImage, vol. 25, no. 1, pp. 294 – 311, 2005.
M. De Vos, D. Nion, S. Van Huffel, and L. De Lathauwer,
“A combination of parallel factor and independent component
analysis,” Signal Processing, vol. 92, no. 12, pp. 2990 – 2999,
2012.
G. Zhou and A. Cichocki, “Canonical polyadic decomposition
based on a single mode blind source separation,” IEEE Signal
Processing Letters, vol. 19, no. 8, pp. 523 –526, Aug. 2012.
M. Niknazar, H. Becker, B. Rivet, C. Jutten, and P. Comon,
“Blind source separation of underdetermined mixtures of
event-related sources,” Signal Processing, vol. 101, no. 0, pp. 52
– 64, 2014.
L.-H. Lim and P. Comon, “Nonnegative approximations of
nonnegative tensors,” Journal of Chemometrics, vol. 23, no. 7-8,
pp. 432–441, 2009.
Y. Qi, P. Comon, and L.-H. Lim, “Uniqueness of nonnegative tensor approximations,” eprint arXiv:1410.8129, p.
http://arxiv.org/abs/1410.8129, 2015.
S. Nielsen and M. Morup, “Non-negative tensor factorization
with missing data for the modeling of gene expressions in
the human brain,” in IEEE International Workshop on Machine
Learning for Signal Processing (MLSP). IEEE, 2014, pp. 1–6.
Y. Liang, J. Harris, S. M. Naqvi, G. Chen, and J. A. Chambers,
“Independent vector analysis with a generalized multivariate
Gaussian source prior for frequency domain blind source separation,” Signal Processing, vol. 105, no. 0, pp. 175 – 184, 2014.
V. Calhoun, J. Liu, and T. Adali, “A review of group ICA for
fMRI data and ICA for joint inference of imaging, genetic, and
ERP data,” NeuroImage, vol. 45, no. 1, pp. 163 – 172, 2009.
V. Calhoun, T. Adali, G. Pearlson, and K. Kiehl, “Neuronal
chronometry of target detection: Fusion of hemodynamic and
event-related potential data,” NeuroImage, vol. 30, no. 2, pp. 544
– 553, 2006.
V. Calhoun, T. Adali, G. Pearlson, and J. Pekar, “A method for
making group inferences from functional MRI data using independent component analysis,” Human Brain Mapping, vol. 14,
no. 3, pp. 140–151, 2001.
N. Correa, T. Adali, Y.-O. Li, and V. Calhoun, “Canonical
correlation analysis for data fusion and group inferences,” IEEE
Signal Processing Magazine, vol. 27, no. 4, pp. 39–50, 2010.
J. R. Kettenring, “Canonical analysis of several sets of variables,” Biometrika, vol. 58, no. 3, pp. 433–451, 1971.
A. Groves, C. Beckmann, S. Smith, and M. Woolrich, “Linked
independent component analysis for multimodal data fusion,”
NeuroImage, vol. 54, no. 1, pp. 2198 – 21 217, 2011.
N. Helwig and S. Hong, “A critique of tensor probabilistic
independent component analysis: implications and recommendations for multi-subject fMRI data analysis,” Journal of neuroscience methods, vol. 213, no. 2, pp. 263–273, 2013.
R. S. Schaefer, P. Desain, and J. Farquhar, “Shared processing
of perception and imagery of music in decomposed {EEG},”
NeuroImage, vol. 70, no. 0, pp. 317 – 326, 2013.
C. Crainiceanu, B. S. Caffo, S. Luo, V. M. Zipunnikov, and N. M.
Punjabi, “Population value decomposition, a framework for the
analysis of image populations,” Journal of the American Statistical
Association, vol. 106, no. 495, pp. 775–790, 2011.
G. Zhou, A. Cichocki, and S. Xie, “Accelerated canonical
polyadic decomposition by using mode reduction,” IEEE Trans-

20

[79]

[80]
[81]

[82]

[83]

[84]

[85]

[86]

[87]

[88]

[89]

[90]

[91]

[92]

[93]

[94]

[95]
[96]
[97]

PROCEEDINGS OF THE IEEE

actions on Neural Networks and Learning Systems, vol. 24, no. 12,
pp. 2051–2062, Dec. 2013.
M. Anderson, G. Fu, R. Phlypo, and T. Adali, “Independent vector analysis: Identification conditions and performance
bounds,” IEEE Transactions on Signal Processing, vol. 62, no. 17,
pp. 4399–4410, Sept 2014.
H. Hotelling, “Relations between two sets of variates,”
Biometrika, vol. 28, no. 3-4, pp. 321–377, Dec 1936.
Y.-O. Li, T. Adali, W. Wang, and V. Calhoun, “Joint blind source
separation by multiset canonical correlation analysis,” IEEE
Transactions on Signal Processing, vol. 57, no. 10, pp. 3918 – 3929,
October 2009.
Y. Zhang, G. Zhou, J. Jin, X. Wang, and A. Cichocki, “Frequency
recognition in SSVEP-based BCI using multiset canonical correlation analysis,” International Journal of Neural Systems, vol. 4,
no. 24, pp. 1 450 013(1–14), 2014.
P. Acharjee, R. Phlypo, L. Wu, V. Calhoun, and T. Adali,
“Independent vector analysis for gradient artifact removal in
concurrent EEG-fMRI data,” IEEE Transactions on Biomedical
Engineering, vol. to appear, no. 99, pp. 1–1, 2015.
Q. Zhao, G. Zhou, T. Adali, L. Zhang, and A. Cichocki, “Kernelization of tensor-based models for multiway data analysis:
Processing of multidimensional structured data,” IEEE Signal
Processing Magazine, vol. 30, no. 4, pp. 137–148, 2013.
T. Adali, Y. Levin-Schwartz, and V. Calhoun, “Multimodal data
fusion using source separation: Two effective models based on
ICA and IVA and their properties,” Proceedings of the IEEE,
vol. PP, no. 99, pp. 1–16, 2015.
A. Michael, M. Anderson, R. Miller, T. Adali, and V. D. Calhoun, “Preserving subject variability in group fMRI analysis:
Performance evaluation of GICA versus IVA,” Frontiers in
Systems Neuroscience, vol. 8, no. 106, pp. 1–18, 2014.
S. Ma, R. Phlypo, V. Calhoun, and T. Adali, “Capturing group
variability using IVA: A simulation study and graph-theoretical
analysis,” in IEEE 38th International Conference on Acoustics,
Speech and Signal Processing (ICASSP), May 2013, pp. 3128–3132.
J. Laney, K. P. Westlake, S. Ma, E. Woytowicz, V. Calhoun, and
T. Adali, “Capturing subject variability in fMRI data: A graphtheoretical analysis of GICA vs. IVA,” Journal of Neuroscience
Methods, vol. 247, pp. 32–40, 2015.
S. Gopal, R. Miller, A. Michael, T. Adali, S. Baum, and V. Calhoun, “A study of spatial variation in fMRI brain networks via
independent vector analysis: Application to schizophrenia,” in
International Workshop on Pattern Recognition in Neuroimaging,
June 2014, pp. 1–4.
Q. Zhao, C. Caiafa, D. Mandic, Z. Chao, Y. Nagasaka, N. Fujii,
L. Zhang, and A. Cichocki, “Higher order partial least squares
(HOPLS): A generalized multilinear regression method,” IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 35,
no. 7, pp. 1660–1673, July 2013.
E. Lock, K. Hoadley, J. S. Marron, and A. Nobel, “Joint and
individual variation explained (JIVE) for integrated analysis
of multiple data types,” The Annals of Applied Statistics, vol. 7,
no. 1, pp. 523–542, 2013.
G. Zhou, A. Cichocki, and D. Mandic, “Common components
analysis via linked blind source separation,” in IEEE 40th
International Conference on Acoustics, Speech and Signal Processing
(ICASSP). IEEE, April 2015.
G. Zhou, A. Cichocki, Y. Zhang, and D. Mandic, “Group
component analysis from multi-block data: Common and
individual feature extraction,” IEEE Transactions on Neural
Networks and Learning Systems, Accepted. [Online]. Available:
http://arxiv.org/abs/1212.3913
Z. He, A. Cichocki., S. Xie, and K. Choi, “Detecting the number
of clusters in N -way probabilistic clustering,” IEEE Transactions
on Pattern Analysis and Machine Intelligence, vol. 32, no. 11, pp.
2006–2021, 2010.
E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal
component analysis?” Journal of the ACM (JACM), vol. 58, no. 3,
pp. 11:1–11:37, 2011.
E. Candes and Y. Plan, “Matrix completion with noise,” Proceedings of the IEEE, vol. 98, no. 6, pp. 925–936, June 2010.
N. Shahid, V. Kalofolias, X. Bresson, M. Bronstein, and P. Vandergheynst, “Robust principal component analysis on graphs,”
ArXiv e-prints: http://arxiv.org/abs/1504.06151, 2015.

[98]

[99]

[100]
[101]

[102]

[103]
[104]

[105]
[106]

[107]

[108]
[109]

[110]
[111]
[112]
[113]

[114]
[115]
[116]

[117]

[118]

Q. Zhao, L. Zhang, and A. Cichocki, “Bayesian CP factorization
of incomplete tensors with automatic rank determination,”
IEEE Transactions on Pattern Analysis and Machine Intelligence,
vol. 37, no. 9, pp. 1751–1763, 2015.
Q. Zhao, G. Zhou, L. Zhang, A. Cichocki, and S. Amari,
“Bayesian robust tensor factorization for incomplete multiway
data,” IEEE Transactions on Neural Networks and Learning Systems, vol. PP, no. 99, pp. 1–1, 2015.
Q. Zhao, L. Zhang, and A. Cichocki, “Bayesian sparse Tucker
models for dimension reduction and tensor completion,” arXiv
preprint arXiv:1505.02343, 2015.
B. Jiang, C. Ding, B. Luo, and J. Tang, “Graph-Laplacian PCA:
Closed-form solution and robustness,” in IEEE Conference on
Computer Vision and Pattern Recognition (CVPR 2013), June 2013,
pp. 3492–3498.
J. Liu, P. Musialski, P. Wonka, and J. Ye, “Tensor completion for
estimating missing values in visual data,” IEEE Transactions on
Pattern Analysis and Machine Intelligence, vol. 35, no. 1, pp. 208–
220, Jan 2013.
D. Goldfarb and Z. T. Qin, “Robust low-rank tensor recovery:
Models and algorithms,” SIAM Journal on Matrix Analysis and
Applications, vol. 35, no. 1, pp. 225–253, 2014.
Z. Xu, F. Yan, and A. Qi, “Infinite Tucker decomposition:
Nonparametric bayesian models for multiway data analysis,”
in 29th International Conference on Machine Learning (ICML),
J. Langford and J. Pineau, Eds. New York, NY, USA: ACM,
2012, pp. 1023–1030.
E. Acar, D. M. Dunlavy, T. G. Kolda, and M. Mørup, “Scalable
tensor factorizations for incomplete data,” Chemometrics and
Intelligent Laboratory Systems, vol. 106, no. 1, pp. 41–56, 2011.
L. Sorber, M. V. Barel, and L. D. Lathauwer, “Optimizationbased algorithms for tensor decompositions: Canonical
polyadic decomposition, decomposition in rank-$(l r,l r,1)$
terms, and a new generalization,” SIAM Journal on Optimization,
vol. 23, no. 2, pp. 695–720, 2013.
P. Rai, Y. Wang, S. Guo, G. Chen, D. Dunson, and L. Carin,
“Scalable Bayesian low-rank decomposition of incomplete multiway tensors,” in 31st International Conference on Machine Learning (ICML), 2014, pp. 1800–1808.
Y. Ma, P. Niyogi, G. Sapiro, and R. Vidal, “Dimensionality
reduction via subspace and submanifold learning,” IEEE Signal
Processing Magazine, vol. 28, no. 2, pp. 14–126, March 2011.
G. Zhou, A. Cichocki, and S. Xie, “Fast nonnegative matrix/tensor factorization based on low-rank approximation,”
IEEE Transactions on Signal Processing, vol. 60, no. 6, pp. 2928
–2940, 2012.
A. Cichocki, “Tensor networks for big data analytics and largescale optimization problems,” ArXiv e-prints, Jul. 2014.
A. Cichocki, “Era of big data processing: A new approach via
tensor networks and tensor decompositions,” ArXiv e-prints,
2014. [Online]. Available: http://arxiv.org/abs/1403.2048
M. W. Mahoney, “Randomized algorithms for matrices and
data,” Found. Trends Mach. Learn., vol. 3, no. 2, pp. 123–224,
2011.
N. Halko, P. Martinsson, and J. Tropp, “Finding structure with
randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM Review, vol. 53, no. 2, pp.
217–288, 2011.
C. F. Caiafa and A. Cichocki, “Generalizing the column-row
matrix decomposition to multi-way arrays,” Linear Algebra and
its Applications, vol. 433, no. 3, pp. 557 – 573, 2010.
N. H. Nguyen, P. Drineas, and T. D. Tran, “Tensor sparsification via a bound on the spectral norm of random tensors,”
Information and Inference, pp. 1–35, 2015.
V. Cevher, S. Becker, and M. Schmidt, “Convex optimization for
big data: Scalable, randomized, and parallel algorithms for big
data analytics,” IEEE Signal Processing Magazine, vol. 31, no. 5,
pp. 32–43, Sept 2014.
Z. Lin, C. Zhang, W. Wu, and X. Gao, “Frequency recognition
based on canonical correlation analysis for SSVEP-based BCIs,”
IEEE Transactions on Biomedical Engineering, vol. 54, no. 6, pp.
1172–1176, June 2007.
Y. Zhang, G. Zhou, J. Jin, M. Wang, X. Wang, and A. Cichocki,
“L1-regularized multiway canonical correlation analysis for

ZHOU et al.: LINKED COMPONENT ANALYSIS AND ITS APPLICATIONS TO BIOMEDICAL DATA

[119]

[120]

[121]

[122]
[123]

[124]

[125]

[126]

[127]

SSVEP-based BCI,” IEEE Transactions on Neural Systems and
Rehabilitation Engineering, vol. 21, no. 6, pp. 887–896, 2013.
D. Emge, F.-B. Vialatte, G. Dreyfus, and T. Adali, “Independent
vector analysis for SSVEP signal enhancement,” in 49th Annual
Conference on Information Sciences and Systems (CISS), March
2015, pp. 1–6.
X. Zhang, Z. Xu, N. Jia, W. Yang, Q. Feng, W. Chen, and Y. Feng,
“Denoising of 3D magnetic resonance images by using higherorder singular value decomposition,” Medical Image Analysis,
vol. 19, no. 1, pp. 75–86, 2015.
Y. Liu, F. Shang, W. Fan, J. Cheng, and H. Cheng, “Generalized
higher-order orthogonal iteration for tensor decomposition and
completion,” in Advances in Neural Information Processing Systems (NIPS), 2014, pp. 1763–1771.
M. Filipović and A. Jukić, “Tucker factorization with missing
data with application to low-n-rank tensor completion,” Multidimensional Systems and Signal Processing, pp. 1–16, 2013.
E. Papalexakis, T. Mitchell, N. Sidiropoulos, C. Faloutsos, P. P.
Talukdar, and B. Murphy, “Turbo-SMT: Accelerating coupled
sparse matrix-tensor factorizations by 200x,” in SIAM International Conference on Data Mining, 2014.
D. Lahat and C. Jutten, “Joint blind source separation of
multidimensional components: Model and algorithm,” in 22nd
European Signal Processing Conference (EUSIPCO), Sept 2014, pp.
1417–1421.
R. Silva, S. Plis, T. Adali, and V. Calhoun, “Multidataset
independent subspace analysis extends independent vector
analysis,” in IEEE International Conference on Image Processing
(ICIP), Oct 2014, pp. 2864–2868.
Y. Song, P. Schreier, and N. J. Roseveare, “Determining the
number of correlated signals between two data sets using PCACCA when sample support is extremely small,” in IEEE 40th
International Conference on Acoustics, Speech and Signal Processing
(ICASSP). IEEE, April 2015.
R. Orús, “A practical introduction to tensor networks: Matrix
product states and projected entangled pair states,” Annals of
Physics, vol. 349, no. 0, pp. 117 – 158, 2014.

21

