Lecture Notes in Computer Science

3.2   Hemodynamic Response of Cerebral Blood Flow with Respect to EEG 

In order to investigate the dynamical relation between EEG and cerebral blood flow, 

we considered the linear autoregressive model (AR), the linear autoregressive model 

with exogenous input (ARX) and nonlinear nonparametric term (ARNX). The defini-

tions of these models are as follows; 

AR                                          

1

t

t

t

B

B

e

α

−

+

=

,                                                     (3) 

ARX                                       

0

1

m

t

t

k

t k

t

k

p

B

B

x

e

β

β

−

−

∈

=

+

+

∑

U

,                              (4) 

ARNX    

(

)

{

}

1

2

0

1

1

2

,

,

;    

,

,

m

t

t

k

t k

t j

t j

t

k

p

np

B

B

x

f x

x

e

j j

β

β

−

−

−

−

∈

=

+

+

+

⊆

∑

…

…

U

U

   (5) 

 

where the

( )

.

f

 is a nonparametric Nadaraya-Watson type kernel estimator and the 

sequence of 

t

e

 is prediction error. 

The indices of those columns of the design matrix that belong to the parametric or 

nonparametric parts were fixed for ARX and ARNX at the same set of values, 

{

}

1, 2, 3, 4, 5, 6

p

np

=

=

U

U

. The reason for choosing a maximum lag of 6 is given by 

the fact that the effect of the neural electric activity on the BOLD signal has been 

experimentally found to decay to zero level in approximately 30 seconds.  

In ARX and ARNX there is an autoregressive (AR) model term, with parame-

ter

0

β

, describing intrinsic dynamics of the BOLD signal, in addition to the terms 

representing the influence of the neural electric activity. 

The parameters of the parametric and the nonparametric model part, as well as the 

kernel width parameter, can be estimated by minimizing the modified Global Cross 

Validation (GCV) criterion (Hong, 1999). It can also be used for model comparison. 

Consider the difference between the values of GCV for two models, A and B: 

808 

F. Miwakeichi et al. 

( )

( )

GCV

-GCV

D

A

B

=

                                            (6) 

If for a particular voxel D assumes a positive value, model B is superior to A. 

If for a particular voxel the GCV value for ARX or ARNX is smaller than for 

AR, it can be said that this voxel is influenced by the EEG. From now on, we will 

focus only on the set of these voxels. We expect that among these voxels there is a 

subset of voxels for which the nonlinear term of ARNX improves the model, as 

compared to ARX. For extracting these voxels, 

(

)

(

)

GCV ARX -GCV ARNX

D

=

 

was computed. If for a particular voxel, D  assumes a positive value, this indicates 

the necessity of employing a nonlinear response to the EEG for modeling the BOLD 

signal of this voxel. On the other hand, if D assumes a negative value, it is sufficient 

to model the BOLD signal of this voxel by ARX. 

Regions where ARX or ARNX outperforms AR, are shown in Fig.3 (a) for the al-

pha atom and in Fig.3 (b) for the theta atom; note that these regions contain only 7% 

(alpha atom) and 4% (theta atom) of all gray-matter voxels. Regions where ARNX 

outperforms ARX (as measured by GCV) are shown by green color; these regions 

contain about 3% (alpha atom) and 2% (theta atom) of all gray-matter voxels. The 

opposite case of regions where ARNX fails to outperform ARX is shown by orange 

color; these regions contain about 4% (alpha atom) and 2% (theta atom) of all gray-

matter voxels.  

 

Fig. 3.

 

Regions where the ARX and ARNX models outperform the AR model, for the alpha 

atom (a) and the theta atom (b). Green color denotes regions where ARNX outperforms ARX, 

while orange color denotes the opposite case (thresholded at p<0.05 FDR). (See colored figure 

in CD-R). 

4   Discussion 

For the analysis of concurrently recorded EEG/fMRI data set, essentially it needs 

multivariate versus multivariate analysis. However, there has been no established 

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

809 

effective method for this purpose. Therefore this problem has been resolved into mul-

tivariate versus unitivariate analysis using representative single electrode or averaged 

EEG of all or selected electrodes.  

Though temporal signature of atom is univariate, it is summarization using all in-

formation across all electrodes. In this study, it was shown that multivariate EEG data 

can be effectively summarized using PARAFAC. And the temporal signatures of 

atoms can be employed as reference function for investigating dynamical response of 

cerebral blood flow respect to brain electrical activity. The underlying theoretical 

requirement is that of a moderate amount of linear independence for atom topogra-

phies, spectra and time courses. This is a much milder requirement than previous 

models underlying space/time atomic decompositions (PCA or ICA).  

From the physiological point of view, we note that those regions, where the ARX 

or ARNX models outperformed the AR model for the alpha atom, approximately 

correspond to areas where the BOLD response was found to be negatively correlated 

to the alpha rhythm [16,17,18]; this situation indicates a decrease of cortical hemody-

namic activity during the resting state. Previous studies reported bilateral premotor, 

posterior parietal, and prefrontal activation during serial subtraction [19,20]. This 

result agrees well with the finding of the present study, according to which the voxels 

in these areas were correlated to the theta atom. In summary we conclude that inclu-

sion of linear and nonlinear EEG input terms improves modeling of brain activation in 

these regions during both rest and mental arithmetic. 

References 

1.  Lopes da Silva, F.: EEG Analysis: Theory and Practice. In: Neidermeyer, E., Lopes da Silva, 

F., Electroencephalography. Urban and Schwartzenberg (1987) 

2.  Dahlhaus, R.: Fitting time series models to non-stationary processes. Annals of Statistics 25, 

1–37 (1997) 

3.  Nunez, P.L.,: Electric Fields of the Brian: The Neurophysics of Eeg. Oxford University 

Press, Oxford (1993) 

4.  Soong, A.C., Koles, Z.J.: Principal-component localization of the sources of the background 

EEG. IEEE Transactions on Biomedical Engineering 42, 59–67 (1995) 

5.  Lagerlund, T.D., Sharbrough, F.W., Busacker, N.E.: Spatial Filtering of Multichannel Elec-

troencephalographic Recordings Through Principal Component Analysis by Singular Value 

Decomposition. Journal of Clinical Neurophysiology 14, 73–83 (1997) 

6.  Hyvarinen, A., Karhunen, J., Oja, E.: Independent Component Analysis. John Wiley & Sons, 

Inc., Chichester (2001) 

7.  Cichocki, A., Amari, S.: Adaptive Blind Signal and Image Processing. John Wiley & Sons, 

LTD., Chichester (2002) 

8.  Makeig, S.: Auditory event-related dynamics of the EEG spectrum and effects of exposure to 

tones. Electroencephalography and Clinical Neurophysiology 86, 283–293 (1993) 

9.  Bertrand, O., Bohorquez, J., Pernier, J.: Time-Frequency Digital Filtering Based on an In-

vertible Wavelet Transform: An Application to Evoked Potential. IEEE Transactions on 

Biomedical Engineering 41, 77–88 (1994) 

10.  Tallon-Baudry, C., Bertrand, O., Delpuech, C., Pernier, J.: Oscillatory g - Band (30-70Hz) 

Activity Induced by a Visual Search Task in Humans. The Journal of Neuroscience 17, 722–

734 (1997) 

810 

F. Miwakeichi et al. 

11.  Chen, S., Donoho, D.: Atomic decomposition by basis pursuit. SIAM Review 43, 129–159 

(2001) 

12.  Sidiropoulos, N.D., Bro, R.: On the uniqueness of multilinear decomposition of N-way ar-

rays. Journal of Chemometrics 14, 229–239 (2000) 

13.  Bro, R.: Multi-way Analysis in the Food Industry: Models, Algorithms and Applications. 

Ph.D. Thesis. University of Amsterdam (NL) & Royal Veterinary and Agricultural Univer-

sity (DK) (1998) 

14.  Miwakeichi, F., Martinez-Montes, E., Valdes-Sosa, P.A., Nishiyama, N., Mizuhara, H., Ya-

maguchi, Y.: Decomposing EEG data into space/time/frequency components using Parallel 

Factor Analysis. NeuroImage 22, 1035–1045 (2004) 

15.  Martinez-Montes, E., Valdes-Sosa, P.A., Miwakeichi, F., Goldman, R.I., Cohen, M.S.: Con-

current EEG/fMRI analysis by multiway Partial Least Squares. NeuroImage 22, 1023–1034 

(2004) 

16.  Goldman, R.I., Stern, J.M., Engel, J.J., Cohen, M.S., Simultaneous, E.E.G.: fMRI of the al-

pha rhythm. NeuroReport 13, 2487–2492 (2002) 

17.  Laufs, H., Kleinschmidt, A., Beyerle, A., Eger, E., Salek-Haddadi, A., Preibisch, C., Krakow, 

K.: EEG-correlated fMRI of human alpha activity. NeuroImage 19, 1463–1476 (2003) 

18.  Laufs, H., Krakow, K., Sterzer, P., Eger, E., Beyerle, A., Salek-Haddadi, A., Kleinschmidt, 

A.: Electroencephalographic signatures of attentional and cognitive default modes in sponta-

neous brain activity fluctuations at rest. Proc.Natl.Acad.Sci. 100, 11053–11058 (2003) 

19.  Roland, P.E., Friberg, L.: Localization of Cortical Areas Activated By Thinking. 

J.Neurophysiol, 53 (1985) 

20.  Rueckert, L., Lange, N., Partiot, A., Appollonio, I., Litvan, I., Bihan, D.L., Grafman, J.: 

Visualizing Cortical Activation during Mental Calculation with Functional MRI. Neuroi-

mage 3, 97–103 (2006) 

Flexible Component Analysis for Sparse,

Smooth, Nonnegative Coding or Representation

Andrzej Cichocki , Anh Huy Phan, Rafal Zdunek , and Li-Qing Zhang

RIKEN Brain Science Institute, Wako-shi, Saitama, Japan

[e-mail]

Abstract. In the paper, we present a new approach to multi-way Blind

Source Separation (BSS) and corresponding 3D tensor factorization that

has many potential applications in neuroscience and multi-sensory or

multidimensional data analysis, and neural sparse coding. We propose

to use a set of local cost functions with flexible penalty and regularization

terms whose simultaneous or sequential (one by one) minimization via

a projected gradient technique leads to simple Hebbian-like local algo-

rithms that work well not only for an over-determined case but also (un-

der some weak conditions) for an under-determined case (i.e., a system

which has less sensors than sources). The experimental results confirm

the validity and high performance of the developed algorithms, especially

with usage of the multi-layer hierarchical approach.

1

Introduction – Problem Formulation

Parallel Factor analysis (PARAFAC) or multi-way factorization models with

sparsity and/or non-negativity constraints have been proposed as promising and

quite efficient tools for processing of signals, images, or general data

[1,2,3,4,5,6,7,8,9]. In this paper, we propose new hierarchical alternating algo-

rithms referred to as the Flexible Component Analysis (FCA) for BSS, including

as special cases: Nonnegative Matrix/Tensor Factorization (NMF/NTF), SCA

(Sparse Components Analysis), SmoCA (Smooth Component Analysis). The

proposed approach can be also considered as an extension of Morphological

Component Analysis (MoCA) [10]. By incorporating nonlinear projection or fil-

tering and/or by adding regularization and penalty terms to the local squared

Euclidean distances, we are able to achieve nonnegative and/or sparse and/or

smooth representations of the desired solution, and to alleviate a problem of

getting stuck in local minima.

In this paper, we consider quite a general factorization related to the 3D

PARAFAC2 model [1,5] (see Fig.1)

Y

q

=

AD

q

X

q

+

N

q

=

AX

q

+

N

q

,

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

(1)

Dr. A. Cichocki is also with Systems Research Institute (SRI), Polish Academy of

Science (PAN), and Warsaw University of Technology, Dept. of EE, Warsaw, Poland.

Dr. R. Zdunek is also with Institute of Telecommunications, Teleinformatics and

Acoustics, Wroclaw University of Technology, Poland.

Dr. L.-Q. Zhang is with the Shanghai Jiaotong University, China.

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

c Springer-Verlag Berlin Heidelberg 2008

812

A. Cichocki et al.

where

Y

q

= [y

itq

]

∈ R

I×T

is the q-th frontal slice (matrix) of the observed

(known) 3D tensor data or signals

Y ∈ R

I×T ×Q

,

D

q

∈ R

J×J

+

is a diagonal scal-

ing matrix that holds the q-th row of the matrix

D ∈ R

Q×J

,

A = [a

ij

] =

[

a

1

, a

2

, . . . , a

J

]

∈ R

I×J

is a mixing or basis matrix,

X

q

= [˜

x

jtq

]

∈ R

J×T

represents unknown normalized sources or hidden components in q-th slice,

X

q

=

D

q

X

q

= [x

jtq

]

∈ R

J×T

represents re-normalized (differently scaled)

sources, and

N

q

= [n

itq

]

∈ R

I×T

represents the q-th frontal slice of the tensor

N ∈ R

I×T ×Q

representing noise or errors. Our objective is to estimate the set

of all matrices:

A, D

q

, X

q

, subject to some natural constraints such as non-

negativity, sparsity or smoothness. Usually, the common factors, i.e., matrices

A and X

q

are normalized to unit length column vectors and rows, respectively,

and are often enforced to be as independent and/or as sparse as possible.

The above system of linear equations can be represented in an equivalent

scalar form as follows y

itq

=

j

a

ij

x

jtq

+ n

itq

, or equivalently in the vector form

Y

q

=

j

a

j

x

jq

+

N

q

, where

x

jq

= [x

j1q

, x

j2q

, . . . , x

jT q

] are the rows of

X

q

,

and

a

j

are the columns of

A (j = 1, 2, . . . , J). Moreover, using the row-wise

unfolding, the model (1) can be represented by one single matrix equation:

Y = AX + N,

(2)

where

Y = [Y

1

, Y

2

, . . . , Y

Q

]

∈ R

I×QT

,

X = [X

1

, X

2

, . . . , X

Q

]

∈ R

J×QT

, and

N = [N

1

, N

2

, . . . , N

Q

]

∈ R

I×QT

are block row-wise unfolded matrices

1

. In the

special case, for Q = 1 the model simplifies to the standard BSS model used

in ICA, NMF, and SCA . The majority of the well-known algorithms for the

PARAFAC models work only if the following assumption T >> I

≥ J is held,

where J is known or can be estimated using PCA/SVD. In the paper, we pro-

pose a family of algorithms that can work also for an under-determined case,

i.e., for T >> J > I, if sources are enough sparse and/or smooth. Our pri-

mary objective is to estimate the mixing (basis) matrix

A and the sources X

q

,

subject to additional natural constraints such as nonnegativity, sparsity and/or

smoothness constraints. To deal with the factorization problem (1) efficiently, we

adopt several approaches from constrained optimization, regularization theory,

multi-criteria optimization, and projected gradient techniques. We minimize si-

multaneously or sequentially several cost functions with the desired constraints

using switching between two sets of parameters:

{A} and {X

q

}.

1

It should be noted that the 2D unfolded model, in a general case, is not exactly

equivalent to the PARAFAC2 model (in respect to sources

X

q

), since we usually

need to impose different additional constraints for each slice

q. In other words,

the PARAFAC2 model should not be considered as a 2-D model with the single

2-D unfolded matrix

X. Profiles of the augmented (row-wise unfolded) X can

only be treated as a single profile, while we need to impose individual constraints

independently to each slice

X

q

or even to each row of

X

q

. Moreover, the 3D tensor

factorization is considered as a dynamical process or a multi-stage process, where

the data analysis is performed several times under different conditions (initial

estimates, selected natural constraints, post-processing) to get full information

about the available data and/or discover some inner structures in the data, or to

extract physically meaningful components.

FCA for Sparse, Smooth, Nonnegative Coding or Representation

813

Q

Q

N

A

Q

J

I

T

T

=

+

...

...

Y

1

i

q

1

T

t

1

1

1

1

J

D

(

)

I T Q

x

x

(

)

I J

x

(

)

J T Q

x

x

(

)

I T Q

x

x

J

I

I

T

q

X

~

Q

X

~

1

X

~

Fig. 1. Modified PARAFAC2 model illustrating factorization of 3D tensor into a set

of fundamental matrices:

A, D, {X

q

}. In the special case, the model is reduced to

standard PARAFAC for

X

q

=

X

1

, ∀q, or tri-NMF model for Q = 1.

2

Projected Gradient Local Least Squares Regularized

Algorithm

Many algorithms for the PARAFAC model are based on Alternating Least

Square (ALS) minimization of the squared Euclidean distance [1,4,5]. In par-

ticular, we can attempt to minimize a set of the following cost functions:

D

F q

(

Y

q

||AX

q

) =

1

2

Y

q

− AX

q

2

F

+ α

A

J

A

(

A) + α

X

J

X

q

(

X

q

),

(3)

subject to some additional constraints, where J

A

(

A), J

X

q

(

X

q

) are penalty or

regularization functions, and α

A

and α

X

are nonnegative coefficients controlling

a tradeoff between data fidelity and

a priori knowledge on the components to

be estimated. A choice of the regularization terms can be critical for attaining

a desired solution and noise suppression. Some of the candidates include the

entropy, l

p

-quasi norm and more complex, possibly non-convex or non-smooth

regularization terms [11]. In such a case a basic approach to the above formulated

optimization problem is alternating minimization or alternating projection: the

specified cost function is alternately minimized with respect to two sets of the

parameters

{x

jtq

} and {a

ij

}, each time optimizing one set of arguments while

keeping the other one fixed [6,7].

In this paper, we consider a different approach: instead of minimizing only one

global cost function, we perform sequential minimization of the set of local cost

functions composed of the squared Euclidean terms and regularization terms:

D

(

j)

F q

(

Y

(

j)

q

||a

j

x

jq

) =

1

2

(

Y

(

j)

q

− a

j

x

jq

)

2

F

+ α

(

j)

a

J

a

(

a

j

) + α

(

j)

X

q

J

x

(

x

jq

),

(4)

for j = 1, 2, . . . , J, q = 1, 2, . . . , Q, subject to additional constraints, where

Y

(

j)

q

=

Y

q

−

r=j

a

r

x

rq

=

Y

q

− AX

q

+

a

j

x

jq

,

(5)

814

A. Cichocki et al.

a

j

∈ R

I×1

are the columns of the basis mixing matrix

A, x

jq

∈ R

1

×T

are the

rows of

X

q

which represent unknown source signals, J

a

(

a

j

) and J

x

(

x

jq

) are local

penalty regularization terms which impose specific constraints for the estimated

parameters, and α

(

j)

a

≥ 0 and α

(

j)

X

q

≥ 0 are nonnegative regularization parameters

that control a tradeoff between data-fidelity and the imposed constraints.

The construction of such a set of local cost functions follows from the sim-

ple observation that the observed data can be decomposed as follows

Y

q

=

J

j=1

a

j

x

jq

+

N

q

, ∀ q. We are motivated to use of such a representation and

decomposition, because

x

jq

have physically meaningful interpretation as sources

with specific temporal and morphological properties.

The penalty terms may take different forms depending on properties of the

estimated sources. For example, if the sources are sparse, we can apply the l

p

-

quasi norm: J

x

(

x

jq

) =

||x

jq

||

p

p

= (

t

|x

jtq

|

p

)

1

/p

with 0 < p

≤ 1, or alternatively

we can use the smooth approximation J

x

(

x

jq

) =

T

t

|x

jtq

|

2

+ ε

p/2

, where

ε ≥ 0 is a small constant. In order to impose local smoothing of signals, we can

apply the total variation (TV) J

x

(

x

jq

) =

T −1

t=1

|x

j,t+1,q

− x

j,t,q

|, or if we wish

to achieve a smoother solution: J

x

(

x

jq

) =

T −1

t=1

|x

j,t+1,q

− x

j,t,q

|

2

+ ε, [12].

The gradients of the local cost function (4) with respect to the unknown

vectors

a

j

and

x

jq

are expressed by

∂D

(

j)

F q

(

Y

(

j)

q

||a

j

x

jq

)

∂x

jq

=

a

T

j

a

j

x

jq

− a

T

j

Y

(

j)

q

+ α

(

j)

X

q

Ψ

x

(

x

jq

),

(6)

∂D

(

j)

F q

(

Y

(

j)

q

||a

j

x

jq

)

∂a

j

=

a

j

x

jq

x

T

jq

− Y

(

j)

q

x

T

jq

+ α

(

j)

a

Ψ

a

(

a

j

),

(7)

where the matrix functions Ψ

a

(

a

j

) and Ψ

x

(

x

jq

) are defined as

2

Ψ

a

(

a

j

) =

∂J

(

j)

a

(

a

j

)

∂a

j

,

Ψ

x

(

x

jq

) =

∂J

(

j)

x

(

x

jq

)

∂x

jq

.

(8)

By equating the gradient components to zero, we obtain a new set of local

learning rules:

x

jq

←

1

a

T

j

a

j

a

T

j

Y

(

j)

q

− α

(

j)

X

q

Ψ

x

(

x

jq

) ,

(9)

a

j

←

1

x

jq

x

T

jq

Y

(

j)

q

x

T

jq

− α

(

j)

a

Ψ

a

(

a

j

) ,

(10)

for j = 1, 2, . . . , J and q = 1, 2, . . . , Q.

However, it should be noted that the above algorithm provides only a reg-

ularized least squares solution, and this is not sufficient to extract the desired

2

If the penalty functions are non-smooth, we can use sub-gradient instead of the

gradient.

FCA for Sparse, Smooth, Nonnegative Coding or Representation

815

sources, especially for an under-determined case since the problem may have

many solutions. To solve this problem, we need additionally to impose nonlinear

projections P

Ω

j

(

x

jq

) or filtering after each iteration or each epoch in order to

enforce that individual estimated sources

x

jq

satisfy the desired constraints. All

such projections or filtering can be imposed individually for the each source

x

jq

depending on morphological properties of the source signals. The similar nonlin-

ear projection P

Ω

j

(

a

j

) can be applied, if necessary, individually for each vector

a

j

of the mixing matrix

A. Hence, using the projected gradient approach, our

algorithm can take the following more general and flexible form:

x

jq

←

1

a

T

j

a

j

(

a

T

j

Y

(

j)

q

− α

(

j)

X

q

Ψ

x

(

x

jq

)),

x

jq

← P

Ω

j

{x

jq

},

(11)

a

j

←

1

x

jq

x

T

jq

(

Y

(

j)

q

x

T

jq

− α

(

j)

a

Ψ

a

(

a

j

)),

a

j

← P

Ω

j

{a

j

};

(12)

where P

Ω

j

{x

jq

} means generally a nonlinear projection, filtering, transforma-

tion, local interpolation/extrapolation, inpainting, smoothing of the row vector

x

jq

. Such projections or transformations can take many different forms depend-

ing on required properties of the estimated sources (see the next section for more

details).

Remark 1. In practice, it is necessary to normalize the column vectors a

j

or

the row vectors

x

jq

to unit length vectors (in the sense of the l

p

norm (p =

1, 2, ...,

∞)) in each iterative step. In the special case of the l

2

norm, the above

algorithm can be further simplified by neglecting the denominator in (11) or in

(12), respectively. After estimating the normalized matrices

A and X

q

(i.e., the

normalized

X

q

to unit-length rows), we can estimate the diagonal matrices, if

necessary, as follows:

D

q

= diag

{A

+

Y

q

X

+

q

},

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

(13)

3

Flexible Component Analysis (FCA) – Possible

Extensions and Practical Implementations

The above simple algorithm can be further extended or improved (in respect to

a convergence rate and performance). First of all, different cost functions can be

used for estimating the rows of the matrices

X

q

(q = 1, 2, . . . , Q) and the columns

of the matrix

A. Furthermore, the columns of A can be estimated simultaneously,

instead one by one. For example, by minimizing the set of cost functions in (4)

with respect to

x

jq

, and simultaneously the cost function (3) with normalization

of the columns

a

j

to an unit l

2

-norm, we obtain a new FCA learning algorithm

in which the individual rows of

X

q

are updated locally (row by row) and the

matrix

A is updated globally (all the columns a

j

simultaneously):

x

jq

← a

T

j

Y

(

j)

q

− α

(

j)

X

q

Ψ

x

(x

jq

),

x

jq

← P

Ω

j

{x

jq

},

(j = 1, . . . , J ), (14)

A ← (Y

q

X

T

q

− α

A

Ψ

A

(

A))(X

q

X

T

q

)

−1

, A ← P

Ω

(

A), (q = 1, . . . , Q),(15)

816

A. Cichocki et al.

with the normalization (scaling) of the columns of

A to an unit length in the

sense of the l

2

norm, where Ψ

A

(

A) = ∂J

A

(

A)/∂A.

In order to estimate the basis matrix

A, we can use alternatively the following

global cost function (see Eq. (2)): D

F

(

Y ||AX) = (1/2) Y −AX

2

F

+α

A

J

A

(

A).

The minimization of the cost function for a fixed

X leads to the updating rule

A ← Y X

T

− α

A

Ψ

A

(

A) (XX

T

)

−1

.

(16)

3.1

Nonnegative Matrix/Tensor Factorization

In order to enforce sparsity and nonnegativity constraints for all the parameters:

a

ij

≥ 0, x

jtq

≥ 0, ∀ i, t, q, we can apply the ”half-way rectifying” element-wise

projection: [

x]

+

= max

{ε, x}, where ε is a small constant to avoid numerical

instabilities and remove background noise (typically, ε = [10

−2

− 10

−9

]). Si-

multaneously, we can impose weak sparsity constriants by using the l

1

-norm

penalty functions: J

A

(

A) = ||A||

1

=

ij

a

ij

and J

(

j)

x

(

x

jq

) =

||x

jq

||

1

=

t

x

jtq

.

In such a case, the FCA algorithm for the 3D NTF2 (i.e., the PARAFAC2 with

nonnegativity constraints) will take the following form:

x

jq

← a

T

j

Y

(

j)

q

− α

(

j)

X

q

1

+

,

(j = 1, . . . , J ), (q = 1, . . . , Q),

(17)

A ← (Y X

T

− α

A

1)(

XX

T

)

−1

+

,

(18)

with normalization of the columns of

A in each iterative step to a unit length

with the l

2

norm, where 1 means a matrix of all ones of appropriate size.

It should be noted the above algorithm can be easily extended to semi-NMF

or semi-NTF in which only some sources

x

jq

are nonnegative and/or the mixing

matrix

A is bipolar, by simple removing the corresponding ”half-wave rectifying”

projections. Moreover, the similar algorithm can be used for arbitrary bounded

sources with known lower and/or upper bounds (or supports), i.e l

jq

≤ x

jtq

≤

u

iq

, ∀t, rather than x

jtq

≥ 0, by using suitably chosen nonlinear projections

which bound the solutions.

3.2

Smooth Component Analysis (SmoCA)

In order to enforce smooth estimation of the sources

x

jq

for all or some pre-

selected indexes j and q, we may apply after each iteration (epoch) the local

smoothing or filtering of the estimated sources, such as the MA (Moving Aver-

age), EMA, SAR or ARMA models.

A quite efficient way of smoothing and denoising can be achieved by minimiz-

ing the following cost function (which satisfies multi-resolution criterion):

J(x

jq

) =

T

t=1

(x

jtq

− x

jtq

)

2

+

T −1

t=1

λ

jtq

g

t

(x

j,t+1,q

− x

jtq

) ,

(19)

where x

jtq

is a smoothed version of the actually estimated (noisy) x

jtq

, g

t

(u) is

a convex continuously differentiable function with a global minimum at u = 0,

and λ

jtq

are parameters that are data driven and chosen automatically.

FCA for Sparse, Smooth, Nonnegative Coding or Representation

817

3.3

Multi-way Sparse Component Analysis (MSCA)

In the sparse component analysis an objective is to estimate the sources

x

jq

which are sparse and usually with a prescribed or specified sparsification profile,

possibly with additional constraints like local smoothness. In order to enforce

that the estimated sources are sufficiently sparse, we need to apply a suitable

nonlinear projection or filtering which allows us adaptively to sparsify the data.

The simplest nonlinear projection which enforces some sparsity to the normalized

data is to apply the following weakly nonlinear element-wise projection:

P

Ω

j

(x

jtq

) = sign(x

jtq

)

|x

jtq

|

1+

α

jq

(20)

where α

jq

is a small parameter which controls sparsity. Such nonlinear projection

can be considered as a simple (trivial) shrinking. Alternatively, we may use

more sophisticated adaptive local soft or hard shrinking in order to sparsify

individual sources. Usually, we have the three-steps procedure: First, we perform

the linear transformation:

x

w

=

xW , then, the nonlinear shrinking (adaptive

thresholding), e.g., the soft element-wise shrinking: P

Ω

(

x

w

) = sign(

x

w

) [

|x

w

| −

δ]

1+

δ

+

, and finally the inverse transform:

x = P

Ω

(

x

w

)

W

−1

. The threshold δ > 0

is usually not fixed but it is adaptively (data-driven) selected or it gradually

decreases to zero with iterations. The optimal choice for a shrinkage function

depends on a distribution of data. We have tested various shrinkage functions

with gradually decreasing δ: the hard thresholding rule, soft thresholding rule,

non-negative Garrotte rule, n-degree Garotte, and Posterior median shrinkage

rule [13]. For all of them, we have obtained the promising results, and usually

the best performance appears for the simple hard rule.

Our method is somewhat related to the MoCA and SCA algorithms, pro-

posed recently by Bobin et al., Daubechies et al., Elad et al., and many others

[10,14,11]. However, in contrast to these approaches our algorithms are local

and more flexible. Moreover, the proposed FCA is more general than the SCA,

since it is not limited only to a sparse representation via shrinking and linear

transformation but allows us to impose general and flexible (soft and hard) con-

straints, nonlinear projections, transformations, and filtering

3

. Furthermore, in

the contrast to many alternative algorithms which process the columns of

X

q

,

we process their rows which represent directly the source signals.

We can outline the FCA algorithm as follows:

1. Set the initial values of the matrix

A and the matrices X

q

, and normalize

the vectors

a

j

to an unit l

2

-norm length.

2. Calculate the new estimate

x

jq

of the matrices

X

q

using the iterative formula

in (14).

3. If necessary, enforce the nonlinear projection or filtering by imposing natural

constraints on individual sources (the rows of

X

q

, (q = 1, 2, . . . , Q)), such

as nonnegativity, boundness, smoothness, and/or sparsity.

3

In this paper, in fact, we use two kinds of constraints: the soft (or weak) con-

straints via penalty and regularization terms in the local cost functions, and the

hard (strong) constraints via iteratively adaptive postprocessing using nonlinear

projections or filtering.

818

A. Cichocki et al.

4. Calculate the new estimate of

A from (16), normalize each column of A to

an unit length, and impose the additional constraints on

A, if necessary.

5. Repeat the steps (2) and (4) until the convergence criterion is satisfied.

3.4

Multi-layer Blind Identification

In order to improve the performance of the FCA algorithms proposed in this pa-

per, especially for ill-conditioned and badly-scaled data and also to reduce the

risk of getting stuck in local minima in non-convex alternating minimization,

we have developed the simple hierarchical multi-stage procedure [15] combined

together with multi-start initializations, in which we perform a sequential de-

composition of matrices as follows. In the first step, we perform the basic de-

composition (factorization)

Y

q

≈ A

(1)

X

(1)

q

using any suitable FCA algorithm

presented in this paper. In the second stage, the results obtained from the first

stage are used to build up a new tensor

X

1

from the estimated frontal slices de-

fined as

Y

(1)

q

=

X

(1)

q

, (q = 1, 2, . . . , Q). In the next step, we perform the similar

decomposition for the new available frontal slices:

Y

(1)

q

=

X

(1)

q

≈ A

(2)

X

(2)

q

,

using the same or different update rules. We continue our decomposition taking

into account only the last achieved components. The process can be repeated

arbitrarily many times until some stopping criteria are satisfied. In each step, we

usually obtain gradual improvements of the performance. Thus, our FCA model

has the following form:

Y

q

≈ A

(1)

A

(2)

· · · A

(

L)

X

(

L)

q

, (q = 1, 2, . . . , Q) with the

final components

A = A

(1)

A

(2)

· · · A

(

L)

and

X

q

=

X

(

L)

q

.

Physically, this means that we build up a system that has many layers or

cascade connections of L mixing subsystems. The key point in our approach is

that the learning (update) process to find the matrices

X

(

l)

q

and

A

(

l)

is performed

sequentially, i.e. layer by layer, where each layer is initialized randomly. In fact,

we found that the hierarchical multi-layer approach plays a key role, and it is

necessary to apply in order to achieve satisfactory performance for the proposed

algorithms.

4

Simulation Results

The algorithms presented in this paper have been tested for many difficult bench-

marks for signals and images with various temporal and morphological properties

of signals and additive noise. Due to space limitation we present here only one

illustrative example. The sparse nonnegative signals with different sparsity and

smoothness profiles are collected in with 10 slices

X

q

(Q = 10) under the form

of the tensor

X ∈ R

5

×1000×10

. The observed (mixed) 3D data

Y ∈ R

4

×1000×10

are obtained by multiplying the randomly generated mixing matrix

A ∈ R

4

×5

by

X. The simulation results are illustrated in Fig. 2 (only for one frontal slice

q = 1).

FCA for Sparse, Smooth, Nonnegative Coding or Representation

819

Fig. 2. (left) Original 5 spectra (representing X

1

); (middle) observed 4 mixed spec-

tra

Y

1

generated by random matrix

A ∈ R

4×5

(under-determined case); (right) es-

timated 5 spectra

X

1

with our algorithm given by (17)–(18), using 10 layers, and

α

A

=

α

(j)

X

1

= 0

.05. Signal-to-Interference Ratios (SIR) for A and X

1

are as fol-

lows:

SIR

A

= 31

.6, 34, 31.5, 29.9, 23.1[dB] and SIR

X

1

= 28

.5, 19.1, 29.3, 20.3, 23.2[dB],

respectively

5

Conclusions and Discussion

The main objective and motivations of this paper is to derive simple algo-

rithms which are suitable both for under-determined (over-complete) and over-

determined cases. We have applied the simple local cost functions with flexible

penalty or regularization terms, which allows us to derive a family of robust FCA

algorithms, where the sources may have different temporal and morphological

properties or different sparsity profiles. Exploiting these properties and

a pri-

ori knowledge about the character of the sources we have proposed a family of

efficient algorithms for estimating sparse, smooth, and/or nonnegative sources,

even if the number of sensors is smaller than the number of hidden components,

under the assumption that the some information about morphological or desired

properties of the sources is accessible.

This is an original extension of the standard MoCA and NMF/NTF algo-

rithms, and to the authors’ best knowledge, the first time such algorithms have

been applied to the multi-way PARAFAC models. In comparison to the ordinary

BSS algorithms, the proposed algorithms are shown to be superior in terms of

the performance, speed, and convergence properties. We implemented the dis-

cussed algorithms in MATLAB [16]. The approach can be extended for other

applications, such as dynamic MRI imaging, and it can be used as an alterna-

tive or improved reconstruction method to: the k-t BLAST, k-t SENSE or k-t

SPARSE, because our approach relaxes the problem of getting stuck to in local

minima, and provides usually better performance than the standard FOCUSS

algorithms.

This research is motivated by potential applications of the proposed models

and algorithms in three areas of the data analysis (especially, EEG and fMRI)

and signal/image processing: (i) multi-way blind source separation, (ii) model

820

A. Cichocki et al.

reduction and selection, and (iii) sparse image coding. Our preliminary experi-

mental results are promising. The models can be further extended by imposing

additional, natural constraints such as orthogonality, continuity, closure, uni-

modality, local rank - selectivity, and/or by taking into account a prior knowledge

about the specific components.

References

1. Smilde, A., Bro, R., Geladi, P.: Multi-way Analysis: Applications in the Chemical

Sciences. John Wiley and Sons, New York (2004)

2. Hazan, T., Polak, S., Shashua, A.: Sparse image coding using a 3D non-negative

tensor factorization. In: International Conference of Computer Vision (ICCV), pp.

50–57 (2005)

3. Heiler, M., Schnoerr, C.: Controlling sparseness in non-negative tensor factoriza-

tion. In: Leonardis, A., Bischof, H., Pinz, A. (eds.) ECCV 2006. LNCS, vol. 3951,

pp. 56–67. Springer, Heidelberg (2006)

4. Miwakeichi, F., Martnez-Montes, E., Valds-Sosa, P., Nishiyama, N., Mizuhara, H.,

Yamaguchi, Y.: Decomposing EEG data into space

−time−frequency components

using parallel factor analysi. NeuroImage 22, 1035–1045 (2004)

5. Mørup, M., Hansen, L.K., Herrmann, C.S., Parnas, J., Arnfred, S.M.: Parallel

factor analysis as an exploratory tool for wavelet transformed event-related EEG.

NeuroImage 29, 938–947 (2006)

6. Lee, D.D., Seung, H.S.: Learning the parts of objects by nonnegative matrix fac-

torization. Nature 401, 788–791 (1999)

7. Cichocki, A., Amari, S.: Adaptive Blind Signal and Image Processing (New revised

and improved edition). John Wiley, New York (2003)

8. Dhillon, I., Sra, S.: Generalized nonnegative matrix approximations with Bregman

divergences. In: Neural Information Proc. Systems, Vancouver, Canada, pp. 283–

290 (2005)

9. Berry, M., Browne, M., Langville, A., Pauca, P., Plemmons, R.: Algorithms and ap-

plications for approximate nonnegative matrix factorization. Computational Statis-

tics and Data Analysis (in press, 2006)

10. Bobin, J., Starck, J.L., Fadili, J., Moudden, Y., Donoho, D.L.: Morphological com-

ponent analysis: An adaptive thresholding strategy. IEEE Transactions on Image

Processing (in print, 2007)

11. Elad, M.: Why simple shrinkage is still relevant for redundant representations?

IEEE Trans. On Information Theory 52, 5559–5569 (2006)

12. Kovac, A.: Smooth functions and local extreme values. Computational Statistics

and Data Analysis 51, 5155–5171 (2007)

13. Tao, T., Vidakovic, B.: Almost everywhere behavior of general wavelet shrinkage

operators. Applied and Computational Harmonic Analysis 9, 72–82 (2000)

14. Daubechies, I., Defrrise, M., Mol, C.D.: An iterative thresholding algorithm for lin-

ear inverse problems with a sparsity constraint. Pure and Applied Mathematics 57,

1413–1457 (2004)

15. Cichocki, A., Zdunek, R.: Multilayer nonnegative matrix factorization. Electronics

Letters 42, 947–948 (2006)

16. Cichocki, A., Zdunek, R.: NTFLAB for Signal Processing. Technical report, Labo-

ratory for Advanced Brain Signal Processing, BSI, RIKEN, Saitama, Japan (2006)

Appearance Models for Medical Volumes with

Few Samples by Generalized 3D-PCA

Rui Xu

1

and Yen-Wei Chen

2,3

1

Graduate School of Engineering and Science, Ritsumeikan University, Japan

[e-mail]

2

College of Information Science and Engineering, Ritsumeikan University, Japan

[e-mail]

3

School of Electronic and Information Engineering,

Dalian University of Technology, China

Abstract. Appearance models is important for the task of medical im-

age analysis, such as segmentation. Principal component analysis (PCA)

is an efficient method to build the appearance models; however the 3D

medical volumes should be first unfolded to form the 1D long vectors

before the PCA is used. For large medical volumes, such a unfolding

preprocessing causes two problems. One is the huge burden of comput-

ing cost and the other is bad performance on generalization. A method

named as generalized 3D-PCA is proposed to build the appearance mod-

els for medical volumes in this paper. In our method, the volumes are

directly treated as the third order tensor in the building of the model

without the unfolding preprocessing. The output of our method is three

matrices whose columns are formed by the orthogonal bases in the three

mode subspaces. With the help of these matrices, the bases in the third

order tensor space can be constructed. According to these properties,

our method is not suffered from the two problems of the PCA-based

methods. Eighteen 256

× 256 × 26 MR brain volumes are used in the

experiments of building appearance models. The leave-one-out testing

shows that our method has good performance in building the appear-

ance models for medical volumes even when few samples are used for

training.

1

Introduction

In recent years, a lot of research has been focused on how to extract the sta-

tistical information from 3D medical volumes. Inspired from 2D active shape

models (ASMs) [1], 3D ASMs[2][3] has been proposed to extract the statistical

shape information from 3D medical volumes. The 3D ASMs is a kind of im-

portant method for segmentation but this method does not consider about the

This work was supported in part by the Grand-in Aid for Scientific Research from

the Japanese Ministry for Education, Science, Culture and Sports under the Grand

No. 19500161, and the Strategic Information and Communications R&D Promotion

Programme (SCOPE) under the Grand No. 051307017.

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

c Springer-Verlag Berlin Heidelberg 2008

822

R. Xu and Y.-W. Chen

information of the voxel’s intensity. So, the authors in [4] extend the works of

2D active appearance model (AAMs) [5] and propose the 3DAAMs to combine

the shape models and appearance models together for 3D medical data. In the

3DAAMs, the appearance models for the volumes is built by the method based

on the principal component analysis (PCA).

PCA is a widely-used method to extract statistical information and success-

fully applied into many fields, such as image representation and face recognition.

When the PCA-based methods are applied to build the appearance models for

medical volumes, it is required to first unfold the 3D volumes to 1D long vec-

tors in order to make the data to be processed by PCA. There are two types of

ways to implement the PCA. In the first way, the bases of the unfolding vector

space are obtained by calculating eigenvectors directly from the covariance ma-

trix Cov, where Cov = A

· A

T

and the columns of the matrix A are formed by

the unfolded vectors. The dimension of the covariance matrix is depended on the

dimension of the unfolded vectors. If the 3D volumes have large dimensions, the

dimensions of covariance matrix are also very large. So we encounter the problem

of huge burden of computing cost when we calculate the eigenvectors. The other

way to implement PCA is similar to the one used in the work of eigenfaces [6].

Instead of calculating the the spaces directly from Cov, the eigenvectors u

i

is

first calculated from the matrix Cov = A

T

· A. Then the bases of the unfold-

ing vector space v

i

can be calculated by v

i

= A

· u

i

. For the second way, the

dimension of Cov is depended on the number of the training samples which is

much fewer than the dimension of the unfolded vectors. Therefore, the second

way to implement PCA does not suffer from the problem of huge computing

cost. However, only M

− 1 bases are available in the second way, where M is

the number of training samples. In the medical field, the typical dimension of

unfolding vectors is several millions; however the typical number of the medical

data is only several tens. So the obtained bases are not enough to represent

an untrained vector in the huge vector space. Therefore, the second way has

the problem of bad performance on generalization when the large 3D volumes

are used. The generalization problem may be overcome by increasing the train-

ing samples; however since the dimension of the unfolding vectors is huge, it is

difficult to collect enough medical data for the training. Additionally, with the

development of medical imaging techniques, the medical volumes become larger

and larger. The conflict between the large dimensions of medical volumes and

small number of training samples will become more prominent. Therefore, it is

desirable to research on how to build the appearance model for medical volumes

when the training samples are relative few.

In this paper, a method called generalized 3D-PCA (G3D-PCA) is proposed

to build the appearance models for the medical volumes. In our method, the

3D volumes are treated as the third order tensors in stead of unfolding them to

be the long 1D vectors. The bases of the tensor space are constructed by three

matrices whose columns are set by the leading eigenvectors calculated from the

covariance matrices on the three mode subspaces. Since the dimensions of these

covariance matrices are only depended by the dimensions of the corresponding

Appearance Models for Medical Volumes with Few Samples

823

mode subspaces, our method is able to overcome the problem of huge computing

cost. Additionally, the base number is not limited by the training sample number

in our method, so our method can also overcome the generalization problem.

The proposed method is based on the multilinear algebra [7] [8] which is the

mathematical tool for higher order tensors. Recently, the multilinear algebra

based method is applied to field of computer vision and pattern recognition and

achieves good results. The authors in [9] propose tensorfaces to build statistical

model of faces in different conditions, such as illumination and viewpoint. A

method called concurrent subspaces analysis (CSA) [10] is proposed for video

sequence reconstruction and digital number recognition. However, we does not

find the research to apply the multilinear algebra based method in the field of

medical image analysis.

2

Background Knowledge of Multilinear Algebra

We introduce some background knowledge about multilinear algebra in this sec-

tion. For more detailed knowledge of multilinear algebra, the readers can refer to

[7] [8]. The notations of the multilinear algebra in this paper are listed as follows.

Scalers are denoted by italic-shape letters, i.e. (a, b, ...) or (A, B, ...). Bold lower

case letters, i.e. (a, b, ...), are used to represent vectors. Matrices are denoted by

bold upper case letters, i.e. (A, B, ...); and higher-order tensors are denoted by

calligraphic upper case letters, i.e. (

A, B, ...).

An N -th order tensor

A is defined as a multiway array with N indices, where

A ∈ R

I

1

×I

2

×...×I

N

. Elements of the tensor

A are denoted as a

i

1

...i

n

...i

N

, where

1

i

n

I

n

. The space of the N -th order tensor is comprised by the N mode

subspaces. In the viewpoint of tensor, scales, vectors and matrices can be seen

as zeroth-order, first order and second order tensors respectively.

Varying the n-th index i

n

and keeping the other indices fixed, we can ob-

tain the ”mode-n vectors” of the tensor

A. The ”mode-n matrix” A

(n)

can be

formed by arranging all the mode-n vectors sequentially as its column vectors,

A

(n)

∈ R

I

n

×(I

1

·...I

n

−1

·I

n+1

·...I

N

)

. The procedure of forming the mode-n matrix

is called unfolding of a tensor. Fig. 1 gives the examples to show how to un-

I

1

I

2

I

3

I

1

I

2

I

2

I

2

I

3