A New Learning Algorithm for Blind
Signal Separation
S. Amari

3

University of Tokyo
Bunkyo-ku, Tokyo 113, JAPANamari@sat.t.u-tokyo.ac.jp

A. Cichocki

Lab. for Arti cial Brain Systems
FRP, RIKEN
Wako-Shi, Saitama, 351-01, JAPAN
cia@kamo.riken.go.jp

H. H. Yang

Lab. for Information Representation
FRP, RIKEN
Wako-Shi, Saitama, 351-01, JAPANhhy@koala.riken.go.jp

Abstract
A new on-line learning algorithm which minimizes a statistical dependency among outputs is derived for blind separation of mixed
signals. The dependency is measured by the average mutual information (MI) of the outputs. The source signals and the mixing
matrix are unknown except for the number of the sources. The
Gram-Charlier expansion instead of the Edgeworth expansion isused in evaluating the MI. The natural gradient approach is used
to minimize the MI. A novel activation function is proposed for theon-line learning algorithm which has an equivariant property andis easily implemented on a neural network like model. The validityof the new learning algorithm are veri ed by computer simulations.

1

INTRODUCTION

The problem of blind signal separation arises in many areas such as speech recognition, data communication, sensor signal processing, and medical science. Several
neural network algorithms [3, 5, 7] have been proposed for solving this problem.
The performance of these algorithms is usually a ected by the selection of the activation functions for the formal neurons in the networks. However, all activation
3

Lab. for Information Representation, FRP, RIKEN, Wako-shi, Saitama, JAPANfunctions attempted are monotonic and the selections of the activation functions
are ad hoc. How should the activation function be determined to minimize the MI?
Is it necessary to use monotonic activation functions for blind signal separation? In
this paper, we shall answer these questions and give an on-line learning algorithm
which uses a non-monotonic activation function selected by the independent component analysis (ICA) [7]. Moreover, we shall show a rigorous way to derive thelearning algorithm which has the equivariant property, i.e., the performance of thealgorithm is independent of the scaling parameters in the noiseless case.

2

PROBLEM

Let us consider unknown source signals si (t); i = 1; 1 1 1 ; n which are mutually independent. It is assumed that the sources si (t) are stationary processes and eachsource has moments of any order with a zero mean. The model for the sensor output
isx(t) = As(t)where A 2 Rn2n is an unknown non-singular mixing matrix, s(t) =[s1 (t); 1 1 1 ; sn (t)]T and x(t) = [x1 (t); 1 1 1 ; xn (t)]T .
Without knowing the source signals and the mixing matrix, we want to recover theoriginal signals from the observations x(t) by the following linear transform:y(t) = W x(t)where y(t) = [y 1 (t); 1 1 1 ; y n (t)]T and W 2 Rn2n is a de-mixing matrix.
It is impossible to obtain the original sources si (t) because they are not identi ablein the statistical sense. However, except for a permutation of indices, it is possible
to obtain ci si (t) where the constants ci are inde nite nonzero scalar factors. The
source signals are identi able in this sense. So our goal is to nd the matrix W suchthat [y 1 ; 1 1 1 ; y n ] coincides with a permutation of [s1 ; 1 1 1 ; sn ] except for the scalar
factors. The solution W is the matrix which nds all independent components inthe outputs. An on-line learning algorithm for W is needed which performs the
ICA. It is possible to nd such a learning algorithm which minimizes the dependency
among the outputs. The algorithm in [6] is based on the Edgeworth expansion[8] forevaluating the marginal negentropy. Both the Gram-Charlier expansion[8] and the
Edgeworth expansion[8] can be used to approximate probability density functions.
We shall use the Gram-Charlier expansion instead of the Edgeworth expansion forevaluating the marginal entropy. We shall explain the reason in section 3.

3

INDEPENDENCE OF SIGNALS

The mathematical framework for the ICA is formulated in [6]. The basic idea of the
ICA is to minimize the dependency among the output components. The dependency
is measured by the Kullback-Leibler divergence between the joint and the productof the marginal distributions of the outputs:

D(W ) =

Z

p(y)
a dy
a=1 pa (y )p(y) log Qn(1)where pa (y a ) is the marginal probability density function (pdf). Note the KullbackLeibler divergence has some invariant properties from the di erential-geometricalpoint of view[1].

It is easy to relate the Kullback-Leibler divergence D(W ) to the average MI of y:D(W ) = 0H (y) +
where

R

n
X
H (ya )
a=1

(2)

H (y) = 0 Rp(y) log p(y)dy,
H (ya ) = 0 pa (ya ) log pa (ya )dya is the marginal entropy.

The minimization of the Kullback-Leibler divergence leads to an ICA algorithm forestimating W in [6] where the Edgeworth expansion is used to evaluate the negentropy. We use the truncated Gram-Charlier expansion to evaluate the KullbackLeibler divergence. The Edgeworth expansion has some advantages over the GramCharlier expansion only for some special distributions. In the case of the Gamma
distribution or the distribution of a random variable which is the sum of iid randomvariables, the coecients of the Edgeworth expansion decrease uniformly. However,there is no such advantage for the mixed output ya in general cases.
To calculate each H (ya ) in (2), we shall apply the Gram-Charlier expansion toapproximate the pdf pa (y a ). Since E [y] = E [WAs] = 0, we have E [ya ] = 0. To
simplify the calculations for the entropy H (y a ) to be carried out later, we assume
ma2 = 1. We use the following truncated Gram-Charlier expansion to approximate
the pdf pa (y a ):pa (ya )  (ya )f1 +

a3

a4

H (ya )g(3)
4! 4
where a3 = ma3 , a4 = ma4 0 3, mak = E [(y a )k ] is the k-th order moment of y a ,y2(y ) = p12 e0 2 , and Hk (y) are Chebyshev-Hermite polynomials de ned by theidentitydk (y)(01)kk = Hk (y ) (y ):3!

H3 (ya ) +dy

We prefer the Gram-Charlier expansion to the Edgeworth expansion because theformer clearly shows how a3 and a4 a ect the approximation of the pdf. The lastterm in (3) characterizes non-Gaussian distributions. To apply (3) to calculate
H (ya ), we need the following integrals:
0 (y)H2 (y) log (y)dy = 41
(4)

Z

Z

(y )(H2 (y))2 H4 (y )dy = 24:(5)These integrals can be obtained easily from the following results for the momentsof a Gaussian random variable N(0,1):Z

y2k+1 (y)dy = 0;Zy2k (y)dy = 1 1 3 1 1 1 (2k 0 1):

By using the expansionlog(1 + y )  y 0

y2

(6)

+ O(y 3 )2and taking account of the orthogonality relations of the Chebyshev-Hermite polynomials and (4)-(5), the entropy H (y a ) is expanded as1(a )2 (a )2 5
1
H (ya )  log(2e) 0 3 0 4 + (a3 )2 a4 + (a4 )3 :(7)2
2 1 3! 2 1 4! 8
16

It is easy to calculateZ1
log(2e):2
From y = W x, we have H (y) = H (x) + log jdet(W )j: Applying (7) and the above
expressions to (2), we have0

D(W) (y ) log (y )dy =0H (x) 0 log jdet(W)j + n2 log(2e) 00 58 (a3 )2 a4 0 161 (a4 )3 ]:

4

n
X
(a )

[

a=1

3

2

2 1 3!

+(a4 )2
2 1 4!
(8)

A NEW LEARNING ALGORITHM

To obtain the gradient descent algorithm to update W recursively, we need to@D where wa is the (a,k) element of W in the a-th row and k-th column.
calculate @w
akk

W . It is not dicult to derive the followings:
jdeta W j = cof (wka ) = (W 0T )ak@wk
det(W )

Let cof (wka ) be the cofactor of wka in
@ log

(

@a3
@wka =@a4@wka =)

3E [(ya )2 xk ]4E [(ya )3 xk ]where (W 0T )ak denotes the (a,k) element of (W T )01 . From (8), we obtain
@D
@wka

where

 0(W 0T )ak + f (a3 ; a4 )E [(ya )2 xk ] + g(a3 ; a4 )E [(ya )3 xk ](9)f (y; z) = 0 21 y + 154 yz, g(y; z ) = 0 16 z + 25 y2 + 43 z 2 .

From (9), we obtain the gradient descent algorithm to update W recursively:
dwka
dt =

=

@D
0(t) @wak0
(t)f(W T )ak 0 f (a3 ; a4 )E [(y a )2 xk ] 0 g(a3 ; a4 )E [(ya )3 xk ]g (10)where (t) is a learning rate function. Replacing the expectation values in (10) by
their instantaneous values, we have the stochastic gradient descent algorithm:
dwka
dt =

(t)f(W 0T )ak 0 f (a3 ; a4 )(ya )2 xk 0 g(a3 ; a4 )(y a )3 xk g:(11)

We need to use the following adaptive algorithm to compute a3 and a4 in (11):da3= 0(t)(a3 0 (y a )3 )dtda4
= 0(t)(a4 0 (y a )4 + 3)
dt

(12)

where (t) is another learning rate function.
The performance of the algorithm (11) relies on the estimation of the third andfourth order cumulants performed by the algorithm (12). Replacing the moments

of the random variables in (11) by their instantaneous values, we obtain the following
algorithm which is a direct but coarse implementation of (11):dwka=  (t)f(W 0T )ak 0 f (ya )xk g
dt
where the activation function f (y ) is de ned by(13)

25 9 14 7 47 5 29 33y 0 y 0 y + y:
(14)
4
4
3
4
4
Note the activation function f (y ) is an odd function, not a monotonic function.
The equation (13) can be written in a matrix form:f (y) = y11 +dW=  (t)fW 0T 0 f (y )xT g:
dt

(15)

This equation can be further simpli ed as following by substituting xT W T = yT :dW
= (t)fI 0 f (y )y T gW 0T(16)dt
where f (y ) = (f (y 1 ); 1 1 1 ; f (yn ))T . The above equation is based on the gradientdescent algorithm (10) with the following matrix form:@DdW= 0 (t)
:
dt
@W

(17)

From information geometry perspective[1], since the mixing matrix A is nonsingular we had better replace the above algorithm by the following natural gradientdescent algorithm:dW@D T= 0 (t)
W W:
(18)

dt

@W

D to (18), we obtain the
Applying the previous approximation of the gradient @@W
following algorithm:dW=  (t)fI 0 f (y )y T gW
(19)dt

which has the same \equivariant" property as the algorithms developed in [4, 5].
Although the on-line learning algorithms (16) and (19) look similar to those in[3, 7] and [5] respectively, the selection of the activation function in this paper isrational, not ad hoc. The activation function (14) is determined by the ICA. It isa non-monotonic activation function di erent from those used in [3, 5, 7].
There is a simple way to justify the stability of the algorithm (19). Let Vec(1)denote an operator on a matrix which cascades the columns of the matrix from theleft to the right and forms a column vector. Note this operator has the following
property:
Vec(ABC ) = (C T A)Vec(B ):(20)
Both the gradient descent algorithm and the natural gradient descent algorithm are
special cases of the following general gradient descent algorithm:
dVec(W )@D= 0 (t)P(21)
dt@ Vec(W )where P is a symmetric and positive de nite matrix. It is trivial that (21) becomes(17) when P = I . When P = W T W I , applying (20) to (21), we obtain@D@D T
dVec(W )= 0 (t)(W T W I )= 0 (t)Vec(
W W)
dt@ Vec(W )@W

and this equation implies (18). So the natural gradient descent algorithm updates

W (t) in the direction of decreasing the dependency D(W ). The information geom-etry theory[1] explains why the natural gradient descent algorithm should be used
to minimize the MI.
Another on-line learning algorithm for blind separation using recurrent network was
proposed in [2]. For this algorithm, the activation function (14) also works well.
In practice, other activation functions such as those proposed in [2]-[6] may also be
used in (19). However, the performance of the algorithm for such functions usuallydepends on the distributions of the sources. The activation function (14) works for
relatively general cases in which the pdf of each source can be approximated by thetruncated Gram-Charlier expansion.

5

SIMULATION

In order to check the validity and performance of the new on-line learning algorithm(19), we simulate it on the computer using synthetic source signals and a random
mixing matrix. The extensive computer simulations have fully con rmed the theory
and the validity of the algorithm (19). Due to the limit of space we present here
only one illustrative example.
Example:Assume that the following three unknown sources are mixed by a random mixing
matrix A:[s1 (t); s2 (t); s3 (t)] = [n(t); 0:1sin(400t)cos(30t); 0:01sign[sin(500t + 9cos(40t))]where n(t) is a noise source uniformly distributed in the range [01; +1], and s2 (t)and s3 (t) are two deterministic source signals. The elements of the mixing matrix
A are randomly chosen in [01; +1]. The learning rate is exponentially decreasing
to zero as  (t) = 250exp(05t).
A simulation result is shown in Figure 1. The rst three signals denoted by X1,
X2 and X3 represent mixing (sensor) signals: x1 (t), x2 (t) and x3 (t). The last
three signals denoted by O1, O2 and O3 represent the output signals: y1 (t), y2 (t),and y 3 (t). By using the proposed learning algorithm, the neural network is able
to extract the deterministic signals from the observations after approximately 500
milliseconds.
The performance index E1 is de ned by
n
n X
n
n X
X
X
j
pij j
E = (
0 1) + ( maxjpijjpj j 0 1)
maxk jpik j
k kj
1

i=1 j =1

where P = (pij ) = WA.

6

j =1 i=1

CONCLUSION

The major contribution of this paper the rigorous derivation of the e ective blindseparation algorithm with equivariant property based on the minimization of the
MI of the outputs. The ICA is a general principle to design algorithms for blind
signal separation. The most diculties in applying this principle are to evaluate
the MI of the outputs and to nd a working algorithm which decreases the MI.
Di erent from the work in [6], we use the Gram-Charlier expansion instead of the
Edgeworth expansion to calculate the marginal entropy in evaluating the MI. Usingthe natural gradient method to minimize the MI, we have found an on-line learningalgorithm to nd a de-mixing matrix. The algorithm has equivariant property andcan be easily implemented on a neural network like model. Our approach provides
a rational selection of the activation function for the formal neurons in the network.
The algorithm has been simulated for separating unknown source signals mixed bya random mixing matrix. Our theory and the validity of the new learning algorithm
are veri ed by the simulations.

2 X1
0
-2
2 X2
0
-2
1 X3
0
-1
2 O1
0
-2
1 O2
0
-1
-2
2 O3
1
0
-1
0.6 E1
0.4
0.2
0

0.1

0.2

0.3

0.4

t

Figure 1: The mixed and separated signals, and the performance index
Acknowledgment

We would like to thank Dr. Xiao Yan SU for the proof-reading of the manuscript.

References[1] S.-I. Amari. Di erential-Geometrical Methods in Statistics, Lecture Notes in
Statistics vol.28. Springer, 1985.
[2] S. Amari, A. Cichocki, and H. H. Yang. Recurrent neural networks for blind separation of sources. In Proceedings 1995 International Symposium on Nonlinear
Theory and Applications, volume I, pages 37{42, December 1995.
[3] A. J. Bell and T. J. Sejnowski. An information-maximisation approach to blindseparation and blind deconvolution. Neural Computation, 7:1129{1159, 1995.
[4] J.-F. Cardoso and Beate Laheld. Equivariant adaptive source separation. To
appear in IEEE Trans. on Signal Processing, 1996.
[5] A. Cichocki, R. Unbehauen, L. Moszczynski, and E. Rummert. A new on-line
adaptive learning algorithm for blind separation of source signals. In ISANN94,
pages 406{411, Taiwan, December 1994.
[6] P. Comon. Independent component analysis, a new concept? Signal Processing,
36:287{314, 1994.
[7] C. Jutten and J. Herault. Blind separation of sources, part i: An adaptivealgorithm based on neuromimetic architecture. Signal Processing, 24:1{10, 1991.
[8] A. Stuart and J. K. Ord. Kendall's Advanced Theory of Statistics. Edward
Arnold, 1994.