Lecture Notes in Computer Science

Fig. 4. The Hinton diagram of the first and the second projection matrices obtained by PTA, as 

shown in the left and the right sub-figures, respectively 

with the following requirements: 

1

2

3

6

, ,

10

l l l

≤

≤

, 

1

2

3

2

, ,

5

l l l

′ ′ ′

≤

≤

, 

1

1

l

l

′ <

, 

2

2

l

l

′ <

, and 

3

3

l

l

′ <

. The total model selection errors are also 0. In every experiment, the value of 

the reconstruction error is very small. 

4   Conclusion 

Vector 

data are normally used for probabilistic graphical models with probabilistic 

inference. However, 

tensor 

data, i.e., multidimensional arrays, are actually natural 

representations of a large amount of data, in data mining, computer vision, and many 

other applications. Aiming at breaking the huge gap between 

vectors 

and 

tensors 

in 

conventional statistical tasks, e.g., model selection, this paper proposes a decoupled 

probabilistic algorithm, named 

probabilistic tensor analysis 

(PTA) with 

Akaike 

information criterion 

(AIC) and 

Bayesian information criterion 

(BIC). PTA 

associated AIC and BIC can select suitable models for 

tensor 

data, as demonstrated 

by empirical studies. 

Acknowledgment 

We authors would like to thank Professor Andrew Blake at the Microsoft Research 

Cambridge for encouragement of developing the tensors probabilistic graphic model. 

This research was supported by the Competitive Research Grants at the Hong Kong 

Polytechnic University (under project number A-PH42 and A-PC0A) and the 

National Natural Science Foundation of China (under grant number 60703037). 

References 

[1]  Bader, B.W., Kolda, T.G.: Efficient MATLAB Computations with Sparse and Factored 

Tensors. Technical Report SAND2006-7592, Sandia National Laboratories, Albuquerque, 

NM and Livermore, CA (2006) 

[2]  Bader, B.W., Kolda, T.G.: MATLAB Tensor Classes for Fast Algorithm Prototyping. 

ACM Transactions on Mathematical Software 32(4) (2006) 

 

Probabilistic Tensor Analysis with Akaike and Bayesian Information Criteria 

801 

[3]  Bishop, C.M.: Neural Networks for Pattern Recognition. Oxford University Press, Oxford 

(1995) 

[4]  Kroonenberg, P., Leeuw, J.D.: Principal Component Analysis of Three-Mode Data by 

Means of Alternating Least Square Algorithms. Psychometrika, 45 (1980) 

[5]  Lathauwer, L.D.: Signal Processing Based on Multilinear Algebra, Ph.D. Thesis. 

Katholike Universiteit Leuven (1997) 

[6]  Sun, J., Tao, D., Faloutsos, C.: Beyond Streams and Graphs: Dynamic Tensor Analysis. 

In: The Twelfth ACM SIGKDD International Conference on Knowledge Discovery and 

Data Mining, Philadelphia, PA, USA, pp. 374–383 (2006) 

[7]  Tao, D., Li, X., Wu, X., Maybank, S.J.: General Tensor Discriminant Analysis and Gabor 

Features for Gait Recognition. IEEE Transactions on Pattern Analysis and Machine 

Intelligence 29(10) (2007) 

[8]  Tipping, M.E., Bishop, C.M.: Probabilistic Principal Component Analysis. Journal of the 

Royal Statistical Society, Series B 21(3), 611–622 (1999) 

[9]  Tucker, L.R.: Some Mathematical Notes on Three-mode Factor Analysis. 

Psychometrika 31(3) (1966) 

[10]  Vasilescu, M.A.O., Terzopoulos, D.: Multilinear Subspace Analysis of Image Ensembles. 

In: IEEE Proc. International Conference on Computer Vision and Pattern Recognition, 

Madison, Wisconsin, USA, vol. 2, pp. 93–99 (2003) 

[11]  Ye, J., Janardan, R., Li, Q.: Two-Dimensional Linear Discriminant Analysis. In: 

Schölkopf, B., Platt, J., Hofmann, T. (eds.) Advances in Neural Information Processing 

Systems, Vancouver, British Columbia, Canada, vol. 17 (2004) 

[12]  Ye, J.: Generalized Low Rank Approximations of Matrices. Machine Learning 61(1-3), 

167–191 (2005) 

M. Ishikawa et al. (Eds.): ICONIP 2007, Part I, LNCS 4984, pp. 802–810, 2008. 

© Springer-Verlag Berlin Heidelberg 2008 

Decomposing EEG Data into Space-Time-Frequency 

Components Using Parallel Factor Analysis and Its 

Relation with Cerebral Blood Flow 

Fumikazu Miwakeichi

1

, Pedro A. Valdes-Sosa

3

, Eduardo Aubert-Vazquez

3

, 

Jorge Bosch Bayard

3

, Jobu Watanabe

4

, Hiroaki Mizuhara

2

, and Yoko Yamaguchi

2

  

1

 Medical System Course, Graduate School of Engineering, Chiba University 

1-33, Yayoi-cho, Inage-ku, Chiba-shi, Chiba, 263-8522 Japan  

[e-mail] 

2 

Laboratory for Dynamics of Emergent Intelligence,  

RIKEN Brain Science Institute, Japan 

3 

Cuban Neuroscience Center, Cuba 

4 

Waseda Institute for Advance Study, Waseda University, Japan

 

Abstract. Finding the means to efficiently summarize electroencephalographic 

data has been a long-standing problem in electrophysiology. Our previous 

works showed that Parallel Factor Analysis (PARAFAC) can effectively per-

form atomic decomposition of the time-varying EEG spectrum in space/ 

frequency/time domain. In this study, we propose to use PARAFAC for extract-

ing significant activities in EEG data that is concurrently recorded with func-

tional Magnetic Resonance Imaging (fMRI), and employ the temporal signature 

of the atom for investigating the relation between brain electrical activity and 

the changing of BOLD signal that reflects cerebral blood flow. We evaluated 

the statistical significance of dynamical effect of BOLD respect to EEG based 

on the modeling of BOLD signal by plain autoregressive model (AR), its AR 

with exogenous EEG input (ARX) and ARX with nonlinear term (ARNX). 

Keywords:  Parallel Factor Analysis, EEG space/frequency/time decomposi-

tion, Nonlinear time series analysis, Concurrent fMRI/EEG data. 

1   Introduction 

The electroencephalogram (EEG) is recorded as multi-channel time-varying data. In 

the history of EEG study, there are so many types of oscillatory phenomena in spon-

taneous and evoked EEG have been observed and reported. A statistical description of 

the oscillatory phenomena of the EEG was carried out first in the frequency domain 

by estimation of the power spectrum for quasi-stationary segments of data [1]. More 

recent characterizations of transient oscillations are carried out by estimation of the 

time-varying (or evolutionary) spectrum in the frequency/time domain [2]. These 

evolutionary spectra of EEG oscillations will have a topographic distribution on the 

sensors that is contingent on the spatial configuration of the neural sources that gener-

ate them as well as the properties of the head as a volume conductor [3].  

  Decomposing EEG Data into Space-Time-Frequency Components Using PARAFAC 

803 

There is a long history of atomic decompositions for the EEG. However, to date, 

atoms have not been defined by the triplet spatial, spectral and temporal signatures 

but rather pairwise combinations of these components.  

Space/time atoms are the basis of both Principal Component Analysis (PCA) and 

Independent Component Analysis (ICA) as applied to multi-channel EEG. PCA has 

been used for artifact removal and to extract significant activities in the EEG [4,5]. A 

basic problem is that atoms defined by only two signatures (space and time) are not 

determined uniquely. In PCA orthogonality is therefore imposed between the corre-

sponding signatures of different atoms. And there is the well-known non-uniqueness 

of PCA that allows the arbitrary choice of rotation of axes (e.g. Varimax and Quarti-

max rotations). More recently, ICA has become a popular tool for space/time atomic 

decomposition [6,7]. It avoids the arbitrary choice of rotation (Jung et al. 2001). 

Uniqueness, however, is achieved at the price of imposing a constraint even stronger 

than orthogonality, namely, statistical independence. In both PCA and ICA the fre-

quency information may be obtained from the temporal signature of the extracted 

atoms in a separate step. 

For the purpose of decomposing of single channel EEG into frequency/time atoms 

the Fast Fourier Transformation (FFT) with sliding window [8] or the wavelet 

transformation [9,10] have been employed. In fact, any of the frequency/time atomic 

decompositions currently available [11] could, in principle, be used for the EEG. How-

ever, these methods do not address the topographic aspects of the EEG time/frequency 

analysis. 

It has long been known, especially in the chemometrics literature, that unique mul-

tilinear decompositions of multi-way arrays of data (more than 2 dimensions) are 

possible under very weak conditions [12]. In fact, this is the basic argument for Paral-

lel Factor Analysis (PARAFAC). This technique recently has been improved by Bro 

[13]. The important difference between PARAFAC and techniques such as PCA or 

ICA, is that the decomposition of multi-way data is unique even without additional 

orthogonality or independence constraints. 

Thus, PARAFAC can be employed for a space/frequency/time atomic decomposition 

of the EEG. This makes use of the fact that multi-channel evolutionary spectra are 

multi-way arrays, indexed by electrode, frequency and time. The inherent uniqueness of 

the PARAFAC solution leads to a topographic time/frequency decomposition with a 

minimum of a priori assumptions. It has been shown that PARAFAC can effectively 

perform a time/frequency/spatial (T/F/S) atomic decomposition which is suitable for 

identification of fundamental modes of oscillatory activity in the EEG [14,15]. 

In this paper, the theory of PARAFAC and its applications to EEG analysis will be 

showed. Moreover the possibility of analysis for investigating the relation between 

extracted EEG atom and cerebral blood flow will be discussed. 

2   Theory 

For the purpose of EEG analysis, we define the 

d

f

t

N

N

N

×

×

 data  matrix 

S

 as  the 

three-way time-varying EEG spectrum array obtained by applying a wavelet trans-

formation, where 

,

d

f

N

N

 and 

t

N

 are the numbers of channels, frequency steps and 

804 

F. Miwakeichi et al. 

time points, respectively. For the wavelet transformation a complex Morlet mother 

function was used. The energy 

d f  t

S

 of  channel 

d

 at  frequency 

f

 and  time 

t

 is 

given by the squared norm of the convolution of a Morlet wavelet with the EEG sig-

nal 

( , )

v d t

,  

2

( , )

( , )

d f  t

S

w t f

v d t

=

∗

,

 

(1) 

where the complex Morlet wavelet 

( , )

w t f

 

is defined by 

2

2

( , )

b

t

i

f t

b

w t f

e

e

σ

π

πσ

⎛

⎞

−

⎜

⎟

⎝

⎠

=

 

with 

b

σ

 denoting the bandwidth parameter. The 

width of the wavelet, 

2

b

m

f

πσ

=

, is set to 7 in this study. 

The basic structural model for a PARAFAC decomposition of the data matrix 

(

)

Nd Nf Nt

× ×

S

 with elements 

 

 

d f t

S

 is defined by 

 

 

 

 

 

 

 

 

 

 

 

1

ˆ

k

N

d f t

d k

f k t k

d f t

d f t

d f t

k

S

a b c

e

S

e

=

=

+

=

+

∑

.

 

(2) 

The problem is to find the loading matrices 

A

, 

B

 and 

C

, the elements of which 

are denoted by 

 

 

,

d k

f k

a

b

 and 

 

t k

c

 in Eq. (2). Here we will refer to components 

k

 as 

“atoms”, and the corresponding loading vectors 

{ }

{ }

{ }

,

,

k

dk

k

fk

k

tk

a

b

c

=

=

=

a

b

c

 will be said to represent the spatial, spectral 

and temporal signatures of the atoms (Fig. 1). The uniqueness of the decomposition 

(2) is guaranteed if 

rank( )+rank( )+rank( )

2(

1)

k

N

≥

+

A

B

C

. As can be seen, this is 

a less stringent condition than either orthogonality or statistical independence [12]. 

 

Fig. 1. Graphical explanation of the PARAFAC model. The multi-channel EEG evolutionary 

spectrum 

S

 is obtained from a channel by channel wavelet transform. 

S

 is a three-way data 

array indicated by channel, frequency and time. PARAFAC decomposes this array into the sum 

of “atoms”. The k-th atom is the trilinear product of loading vectors representing spatial (

k

a

), 

spectral (

k

b

) and temporal (

k

c

) “signatures”. Under these condition PARAFAC can be sum-

marized as finding the matrices 

{ }

k

=

A

a

, 

{ }

k

=

B

b

 and 

{ }

k

=

C

c

 that  explain 

S

 with 

minimal residual error. 

  Decomposing EEG Data into Space-Time-Frequency Components Using PARAFAC 

805 

The decomposition (2) can be obtained by evaluating 

 

 

 

 

 

 

 

 

 

 

2

 

 

ˆ

min

d k

f k

t k

d f t

d f t

a

b

c

S

S

−

. Since 

the 

 

 

d f t

S

 can be regarded as representing spectra, this minimization should be carried 

out under a non-negativity constraint for the loading vectors. The result of the 

PARAFAC decomposition is given by the  (

1)

d

N

×

 vector 

k

a

, representing the topog-

raphical map of the 

k

-th component, the 

(

1)

f

N

×

vector 

k

b

 representing  the  fre-

quency spectrum of the 

k

-th component, and the (

1)

t

N

×

 vector 

k

c

 representing the 

temporal signature of the 

k

-th component.  

3   Results 

3.1   Extracting Relevant Components from EEG

 

The experiment consisted of two conditions, eyes-closed resting condition and mental 

arithmetic condition; for each condition five epochs, each lasting 30 s per condition, 

were recorded. During the mental arithmetic epochs, the subjects were asked to count 

backwards from 1000 in steps of a constant single-digit natural number, while keep-

ing their eyes closed. At the beginning of each mental arithmetic epoch this number 

was randomly chosen by a computer and presented to the subjects through head-

phones. The end of each mental arithmetic epoch was announced by a beeping sound 

(4 kHz, 20 ms). Mental arithmetic and resting epochs were occurring five times 

within each trial and arranged in alternating order; each subject was examined in two 

trials.  

During the experiment, fMRI and EEG were recorded simultaneously. The elec-

trode set consisted of 61 EEG channels, two ECG channels and one EOG channel. 

The reference electrode was at FCz. Raw EEG signals were sampled at a sampling 

frequency of 5 kHz, using a 1Hz high-pass filter provided by the Brain Vision Re-

corder software (Brain Products, Munich, Germany). The fMRI was acquired as 

blood-oxygenation-sensitive (T2*-weighted) echoplanar images, using a 1.5 T MR 

scanner (Staratis II, Hitachi Medico, Japan) and a standard head coil. Thirty slices 

(4mm in thickness, gapless), covering almost the entire cerebrum, were measured, 

under the following conditions: TR 5 s, TA 3.3 s, TE 47.2 ms, FA 901, FoV 240mm, 

matrix size 64×64. In order to minimize head motion, heads of subjects were immobi-

lized by a vacuum pad. SPM99 was used for preprocessing of fMRI images (motion 

correction, slice timing correction). Using bilinear interpolation, images were normal-

ized for each subject to a standard brain defined by Montreal Neurological Institute; 

the normalized images were subsequently smoothed using a Gaussian kernel (full-

width half-maximum: 10 mm). For each voxel, the BOLD signal time series were 

high-pass filtered (>0.01Hz), in order to remove fluctuations with frequencies lower 

than the frequency defined by the switching between resting and task. 

The EEG signal

( , )

v d t

 was subsampled to 500Hz, and, in order to construct a 

three-way (channel/frequency/time) time-varying EEG 

d

f

t

N

N

N

×

×

 data matrix 

S

, a 

 

806 

F. Miwakeichi et al. 

 

Fig. 2. 

(a) Spectral signatures of atoms of Parallel Factor Analysis (PARAFAC) for a typical 

subject. Note the recurrent appearance of frequency peaks in the theta and alpha bands. The 

horizontal axis represents frequency in Hz, the vertical axis represents the normalized ampli-

tude. (b) Spatial signatures of atoms displayed as a topographic map, for the theta, alpha and 

high alpha atoms (from above) of Parallel Factor Analysis (PARAFAC) for the same subject. 

(c) Temporal signatures of atoms, same order as in (a). The horizontal axis represents time in 

units of scans (time between scans is five seconds), the vertical axis represents normalized 

intensity. The black and red colored lines are corresponding to resting and task stages, respec-

tively. (See colored figure in CD-R). 

wavelet transform was applied in the frequency range from 0.5 to 40.0 Hz with 0.5 Hz 

steps, using the complex Morlet wavelet as mother function. For adjusting the temporal 

resolution to the resolution of the fMRI data, the wavelet-transformed EEG was aver-

aged over consecutive 5 seconds intervals. Then, by applying PARAFAC to this three-

way data set, major signal components in the frequency range from 3.0 to 30.0 Hz were 

  Decomposing EEG Data into Space-Time-Frequency Components Using PARAFAC 

807 

extracted. The reason for choosing this frequency band was that at lower frequency 

EEG data typically contains eye movement artifacts, while at higher frequency there is 

no relevant activity in the data set analyzed in this study. In Fig.2(a) spectral signatures 

of the identified atoms are shown: theta (around 6Hz), alpha (around 9Hz), and high 

alpha (around 12Hz) atoms were found; in this study we denote the frequency of the last 

atom as “high alpha”. We find that if the order of the decomposition is increased beyond 

three, two atoms with overlapping frequency peaks appear. Since this indicates overfit-

ting, we select three as the order of the decomposition. In Fig. 2(b) spatial signatures of 

atoms are shown. It can be seen that the main power of theta and alpha atoms is focused 

in frontal and occipital areas, respectively. In Fig.2(c) temporal signatures of atoms are 

shown; red color refers to task condition, black color to resting condition. By comparing 

the amplitudes of temporal signatures during both conditions by a standard T-test, we 

find that for the theta atom the amplitudes of the temporal signature are higher during 

the task condition, while for the alpha atom they are higher during the resting condition. 

For the high alpha atom no significant result is found with this test, therefore this atom 

will be omitted from further study.