Independent Component Analysis:
Algorithms and Applications
Aapo Hyvärinen and Erkki Oja
Neural Networks Research Centre
Helsinki University of Technology
P.O. Box 5400, FIN-02015 HUT, Finland
Neural Networks, 13(4-5):411-430, 2000
Abstract
A fundamental problem in neural network research, as well as in many other disciplines, is finding a suitable
representation of multivariate data, i.e. random vectors. For reasons of computational and conceptual simplicity,
the representation is often sought as a linear transformation of the original data. In other words, each componentof the representation is a linear combination of the original variables. Well-known linear transformation methods
include principal component analysis, factor analysis, and projection pursuit. Independent component analysis(ICA) is a recently developed method in which the goal is to find a linear representation of nongaussian data sothat the components are statistically independent, or as independent as possible. Such a representation seems tocapture the essential structure of the data in many applications, including feature extraction and signal separation.
In this paper, we present the basic theory and applications of ICA, and our recent work on the subject.

Keywords: Independent component analysis, projection pursuit, blind signal separation, source separation, factor
analysis, representation

1 Motivation
Imagine that you are in a room where two people are speaking simultaneously. You have two microphones, which
you hold in different locations. The microphones give you two recorded time signals, which we could denote byx1 (t) and x2 (t), with x1 and x2 the amplitudes, and t the time index. Each of these recorded signals is a weighted
sum of the speech signals emitted by the two speakers, which we denote by s1 (t) and s2 (t). We could express thisas a linear equation:
x1 (t) = a11 s1 + a12s2
x2 (t) = a21 s1 + a22s2

(1)
(2)

where a11 , a12 , a21 , and a22 are some parameters that depend on the distances of the microphones from the speakers.
It would be very useful if you could now estimate the two original speech signals s1 (t) and s2 (t), using only therecorded signals x1 (t) and x2 (t). This is called the cocktail-party problem. For the time being, we omit any timedelays or other extra factors from our simplified mixing model.
As an illustration, consider the waveforms in Fig. 1 and Fig. 2. These are, of course, not realistic speech signals,
but suffice for this illustration. The original speech signals could look something like those in Fig. 1 and the mixed
signals could look like those in Fig. 2. The problem is to recover the data in Fig. 1 using only the data in Fig. 2.
Actually, if we knew the parameters ai j , we could solve the linear equation in (1) by classical methods. The
point is, however, that if you don’t know the ai j , the problem is considerably more difficult.
One approach to solving this problem would be to use some information on the statistical properties of thesignals si (t) to estimate the aii . Actually, and perhaps surprisingly, it turns out that it is enough to assume that s1 (t)and s2 (t), at each time instant t, are statistically independent. This is not an unrealistic assumption in many cases,
1and it need not be exactly true in practice. The recently developed technique of Independent Component Analysis,or ICA, can be used to estimate the ai j based on the information of their independence, which allows us to separate
the two original source signals s1 (t) and s2 (t) from their mixtures x1 (t) and x2 (t). Fig. 3 gives the two signals
estimated by the ICA method. As can be seen, these are very close to the original source signals (their signs arereversed, but this has no significance.)
Independent component analysis was originally developed to deal with problems that are closely related to the
cocktail-party problem. Since the recent increase of interest in ICA, it has become clear that this principle has a
lot of other interesting applications as well.
Consider, for example, electrical recordings of brain activity as given by an electroencephalogram (EEG). The
EEG data consists of recordings of electrical potentials in many different locations on the scalp. These potentials
are presumably generated by mixing some underlying components of brain activity. This situation is quite similar
to the cocktail-party problem: we would like to find the original components of brain activity, but we can only
observe mixtures of the components. ICA can reveal interesting information on brain activity by giving access to
its independent components.
Another, very different application of ICA is on feature extraction. A fundamental problem in digital signal
processing is to find suitable representations for image, audio or other kind of data for tasks like compression anddenoising. Data representations are often based on (discrete) linear transformations. Standard linear transformations widely used in image processing are the Fourier, Haar, cosine transforms etc. Each of them has its own
favorable properties (Gonzales and Wintz, 1987).
It would be most useful to estimate the linear transformation from the data itself, in which case the transform
could be ideally adapted to the kind of data that is being processed. Figure 4 shows the basis functions obtained by
ICA from patches of natural images. Each image window in the set of training images would be a superpositionof these windows so that the coefficient in the superposition are independent. Feature extraction by ICA will be
explained in more detail later on.
All of the applications described above can actually be formulated in a unified mathematical framework, thatof ICA. This is a very general-purpose method of signal processing and data analysis.
In this review, we cover the definition and underlying principles of ICA in Sections 2 and 3. Then, startingfrom Section 4, the ICA problem is solved on the basis of minimizing or maximizing certain conrast functions;this transforms the ICA problem to a numerical optimization problem. Many contrast functions are given andthe relations between them are clarified. Section 5 covers a useful preprocessing that greatly helps solving the
ICA problem, and Section 6 reviews one of the most efficient practical learning rules for solving the problem, the
FastICA algorithm. Then, in Section 7, typical applications of ICA are covered: removing artefacts from brain
signal recordings, finding hidden factors in financial time series, and reducing noise in natural images. Section 8
concludes the text.

2 Independent Component Analysis
2.1 Definition of ICA
To rigorously define ICA (Jutten and Hérault, 1991; Comon, 1994), we can use a statistical “latent variables”
model. Assume that we observe n linear mixtures x1 , ..., xn of n independent components
x j = a j1 s1 + a j2 s2 + ... + a jnsn , for all j.

(3)

We have now dropped the time index t; in the ICA model, we assume that each mixture x j as well as eachindependent component sk is a random variable, instead of a proper time signal. The observed values x j (t), e.g.,
the microphone signals in the cocktail party problem, are then a sample of this random variable. Without loss ofgenerality, we can assume that both the mixture variables and the independent components have zero mean: If thisis not true, then the observable variables xi can always be centered by subtracting the sample mean, which makes
the model zero-mean.
It is convenient to use vector-matrix notation instead of the sums like in the previous equation. Let us denote by
x the random vector whose elements are the mixtures x1 , ..., xn , and likewise by s the random vector with elements2

s1 , ..., sn . Let us denote by A the matrix with elements ai j . Generally, bold lower case letters indicate vectors andbold upper-case letters denote matrices. All vectors are understood as column vectors; thus xT , or the transpose ofx, is a row vector. Using this vector-matrix notation, the above mixing model is written as
x = As.

(4)

Sometimes we need the columns of matrix A; denoting them by a j the model can also be written asn

x = ∑ ai si .

(5)

i=1

The statistical model in Eq. 4 is called independent component analysis, or ICA model. The ICA model is a
generative model, which means that it describes how the observed data are generated by a process of mixing thecomponents si . The independent components are latent variables, meaning that they cannot be directly observed.
Also the mixing matrix is assumed to be unknown. All we observe is the random vector x, and we must estimate
both A and s using it. This must be done under as general assumptions as possible.
The starting point for ICA is the very simple assumption that the components si are statistically independent.
Statistical independence will be rigorously defined in Section 3. It will be seen below that we must also assume thatthe independent component must have nongaussian distributions. However, in the basic model we do not assume
these distributions known (if they are known, the problem is considerably simplified.) For simplicity, we are also
assuming that the unknown mixing matrix is square, but this assumption can be sometimes relaxed, as explainedin Section 4.5. Then, after estimating the matrix A, we can compute its inverse, say W, and obtain the independent
component simply by:
s = Wx.
(6)
ICA is very closely related to the method called blind source separation (BSS) or blind signal separation. A“source” means here an original signal, i.e. independent component, like the speaker in a cocktail party problem.
“Blind” means that we no very little, if anything, on the mixing matrix, and make little assumptions on the source
signals. ICA is one method, perhaps the most widely used, for performing blind source separation.
In many applications, it would be more realistic to assume that there is some noise in the measurements (seee.g. (Hyvärinen, 1998a; Hyvärinen, 1999c)), which would mean adding a noise term in the model. For simplicity,
we omit any noise terms, since the estimation of the noise-free model is difficult enough in itself, and seems to be
sufficient for many applications.

2.2 Ambiguities of ICA
In the ICA model in Eq. (4), it is easy to see that the following ambiguities will hold:1. We cannot determine the variances (energies) of the independent components.
The reason is that, both s and A being unknown, any scalar multiplier in one of the sources si could always
be cancelled by dividing the corresponding column ai of A by the same scalar; see eq. (5). As a consequence,
we may quite as well fix the magnitudes of the independent components; as they are random variables, the mostnatural way to do this is to assume that each has unit variance: E{s2i } = 1. Then the matrix A will be adapted in the
ICA solution methods to take into account this restriction. Note that this still leaves the ambiguity of the sign: we
could multiply the an independent component by −1 without affecting the model. This ambiguity is, fortunately,
insignificant in most applications.
2. We cannot determine the order of the independent components.
The reason is that, again both s and A being unknown, we can freely change the order of the terms in the sumin (5), and call any of the independent components the first one. Formally, a permutation matrix P and its inverse
can be substituted in the model to give x = AP−1 Ps. The elements of Ps are the original independent variables s j ,but in another order. The matrix AP−1 is just a new unknown mixing matrix, to be solved by the ICA algorithms.

3

2.3 Illustration of ICA
To illustrate the ICA model in statistical terms, consider two independent components that have the following
uniform distributions:(√1√
if |si | ≤ 3
2
3
p(si ) =(7)0otherwise
The range of values for this uniform distribution were chosen so as to make the mean zero and the variance equal
to one, as was agreed in the previous Section. The joint density of s1 and s2 is then uniform on a square. This
follows from the basic definition that the joint density of two independent variables is just the product of their
marginal densities (see Eq. 10): we need to simply compute the product. The joint density is illustrated in Figure 5
by showing data points randomly drawn from this distribution.
Now let as mix these two independent components. Let us take the following mixing matrix:


2 3
(8)
A0 =
2 1
This gives us two mixed variables, x1 and x2 . It is easily computed that the mixed data has a uniform distribution
on a parallelogram, as shown in Figure 6. Note that the random variables x1 and x2 are not independent any more;an easy way to see this is to consider, whether it is possible to predict the value of one of them, say x2 , from thevalue of the other. Clearly if x1 attains one of its maximum or minimum values, then this completely determines
the value of x2 . They are therefore not independent. (For variables s1 and s2 the situation is different: from Fig. 5
it can be seen that knowing the value of s1 does not in any way help in guessing the value of s2 .)The problem of estimating the data model of ICA is now to estimate the mixing matrix A0 using only information contained in the mixtures x1 and x2 . Actually, you can see from Figure 6 an intuitive way of estimating A:
The edges of the parallelogram are in the directions of the columns of A. This means that we could, in principle,
estimate the ICA model by first estimating the joint density of x1 and x2 , and then locating the edges. So, theproblem seems to have a solution.
In reality, however, this would be a very poor method because it only works with variables that have exactly
uniform distributions. Moreover, it would be computationally quite complicated. What we need is a method thatworks for any distributions of the independent components, and works fast and reliably.
Next we shall consider the exact definition of independence before starting to develop methods for estimationof the ICA model.

3 What is independence?
3.1 Definition and fundamental properties
To define the concept of independence, consider two scalar-valued random variables y1 and y2 . Basically, thevariables y1 and y2 are said to be independent if information on the value of y1 does not give any information on
the value of y2 , and vice versa. Above, we noted that this is the case with the variables s1 , s2 but not with the
mixture variables x1 , x2 .
Technically, independence can be defined by the probability densities. Let us denote by p(y1 , y2 ) the joint
probability density function (pdf) of y1 and y2 . Let us further denote by p1 (y1 ) the marginal pdf of y1 , i.e. the pdf
of y1 when it is considered alone:Z
p1 (y1 ) =p(y1 , y2 )dy2 ,(9)and similarly for y2 . Then we define that y1 and y2 are independent if and only if the joint pdf is factorizable in thefollowing way:p(y1 , y2 ) = p1 (y1 )p2 (y2 ).
(10)This definition extends naturally for any number n of random variables, in which case the joint density must be a
product of n terms.

4

The definition can be used to derive a most important property of independent random variables. Given two
functions, h1 and h2 , we always haveE{h1(y1 )h2 (y2 )} = E{h1 (y1 )}E{h2(y2 )}.

(11)This can be proven as follows:E{h1 (y1 )h2 (y2 )} =
=

Z Z

Z Z

h1 (y1 )h2 (y2 )p(y1 , y2 )dy1 dy2

h1 (y1 )p1 (y1 )h2 (y2 )p2 (y2 )dy1 dy2 =Z

h1 (y1 )p1 (y1 )dy1Zh2 (y2 )p2 (y2 )dy2
= E{h1(y1 )}E{h2 (y2 )}. (12)3.2 Uncorrelated variables are only partly independent
A weaker form of independence is uncorrelatedness. Two random variables y1 and y2 are said to be uncorrelated,
if their covariance is zero:E{y1 y2 } − E{y1}E{y2 } = 0(13)
If the variables are independent, they are uncorrelated, which follows directly from Eq. (11), taking h1 (y1 ) = y1
and h2 (y2 ) = y2 .
On the other hand, uncorrelatedness does not imply independence. For example, assume that (y1 , y2 ) are
discrete valued and follow such a distribution that the pair are with probability 1/4 equal to any of the following
values: (0, 1), (0, −1), (1, 0), (−1, 0). Then y1 and y2 are uncorrelated, as can be simply calculated. On the other
hand,
1
(14)
E{y21 y22 } = 0 6= = E{y21 }E{y22 }.
4
so the condition in Eq. (11) is violated, and the variables cannot be independent.
Since independence implies uncorrelatedness, many ICA methods constrain the estimation procedure so that it
always gives uncorrelated estimates of the independent components. This reduces the number of free parameters,
and simplifies the problem.

3.3 Why Gaussian variables are forbidden
The fundamental restriction in ICA is that the independent components must be nongaussian for ICA to be possible.
To see why gaussian variables make ICA impossible, assume that the mixing matrix is orthogonal and the si
are gaussian. Then x1 and x2 are gaussian, uncorrelated, and of unit variance. Their joint density is given byp(x1 , x2 ) =1x2 + x22
exp(− 1
)
2π
2

(15)

This distribution is illustrated in Fig. 7. The Figure shows that the density is completely symmetric. Therefore, itdoes not contain any information on the directions of the columns of the mixing matrix A. This is why A cannot
be estimated.
More rigorously, one can prove that the distribution of any orthogonal transformation of the gaussian (x1 , x2 )has exactly the same distribution as (x1 , x2 ), and that x1 and x2 are independent. Thus, in the case of gaussianvariables, we can only estimate the ICA model up to an orthogonal transformation. In other words, the matrix A
is not identifiable for gaussian independent components. (Actually, if just one of the independent components is
gaussian, the ICA model can still be estimated.)5

4 Principles of ICA estimation
4.1 “Nongaussian is independent”
Intuitively speaking, the key to estimating the ICA model is nongaussianity. Actually, without nongaussianity theestimation is not possible at all, as mentioned in Sec. 3.3. This is at the same time probably the main reason for
the rather late resurgence of ICA research: In most of classical statistical theory, random variables are assumed tohave gaussian distributions, thus precluding any methods related to ICA.
The Central Limit Theorem, a classical result in probability theory, tells that the distribution of a sum of
independent random variables tends toward a gaussian distribution, under certain conditions. Thus, a sum of two
independent random variables usually has a distribution that is closer to gaussian than any of the two original
random variables.
Let us now assume that the data vector x is distributed according to the ICA data model in Eq. 4, i.e. it
is a mixture of independent components. For simplicity, let us assume in this section that all the independent
components have identical distributions. To estimate one of the independent components, we consider a linear
combination of the xi (see eq. 6); let us denote this by y = wT x = ∑i wi xi , where w is a vector to be determined. If
w were one of the rows of the inverse of A, this linear combination would actually equal one of the independent
components. The question is now: How could we use the Central Limit Theorem to determine w so that it would
equal one of the rows of the inverse of A? In practice, we cannot determine such a w exactly, because we have noknowledge of matrix A, but we can find an estimator that gives a good approximation.
To see how this leads to the basic principle of ICA estimation, let us make a change of variables, definingz = AT w. Then we have y = wT x = wT As = zT s. y is thus a linear combination of si , with weights given by zi .
Since a sum of even two independent random variables is more gaussian than the original variables, zT s is more
gaussian than any of the si and becomes least gaussian when it in fact equals one of the si . In this case, obviously
only one of the elements zi of z is nonzero. (Note that the si were here assumed to have identical distributions.)
Therefore, we could take as w a vector that maximizes the nongaussianity of wT x. Such a vector would
necessarily correspond (in the transformed coordinate system) to a z which has only one nonzero component. Thismeans that wT x = zT s equals one of the independent components!
Maximizing the nongaussianity of wT x thus gives us one of the independent components. In fact, the optimization landscape for nongaussianity in the n-dimensional space of vectors w has 2n local maxima, two for eachindependent component, corresponding to si and −si (recall that the independent components can be estimatedonly up to a multiplicative sign). To find several independent components, we need to find all these local maxima.
This is not difficult, because the different independent components are uncorrelated: We can always constrain the
search to the space that gives estimates uncorrelated with the previous ones. This corresponds to orthogonalization
in a suitably transformed (i.e. whitened) space.
Our approach here is rather heuristic, but it will be seen in the next section and Sec. 4.3 that it has a perfectly
rigorous justification.

4.2 Measures of nongaussianity
To use nongaussianity in ICA estimation, we must have a quantitative measure of nongaussianity of a randomvariable, say y. To simplify things, let us assume that y is centered (zero-mean) and has variance equal to one.
Actually, one of the functions of preprocessing in ICA algorithms, to be covered in Section 5, is to make thissimplification possible.
4.2.1 Kurtosis
The classical measure of nongaussianity is kurtosis or the fourth-order cumulant. The kurtosis of y is classicallydefined bykurt(y) = E{y4 } − 3(E{y2})2(16)
Actually, since we assumed that y is of unit variance, the right-hand side simplifies to E{y4 } − 3. This shows thatkurtosis is simply a normalized version of the fourth moment E{y4 }. For a gaussian y, the fourth moment equals
63(E{y2 })2 . Thus, kurtosis is zero for a gaussian random variable. For most (but not quite all) nongaussian randomvariables, kurtosis is nonzero.
Kurtosis can be both positive or negative. Random variables that have a negative kurtosis are called subgaussian, and those with positive kurtosis are called supergaussian. In statistical literature, the corresponding
expressions platykurtic and leptokurtic are also used. Supergaussian random variables have typically a “spiky”
pdf with heavy tails, i.e. the pdf is relatively large at zero and at large values of the variable, while being small
for intermediate values. A typical example is the Laplace distribution, whose pdf (normalized to unit variance) is
given by√1
(17)p(y) = √ exp( 2|y|)2This pdf is illustrated in Fig. 8. Subgaussian random variables, on the other hand, have typically a “flat” pdf, whichis rather constant near zero, and very small for larger values of the variable. A typical example is the uniform
distibution in eq. (7).
Typically nongaussianity is measured by the absolute value of kurtosis. The square of kurtosis can also be
used. These are zero for a gaussian variable, and greater than zero for most nongaussian random variables. There
are nongaussian random variables that have zero kurtosis, but they can be considered as very rare.
Kurtosis, or rather its absolute value, has been widely used as a measure of nongaussianity in ICA and related
fields. The main reason is its simplicity, both computational and theoretical. Computationally, kurtosis can be
estimated simply by using the fourth moment of the sample data. Theoretical analysis is simplified because of thefollowing linearity property: If x1 and x2 are two independent random variables, it holds

andkurt(x1 + x2 ) = kurt(x1 ) + kurt(x2 )(18)kurt(αx1 ) = α4 kurt(x1 )(19)

where α is a scalar. These properties can be easily proven using the definition.
To illustrate in a simple example what the optimization landscape for kurtosis looks like, and how independent
components could be found by kurtosis minimization or maximization, let us look at a 2-dimensional model x =As. Assume that the independent components s1 , s2 have kurtosis values kurt(s1 ), kurt(s2 ), respectively, bothdifferent from zero. Remember that we assumed that they have unit variances. We seek for one of the independent
components as y = wT x.
Let us again make the transformation z = AT w. Then we have y = wT x = wT As = zT s = z1 s1 + z2 s2 . Now,
based on the additive property of kurtosis, we have kurt(y) = kurt(z1 s1 ) + kurt(z2 s2 ) = z41 kurt(s1 ) + z42 kurt(s2 ).
On the other hand, we made the constraint that the variance of y is equal to 1, based on the same assumption
concerning s1 , s2 . This implies a constraint on z: E{y2 } = z21 + z22 = 1. Geometrically, this means that vector z isconstrained to the unit circle on the 2-dimensional plane. The optimization problem is now: what are the maxima
of the function | kurt(y)| = |z41 kurt(s1 ) + z42 kurt(s2 )| on the unit circle? For simplicity, you may consider that the
kurtosis are of the same sign, in which case it absolute value operators can be omitted. The graph of this functionis the "optimization landscape" for the problem.
It is not hard to show (Delfosse and Loubaton, 1995) that the maxima are at the points when exactly one of theelements of vector z is zero and the other nonzero; because of the unit circle constraint, the nonzero element must
be equal to 1 or -1. But these points are exactly the ones when y equals one of the independent components ±si ,and the problem has been solved.
In practice we would start from some weight vector w, compute the direction in which the kurtosis of y = wT x
is growing most strongly (if kurtosis is positive) or decreasing most strongly (if kurtosis is negative) based on theavailable sample x(1), ..., x(T ) of mixture vector x, and use a gradient method or one of their extensions for finding
a new vector w. The example can be generalized to arbitrary dimensions, showing that kurtosis can theoreticallybe used as an optimization criterion for the ICA problem.
However, kurtosis has also some drawbacks in practice, when its value has to be estimated from a measured
sample. The main problem is that kurtosis can be very sensitive to outliers (Huber, 1985). Its value may depend on
only a few observations in the tails of the distribution, which may be erroneous or irrelevant observations. In otherwords, kurtosis is not a robust measure of nongaussianity.
7

Thus, other measures of nongaussianity might be better than kurtosis in some situations. Below we shallconsider negentropy whose properties are rather opposite to those of kurtosis, and finally introduce approximationsof negentropy that more or less combine the good properties of both measures.
4.2.2 Negentropy
A second very important measure of nongaussianity is given by negentropy. Negentropy is based on the informationtheoretic quantity of (differential) entropy.
Entropy is the basic concept of information theory. The entropy of a random variable can be interpreted as
the degree of information that the observation of the variable gives. The more “random”, i.e. unpredictable andunstructured the variable is, the larger its entropy. More rigorously, entropy is closely related to the coding lengthof the random variable, in fact, under some simplifying assumptions, entropy is the coding length of the randomvariable. For introductions on information theory, see e.g. (Cover and Thomas, 1991; Papoulis, 1991).
Entropy H is defined for a discrete random variable Y as
H(Y ) = − ∑ P(Y = ai ) log P(Y = ai )

(20)

i

where the ai are the possible values of Y . This very well-known definition can be generalized for continuous-valued
random variables and vectors, in which case it is often called differential entropy. The differential entropy H of arandom vector y with density f (y) is defined as (Cover and Thomas, 1991; Papoulis, 1991):
H(y) = −Zf (y) log f (y)dy.

(21)

A fundamental result of information theory is that a gaussian variable has the largest entropy among all
random variables of equal variance. For a proof, see e.g. (Cover and Thomas, 1991; Papoulis, 1991). This meansthat entropy could be used as a measure of nongaussianity. In fact, this shows that the gaussian distribution is
the “most random” or the least structured of all distributions. Entropy is small for distributions that are clearly
concentrated on certain values, i.e., when the variable is clearly clustered, or has a pdf that is very “spiky”.
To obtain a measure of nongaussianity that is zero for a gaussian variable and always nonnegative, one often
uses a slightly modified version of the definition of differential entropy, called negentropy. Negentropy J is defined
as follows
J(y) = H(ygauss ) − H(y)(22)
where ygauss is a Gaussian random variable of the same covariance matrix as y. Due to the above-mentioned
properties, negentropy is always non-negative, and it is zero if and only if y has a Gaussian distribution. Negentropy has the additional interesting property that it is invariant for invertible linear transformations (Comon, 1994;
Hyvärinen, 1999e).
The advantage of using negentropy, or, equivalently, differential entropy, as a measure of nongaussianity is that
it is well justified by statistical theory. In fact, negentropy is in some sense the optimal estimator of nongaussianity,
as far as statistical properties are concerned. The problem in using negentropy is, however, that it is computationallyvery difficult. Estimating negentropy using the definition would require an estimate (possibly nonparametric) of
the pdf. Therefore, simpler approximations of negentropy are very useful, as will be discussed next.
4.2.3 Approximations of negentropy
The estimation of negentropy is difficult, as mentioned above, and therefore this contrast function remains mainly
a theoretical one. In practice, some approximation have to be used. Here we introduce approximations that havevery promising properties, and which will be used in the following to derive an efficient method for ICA.
The classical method of approximating negentropy is using higher-order moments, for example as follows(Jones and Sibson, 1987):11kurt(y)2(23)
J(y) ≈ E{y3 }2 +
12
48

8

The random variable y is assumed to be of zero mean and unit variance. However, the validity of such approximations may be rather limited. In particular, these approximations suffer from the nonrobustness encountered withkurtosis.
To avoid the problems encountered with the preceding approximations of negentropy, new approximations
were developed in (Hyvärinen, 1998b). These approximation were based on the maximum-entropy principle. In
general we obtain the following approximation:p

J(y) ≈ ∑ ki [E{Gi (y)} − E{Gi(ν)}]2 ,

(24)

i=1

where ki are some positive constants, and ν is a Gaussian variable of zero mean and unit variance (i.e., standardized). The variable y is assumed to be of zero mean and unit variance, and the functions Gi are some nonquadratic
functions (Hyvärinen, 1998b). Note that even in cases where this approximation is not very accurate, (24) can beused to construct a measure of nongaussianity that is consistent in the sense that it is always non-negative, andequal to zero if y has a Gaussian distribution.
In the case where we use only one nonquadratic function G, the approximation becomesJ(y) ∝ [E{G(y)} − E{G(ν)}]2(25)

for practically any non-quadratic function G. This is clearly a generalization of the moment-based approximation in(23), if y is symmetric. Indeed, taking G(y) = y4 , one then obtains exactly (23), i.e. a kurtosis-based approximation.
But the point here is that by choosing G wisely, one obtains approximations of negentropy that are muchbetter than the one given by (23). In particular, choosing G that does not grow too fast, one obtains more robust
estimators. The following choices of G have proved very useful:
G1 (u) =1
log cosh a1 u, G2 (u) = − exp(−u2 /2)
a1

(26)

where 1 ≤ a1 ≤ 2 is some suitable constant.
Thus we obtain approximations of negentropy that give a very good compromise between the properties of
the two classical nongaussianity measures given by kurtosis and negentropy. They are conceptually simple, fast
to compute, yet have appealing statistical properties, especially robustness. Therefore, we shall use these contrast
functions in our ICA methods. Since kurtosis can be expressed in this same framework, it can still be used by our
ICA methods. A practical algorithm based on these contrast function will be presented in Section 6.

4.3 Minimization of Mutual Information
Another approach for ICA estimation, inspired by information theory, is minimization of mutual information. We
will explain this approach here, and show that it leads to the same principle of finding most nongaussian directions
as was described above. In particular, this approach gives a rigorous justification for the heuristic principles used
above.
4.3.1 Mutual Information
Using the concept of differential entropy, we define the mutual information I between m (scalar) random variables,
yi , i = 1...m as followsm

I(y1 , y2 , ..., ym ) = ∑ H(yi ) − H(y).

(27)

i=1

Mutual information is a natural measure of the dependence between random variables. In fact, it is equivalent
to the well-known Kullback-Leibler divergence between the joint density f (y) and the product of its marginaldensities; a very natural measure for independence. It is always non-negative, and zero if and only if the variables
are statistically independent. Thus, mutual information takes into account the whole dependence structure of thevariables, and not only the covariance, like PCA and related methods.
9Mutual information can be interpreted by using the interpretation of entropy as code length. The terms H(yi )give the lengths of codes for the yi when these are coded separately, and H(y) gives the code length when y iscoded as a random vector, i.e. all the components are coded in the same code. Mutual information thus shows whatcode length reduction is obtained by coding the whole vector instead of the separate components. In general, better
codes can be obtained by coding the whole vector. However, if the yi are independent, they give no information oneach other, and one could just as well code the variables separately without increasing code length.
An important property of mutual information (Papoulis, 1991; Cover and Thomas, 1991) is that we have for an
invertible linear transformation y = Wx:
I(y1 , y2 , ..., yn ) = ∑ H(yi ) − H(x) − log| det W|.

(28)

i

Now, let us consider what happens if we constrain the yi to be uncorrelated and of unit variance. This means
E{yyT } = WE{xxT }WT = I, which implies
detI = 1 = (det WE{xxT }WT ) = (detW)(det E{xxT })(det WT ),(29)

and this implies that det W must be constant. Moreover, for yi of unit variance, entropy and negentropy differ only
by a constant, and the sign. Thus we obtain,
I(y1 , y2 , ..., yn ) = C − ∑ J(yi ).

(30)

i

where C is a constant that does not depend on W. This shows the fundamental relation between negentropy andmutual information.
4.3.2 Defining ICA by Mutual Information
Since mutual information is the natural information-theoretic measure of the independence of random variables,
we could use it as the criterion for finding the ICA transform. In this approach that is an alternative to the model
estimation approach, we define the ICA of a random vector x as an invertible transformation as in (6), where thematrix W is determined so that the mutual information of the transformed components si is minimized.
It is now obvious from (30) that finding an invertible transformation W that minimizes the mutual information isroughly equivalent to finding directions in which the negentropy is maximized. More precisely, it is roughly equivalent to finding 1-D subspaces such that the projections in those subspaces have maximum negentropy. Rigorously,
speaking, (30) shows that ICA estimation by minimization of mutual information is equivalent to maximizing thesum of nongaussianities of the estimates, when the estimates are constrained to be uncorrelated. The constraintof uncorrelatedness is in fact not necessary, but simplifies the computations considerably, as one can then use thesimpler form in (30) instead of the more complicated form in (28).
Thus, we see that the formulation of ICA as minimization of mutual information gives another rigorous justification of our more heuristically introduced idea of finding maximally nongaussian directions.

4.4 Maximum Likelihood Estimation
4.4.1 The likelihood
A very popular approach for estimating the ICA model is maximum likelihood estimation, which is closely connected to the infomax principle. Here we discuss this approach, and show that it is essentially equivalent tominimization of mutual information.
It is possible to formulate directly the likelihood in the noise-free ICA model, which was done in (Pham et al.,1992), and then estimate the model by a maximum likelihood method. Denoting by W = (w1 , ..., wn )T the matrix
A−1 , the log-likelihood takes the form (Pham et al., 1992):
Tn

L = ∑ ∑ log fi (wTi x(t)) + T log | det W|
t=1 i=110(31)where the fi are the density functions of the si (here assumed to be known), and the x(t),t = 1, ..., T are therealizations of x. The term log | detW| in the likelihood comes from the classic rule for (linearly) transforming
random variables and their densities (Papoulis, 1991): In general, for any random vector x with density px and forany matrix W, the density of y = Wx is given by px (Wx)| detW|.
4.4.2 The Infomax Principle
Another related contrast function was derived from a neural network viewpoint in (Bell and Sejnowski, 1995;
Nadal and Parga, 1994). This was based on maximizing the output entropy (or information flow) of a neuralnetwork with non-linear outputs. Assume that x is the input to the neural network whose outputs are of the form
φi (wTi x), where the φi are some non-linear scalar functions, and the wi are the weight vectors of the neurons. One
then wants to maximize the entropy of the outputs:
L2 = H(φ1 (wT1 x), ..., φn (wTn x)).

(32)

If the φi are well chosen, this framework also enables the estimation of the ICA model. Indeed, several authors, e.g.,(Cardoso, 1997; Pearlmutter and Parra, 1997), proved the surprising result that the principle of network entropy
maximization, or “infomax”, is equivalent to maximum likelihood estimation. This equivalence requires that the
non-linearities φi used in the neural network are chosen as the cumulative distribution functions corresponding to
the densities fi , i.e., φ0i (.) = fi (.).
4.4.3 Connection to mutual information
To see the connection between likelihood and mutual information, consider the expectation of the log-likelihood:n1
E{L} = ∑ E{log fi (wTi x)} + log| det W|.
T
i=1

(33)

Actually, if the fi were equal to the actual distributions of wTi x, the first term would be equal to − ∑i H(wTi x).
Thus the likelihood would be equal, up to an additive constant, to the negative of mutual information as given in
Eq. (28).
Actually, in practice the connection is even stronger. This is because in practice we don’t know the distributionsof the independent components. A reasonable approach would be to estimate the density of wTi x as part of the ML
estimation method, and use this as an approximation of the density of si . In this case, likelihood and mutual
information are, for all practical purposes, equivalent.
Nevertheless, there is a small difference that may be very important in practice. The problem with maximum
likelihood estimation is that the densities fi must be estimated correctly. They need not be estimated with any
great precision: in fact it is enough to estimate whether they are sub- or supergaussian (Cardoso and Laheld, 1996;
Hyvärinen and Oja, 1998; Lee et al., 1999). In many cases, in fact, we have enough prior knowledge on the
independent components, and we don’t need to estimate their nature from the data. In any case, if the informationon the nature of the independent components is not correct, ML estimation will give completely wrong results.
Some care must be taken with ML estimation, therefore. In contrast, using reasonable measures of nongaussianity,
this problem does not usually arise.

4.5 ICA and Projection Pursuit
It is interesting to note how our approach to ICA makes explicit the connection between ICA and projection
pursuit. Projection pursuit (Friedman and Tukey, 1974; Friedman, 1987; Huber, 1985; Jones and Sibson, 1987) isa technique developed in statistics for finding “interesting” projections of multidimensional data. Such projections
can then be used for optimal visualization of the data, and for such purposes as density estimation and regression.
In basic (1-D) projection pursuit, we try to find directions such that the projections of the data in those directions
have interesting distributions, i.e., display some structure. It has been argued by Huber (Huber, 1985) and by Jonesand Sibson (Jones and Sibson, 1987) that the Gaussian distribution is the least interesting one, and that the most11

interesting directions are those that show the least Gaussian distribution. This is exactly what we do to estimate
the ICA model.
The usefulness of finding such projections can be seen in Fig. 9, where the projection on the projection pursuit
direction, which is horizontal, clearly shows the clustered structure of the data. The projection on the first principal
component (vertical), on the other hand, fails to show this structure.
Thus, in the general formulation, ICA can be considered a variant of projection pursuit. All the nongaussianity
measures and the corresponding ICA algorithms presented here could also be called projection pursuit “indices”
and algorithms. In particular, the projection pursuit allows us to tackle the situation where there are less independent components si than original variables xi is. Assuming that those dimensions of the space that are not spannedby the independent components are filled by gaussian noise, we see that computing the nongaussian projection
pursuit directions, we effectively estimate the independent components. When all the nongaussian directions havebeen found, all the independent components have been estimated. Such a procedure can be interpreted as a hybridof projection pursuit and ICA.
However, it should be noted that in the formulation of projection pursuit, no data model or assumption about
independent components is made. If the ICA model holds, optimizing the ICA nongaussianity measures produce
independent components; if the model does not hold, then what we get are the projection pursuit directions.

5 Preprocessing for ICA
In the preceding section, we discussed the statistical principles underlying ICA methods. Practical algorithms
based on these principles will be discussed in the next section. However, before applying an ICA algorithm on thedata, it is usually very useful to do some preprocessing. In this section, we discuss some preprocessing techniquesthat make the problem of ICA estimation simpler and better conditioned.

5.1 Centering
The most basic and necessary preprocessing is to center x, i.e. subtract its mean vector m = E{x} so as to make x
a zero-mean variable. This implies that s is zero-mean as well, as can be seen by taking expectations on both sidesof Eq. (4).
This preprocessing is made solely to simplify the ICA algorithms: It does not mean that the mean could not be
estimated. After estimating the mixing matrix A with centered data, we can complete the estimation by adding themean vector of s back to the centered estimates of s. The mean vector of s is given by A−1 m, where m is the meanthat was subtracted in the preprocessing.

5.2 Whitening
Another useful preprocessing strategy in ICA is to first whiten the observed variables. This means that before
the application of the ICA algorithm (and after centering), we transform the observed vector x linearly so that we
obtain a new vector x̃ which is white, i.e. its components are uncorrelated and their variances equal unity. In otherwords, the covariance matrix of x̃ equals the identity matrix:E{x̃x̃T } = I.

(34)

The whitening transformation is always possible. One popular method for whitening is to use the eigen-value
decomposition (EVD) of the covariance matrix E{xxT } = EDET , where E is the orthogonal matrix of eigenvectorsof E{xxT } and D is the diagonal matrix of its eigenvalues, D = diag(d1 , ..., dn ). Note that E{xxT } can be estimatedin a standard way from the available sample x(1), ..., x(T ). Whitening can now be done byx̃ = ED−1/2 ET x

(35)
−1/2

where the matrix D−1/2 is computed by a simple component-wise operation as D−1/2 = diag(d1is easy to check that now E{x̃x̃T } = I.
12

−1/2

, ..., dn). ItWhitening transforms the mixing matrix into a new one, Ã. We have from (4) and (35):x̃ = ED−1/2 ET As = Ãs(36)

The utility of whitening resides in the fact that the new mixing matrix Ã is orthogonal. This can be seen fromE{x̃x̃T } = ÃE{ssT }ÃT = ÃÃT = I.

(37)

Here we see that whitening reduces the number of parameters to be estimated. Instead of having to estimate then2 parameters that are the elements of the original matrix A, we only need to estimate the new, orthogonal mixing
matrix Ã. An orthogonal matrix contains n(n − 1)/2 degrees of freedom. For example, in two dimensions, anorthogonal transformation is determined by a single angle parameter. In larger dimensions, an orthogonal matrix
contains only about half of the number of parameters of an arbitrary matrix. Thus one can say that whitening
solves half of the problem of ICA. Because whitening is a very simple and standard procedure, much simpler thanany ICA algorithms, it is a good idea to reduce the complexity of the problem this way.
It may also be quite useful to reduce the dimension of the data at the same time as we do the whitening. Then
we look at the eigenvalues d j of E{xxT } and discard those that are too small, as is often done in the statisticaltechnique of principal component analysis. This has often the effect of reducing noise. Moreover, dimension
reduction prevents overlearning, which can sometimes be observed in ICA (Hyvärinen et al., 1999).
A graphical illustration of the effect of whitening can be seen in Figure 10, in which the data in Figure 6
has been whitened. The square defining the distribution is now clearly a rotated version of the original square in
Figure 5. All that is left is the estimation of a single angle that gives the rotation.
In the rest of this paper, we assume that the data has been preprocessed by centering and whitening. For
simplicity of notation, we denote the preprocessed data just by x, and the transformed mixing matrix by A, omitting
the tildes.

5.3 Further preprocessing
The success of ICA for a given data set may depende crucially on performing some application-dependent preprocessing steps. For example, if the data consists of time-signals, some band-pass filtering may be very useful. Notethat if we filter linearly the observed signals xi (t) to obtain new signals, say x∗i (t), the ICA model still holds forx∗i (t), with the same mixing matrix.
This can be seen as follows. Denote by X the matrix that contains the observations x(1), ..., x(T ) as its columns,
and similarly for S. Then the ICA model can be expressed as:
X = AS(38)

Now, time filtering of X corresponds to multiplying X from the right by a matrix, let us call it M. This gives
X∗ = XM = ASM = AS∗ ,(39)which shows that the ICA model remains still valid.

6 The FastICA Algorithm
In the preceding sections, we introduced different measures of nongaussianity, i.e. objective functions for ICA
estimation. In practice, one also needs an algorithm for maximizing the contrast function, for example the one in(25). In this section, we introduce a very efficient method of maximization suited for this task. It is here assumed
that the data is preprocessed by centering and whitening as discussed in the preceding section.

13

6.1 FastICA for one unit
To begin with, we shall show the one-unit version of FastICA. By a "unit" we refer to a computational unit,
eventually an artificial neuron, having a weight vector w that the neuron is able to update by a learning rule. The
FastICA learning rule finds a direction, i.e. a unit vector w such that the projection wT x maximizes nongaussianity.
Nongaussianity is here measured by the approximation of negentropy J(wT x) given in (25). Recall that the varianceof wT x must here be constrained to unity; for whitened data this is equivalent to constraining the norm of w to beunity.
The FastICA is based on a fixed-point iteration scheme for finding a maximum of the nongaussianity of wT x,
as measured in (25), see (Hyvärinen and Oja, 1997; Hyvärinen, 1999a). It can be also derived as an approximative
Newton iteration (Hyvärinen, 1999a). Denote by g the derivative of the nonquadratic function G used in (25); forexample the derivatives of the functions in (26) are:g1 (u) = tanh(a1 u),
g2 (u) = u exp(−u2 /2)(40)

where 1 ≤ a1 ≤ 2 is some suitable constant, often taken as a1 = 1. The basic form of the FastICA algorithm is asfollows:
1. Choose an initial (e.g. random) weight vector w.
2. Let w+ = E{xg(wT x)} − E{g0(wT x)}w
3. Let w = w+ /kw+ k
4. If not converged, go back to 2.
Note that convergence means that the old and new values of w point in the same direction, i.e. their dot-product is(almost) equal to 1. It is not necessary that the vector converges to a single point, since w and −w define the same
direction. This is again because the independent components can be defined only up to a multiplicative sign. Notealso that it is here assumed that the data is prewhitened.
The derivation of FastICA is as follows. First note that the maxima of the approximation of the negentropy ofwT x are obtained at certain optima of E{G(wT x)}. According to the Kuhn-Tucker conditions (Luenberger, 1969),
the optima of E{G(wT x)} under the constraint E{(wT x)2 } = kwk2 = 1 are obtained at points whereE{xg(wT x)} − βw = 0(41)

Let us try to solve this equation by Newton’s method. Denoting the function on the left-hand side of (41) by F, we
obtain its Jacobian matrix JF(w) as
JF(w) = E{xxT g0 (wT x)} − βI
(42)
To simplify the inversion of this matrix, we decide to approximate the first term in (42). Since the data is sphered,a reasonable approximation seems to be E{xxT g0 (wT x)} ≈ E{xxT }E{g0(wT x)} = E{g0 (wT x)}I. Thus the Jacobian matrix becomes diagonal, and can easily be inverted. Thus we obtain the following approximative Newton
iteration:
w+ = w − [E{xg(wT x)} − βw]/[E{g0(wT x)} − β](43)This algorithm can be further simplified by multiplying both sides of (43) by β − E{g0(wT x)}. This gives, after
algebraic simplication, the FastICA iteration.
In practice, the expectations in FastICA must be replaced by their estimates. The natural estimates are ofcourse the corresponding sample means. Ideally, all the data available should be used, but this is often not a goodidea because the computations may become too demanding. Then the averages can be estimated using a smaller
sample, whose size may have a considerable effect on the accuracy of the final estimates. The sample points should
be chosen separately at every iteration. If the convergence is not satisfactory, one may then increase the sample
size.
14

6.2 FastICA for several units
The one-unit algorithm of the preceding subsection estimates just one of the independent components, or one
projection pursuit direction. To estimate several independent components, we need to run the one-unit FastICAalgorithm using several units (e.g. neurons) with weight vectors w1 , ..., wn .
To prevent different vectors from converging to the same maxima we must decorrelate the outputs wT1 x, ..., wTn x
after every iteration. We present here three methods for achieving this.
A simple way of achieving decorrelation is a deflation scheme based on a Gram-Schmidt-like decorrelation.
This means that we estimate the independent components one by one. When we have estimated p independent
components, or p vectors w1 , ..., w p , we run the one-unit fixed-point algorithm for w p+1 , and after every iteration
step subtract from w p+1 the “projections” wTp+1w j w j , j = 1, ..., p of the previously estimated p vectors, and thenrenormalize w p+1 :p
1. Let w p+1 = w p+1 − ∑ j=1 wTp+1 w j w jq(44)
2. Let w p+1 = w p+1 / wTp+1 w p+1
In certain applications, however, it may be desired to use a symmetric decorrelation, in which no vectors are“privileged” over others (Karhunen et al., 1997). This can be accomplished, e.g., by the classical method involving
matrix square roots,
Let W = (WWT )−1/2 W(45)

where W is the matrix (w1 , ..., wn )T of the vectors, and the inverse square root (WWT )−1/2 is obtained from theeigenvalue decomposition of WWT = FDFT as (WWT )−1/2 = FD−1/2FT . A simpler alternative is the following
iterative algorithm (Hyvärinen, 1999a),p
1. Let W = W/ kWWT k(46)
Repeat 2. until convergence:
1
3
T
2. Let W = 2 W − 2 WW W
The norm in step 1 can be almost any ordinary matrix norm, e.g., the 2-norm or the largest absolute row (or column)sum (but not the Frobenius norm).

6.3 FastICA and maximum likelihood
Finally, we give a version of FastICA that shows explicitly the connection to the well-known infomax or maximum
likelihood algorithm introduced in (Amari et al., 1996; Bell and Sejnowski, 1995; Cardoso and Laheld, 1996;
Cichocki and Unbehauen, 1996). If we express FastICA using the intermediate formula in (43), and write it inmatrix form (see (Hyvärinen, 1999b) for details), we see that FastICA takes the following form:
W+ = W + diag(αi )[diag(βi ) + E{g(y)yT }]W.

(47)

where y = Wx, βi = −E{yi g(yi )}, and αi = −1/(βi − E{g0(yi )}). The matrix W needs to be orthogonalized after
every step. In this matrix version, it is natural to orthogonalize W symmetrically.
The above version of FastICA could be compared with the stochastic gradient method for maximizing likelihood (Amari et al., 1996; Bell and Sejnowski, 1995; Cardoso and Laheld, 1996; Cichocki and Unbehauen, 1996):
W+ = W + µ[I + g(y)yT ]W.

(48)

where µ is the learning rate, not necessarily constant in time. Here, g is a function of the pdf of the independent
components: g = fi0 / fi where the fi is the pdf of an independent component.
Comparing (47) and (48), we see that FastICA can be considered as a fixed-point algorithm for maximum
likelihood estimation of the ICA data model. For details, see (Hyvärinen, 1999b). In FastICA, convergence speed
is optimized by the choice of the matrices diag(αi ) and diag(βi ). Another advantage of FastICA is that it canestimate both sub- and super-gaussian independent components, which is in contrast to ordinary ML algorithms,
which only work for a given class of distributions (see Sec. 4.4).
15

6.4 Properties of the FastICA Algorithm
The FastICA algorithm and the underlying contrast functions have a number of desirable properties when comparedwith existing methods for ICA.
1. The convergence is cubic (or at least quadratic), under the assumption of the ICA data model (for a proof, see(Hyvärinen, 1999a)). This is in contrast to ordinary ICA algorithms based on (stochastic) gradient descentmethods, where the convergence is only linear. This means a very fast convergence, as has been confirmedby simulations and experiments on real data (see (Giannakopoulos et al., 1998)).
2. Contrary to gradient-based algorithms, there are no step size parameters to choose. This means that thealgorithm is easy to use.
3. The algorithm finds directly independent components of (practically) any non-Gaussian distribution using
any nonlinearity g. This is in contrast to many algorithms, where some estimate of the probability distributionfunction has to be first available, and the nonlinearity must be chosen accordingly.
4. The performance of the method can be optimized by choosing a suitable nonlinearity g. In particular, one
can obtain algorithms that are robust and/or of minimum variance. In fact, the two nonlinearities in (40)have some optimal properties; for details see (Hyvärinen, 1999a).
5. The independent components can be estimated one by one, which is roughly equivalent to doing projection
pursuit. This es useful in exploratory data analysis, and decreases the computational load of the method incases where only some of the independent components need to be estimated.
6. The FastICA has most of the advantages of neural algorithms: It is parallel, distributed, computationallysimple, and requires little memory space. Stochastic gradient methods seem to be preferable only if fast
adaptivity in a changing environment is required.
A MatlabT M implementation of the FastICA algorithm is available on the World Wide Web free of charge.1

7 Applications of ICA
In this section we review some applications of ICA. The most classical application of ICA, the cocktail-party
problem, was already explained in the opening section of this paper.

7.1 Separation of Artifacts in MEG Data
Magnetoencephalography (MEG) is a noninvasive technique by which the activity or the cortical neurons can bemeasured with very good temporal resolution and moderate spatial resolution. When using a MEG record, as a
research or clinical tool, the investigator may face a problem of extracting the essential features of the neuromagnetic signals in the presence of artifacts. The amplitude of the disturbances may be higher than that of the brain
signals, and the artifacts may resemble pathological signals in shape.
In (Vigário et al., 1998), the authors introduced a new method to separate brain activity from artifacts using
ICA. The approach is based on the assumption that the brain activity and the artifacts, e.g. eye movements or blinks,
or sensor malfunctions, are anatomically and physiologically separate processes, and this separation is reflected in
the statistical independence between the magnetic signals generated by those processes. The approach follows theearlier experiments with EEG signals, reported in (Vigário, 1997). A related approach is that of (Makeig et al.,1996).
The MEG signals were recorded in a magnetically shielded room with a 122-channel whole-scalp Neuromag122 neuromagnetometer. This device collects data at 61 locations over the scalp, using orthogonal double-loop
pick-up coils that couple strongly to a local source just underneath. The test person was asked to blink and makehorizontal saccades, in order to produce typical ocular (eye) artifacts. Moreover, to produce myographic (muscle)
1 WWW address: http://www.cis.hut.fi/projects/ica/fastica/

16

artifacts, the subject was asked to bite his teeth for as long as 20 seconds. Yet another artifact was created byplacing a digital watch one meter away from the helmet into the shielded room.
Figure 11 presents a subset of 12 spontaneous MEG signals xi (t) from the frontal, temporal, and occipitalareas (Vigário et al., 1998). The figure also shows the positions of the corresponding sensors on the helmet. Dueto the dimension of the data (122 magnetic signals were recorded), it is impractical to plot all the MEG signals
xi (t), i = 1, ..., 122. Also two electro-oculogram channels and the electrocardiogram are presented, but they were
not used in computing the ICA.
The signal vector x in the ICA model (4) consists now of the amplitudes xi (t) of the 122 signals at a certain
time point, so the dimensionality is n = 122. In the theoretical model, x is regarded as a random vector, and themeasurements x(t) give a set of realizations of x as time proceeds. Note that in the basic ICA model that we areusing, the temporal correlations in the signals are not utilized at all.
The x(t) vectors were whitened using PCA and the dimensionality was decreased at the same time. Then, using
the FastICA algorithm, a subset of the rows of the separating matrix W of eq. (6) were computed. Once a vector
wi has become available, an ICA signal si (t) can be computed from si (t) = wTi x0 (t) with x0 (t) now denoting thewhitened and lower dimensional signal vector.
Figure 12 shows sections of 9 independent components (IC’s) si (t), i = 1, ..., 9 found from the recorded data
together with the corresponding field patterns (Vigário et al., 1998). The first two IC’s are clearly due to themusclular activity originated from the biting. Their separation into two components seems to correspond, on thebasis of the field patterns, to two different sets of muscles that were activated during the process. IC3 and IC5 areshowing the horizontal eye movements and the eye blinks, respectively. IC4 represents the cardiac artifact that isvery clearly extracted.
To find the remaining artifacts, the data were high-pass filtered, with cutoff frequency at 1 Hz. Next, theindependent component IC8 was found. It shows clearly the artifact originated at the digital watch, located to theright side of the magnetometer. The last independent component IC9 is related to a sensor presenting higher RMS(root mean squared) noise than the others.
The results of Fig. 12 clearly show that using the ICA technique and the FastICA algorithm, it is possible toisolate both eye movement and eye blinking artifacts, as well as cardiac, myographic, and other artifacts from MEG
signals. The FastICA algorithm is an especially suitable tool, because artifact removal is an interactive techniqueand the investigator may freely choose how many of the IC’s he or she wants.
In addition to reducing artifacts, ICA can be used to decompose evoked fields (Vigário et al., 1998), which enables direct access to the underlying brain functioning, which is likely to be of great significance in neuroscientificresearch.

7.2 Finding Hidden Factors in Financial Data
It is a tempting alternative to try ICA on financial data. There are many situations in that application domain inwhich parallel time series are available, such as currency exchange rates or daily returns of stocks, that may havesome common underlying factors. ICA might reveal some driving mechanisms that otherwise remain hidden. In
a recent study of a stock portfolio (Back and Weigend, 1997), it was found that ICA is a complementary tool to
PCA, allowing the underlying structure of the data to be more readily observed.
In (Kiviluoto and Oja, 1998), we applied ICA on a different problem: the cashflow of several stores belonging
to the same retail chain, trying to find the fundamental factors common to all stores that affect the cashflow data.
Thus, the cashflow effect of the factors specific to any particular store, i.e., the effect of the actions taken at theindividual stores and in its local environment could be analyzed.
The assumption of having some underlying independent components in this specific application may not be
unrealistic. For example, factors like seasonal variations due to holidays and annual variations, and factors having
a sudden effect on the purchasing power of the customers like prize changes of various commodities, can be
expected to have an effect on all the retail stores, and such factors can be assumed to be roughly independent ofeach other. Yet, depending on the policy and skills of the individual manager like e.g. advertising efforts, the effectof the factors on the cash flow of specific retail outlets are slightly different. By ICA, it is possible to isolate boththe underlying factors and the effect weights, thus also making it possible to group the stores on the basis of their
managerial policies using only the cash flow time series data.

17

The data consisted of the weekly cash flow in 40 stores that belong to the same retail chain; the cash flow
measurements cover 140 weeks. Some examples of the original data xi (t) are shown in Fig. 13.
The prewhitening was performed so that the original signal vectors were projected to the subspace spannedby their first five principal components and the variances were normalized to 1. Thus the dimension of the signal
space was decreased from 40 to 5. Using the FastICA algorithm, four IC’s si (t), i = 1, ..., 5 were estimated. As
depicted in Fig. 14, the FastICA algorithm has found several clearly different fundamental factors hidden in theoriginal data.
The factors have clearly different interpretations. The upmost two factors follow the sudden changes that arecaused by holidays etc.; the most prominent example is the Christmas time. The factor on the bottom row, on theother hand, reflects the slower seasonal variation, with the effect of the summer holidays clearly visible. The factor
on the third row could represent a still slower variation, something resembling a trend. The last factor, on the fourthrow, is different from the others; it might be that this factor follows mostly the relative competitive position of theretail chain with respect to its competitors, but other interpretations are also possible.
More details on the experiments and their interpretation can be found in (Kiviluoto and Oja, 1998).

7.3 Reducing Noise in Natural Images
The third example deals with finding ICA filters for natural images and, based on the ICA decomposition, removing
noise from images corrupted with additive Gaussian noise.
A set of digitized natural images were used. Denote the vector of pixel gray levels in an image window by
x. Note that, contrary to the other two applications in the previous sections, we are not this time considering
multivalued time series or images changing with time; instead the elements of x are indexed by the location in theimage window or patch. The sample windows were taken at random locations. The 2-D structure of the windows
is of no significance here: row by row scanning was used to turn a square image window into a vector of pixel
values. The independent components of such image windows are represented in Fig. 4. Each window in this Figure
corresponds to one of the columns ai of the mixing matrix A. Thus an observed image window is a superpositionof these windows as in (5), with independent coefficients (Bell and Sejnowski, 1997; Olshausen and Field, 1996).
Now, suppose a noisy image model holds:z = x+n(49)
where n is uncorrelated noise, with elements indexed in the image window in the same way as x, and z is themeasured image window corrupted with noise. Let us further assume that n is Gaussian and x is non-Gaussian.
There are many ways to clean the noise; one example is to make a transformation to spatial frequency space by
DFT, do low-pass filtering, and return to the image space by IDFT (Gonzales and Wintz, 1987). This is not very
efficient, however. A better method is the recently introduced Wavelet Shrinkage method (Donoho et al., 1995) inwhich a transform based on wavelets is used, or methods based on median filtering (Gonzales and Wintz, 1987).
None of these methods is explicitly taking advantage of the image statistics, however.
We have recently introduced another, statistically principled method called Sparse Code Shrinkage (Hyvärinen,
1999d). It is very closely related to independent component analysis. Briefly, if we model the density of x by ICA,
and assume n Gaussian, then the Maximum Likelihood (ML) solution for x given the measurement z can be
developed in the signal model (49).
The ML solution can be simply computed, albeit approximately, by using a decomposition that is an orthogonalized version of ICA. The transform is given by
Wz = Wx + Wn = s + Wn,(50)

where W is here an orthogonal matrix that is the best orthognal approximation of the inverse of the ICA mixing
matrix. The noise term Wn is still Gaussian and white. With a suitably chosen orthogonal transform W, however,the density of Wx = s becomes highly non-Gaussian, e.g., super-Gaussian with a high positive kurtosis. Thisdepends of course on the original x signals, as we are assuming in fact that there exists a model x = WT s forthe signal, such that the “source signals” or elements of s have a positive kurtotic density, in which case the ICA
transform gives highly supergaussian components. This seems to hold at least for image windows of natural scenes(Mallat, 1989).
18

It was shown in (Hyvärinen, 1999d) that, assuming a Laplacian density for si , the ML solution for si is given
by a “shrinkage function” sˆi = g([Wz]i ), or in vector form, ŝ = g(Wz). Function g(.) has a characteristic shape: itis zero close to the origin and then linear after a cutting value depending on the parameters of the Laplacian density
and the Gaussian noise density. Assuming other forms for the densities, other optimal shrinkage functions can bederived (Hyvärinen, 1999d).
In the Sparse Code Shrinkage method, the shrinkage operation is performed in the rotated space, after which
the estimate for the signal in the original space is given by rotating back:x̂ = WT ŝ = WT g(Wz).

(51)

Thus we get the Maximum Likelihood estimate for the image window x in which much of the noise has been
removed.
The rotation operator W is such that the sparsity of the components s = Wx is maximized. This operator canbe learned with a modification of the FastICA algorithm; see (Hyvärinen, 1999d) for details.
A noise cleaning result is shown in Fig. 15. A noiseless image and a noisy version, in which the noise level is
50 % of the signal level, are shown. The results of the Sparse Code Shrinkage method and classic wiener filtering
are given, indicating that Sparse Code Shrinkage may be a promising approach. The noise is reduced without
blurring edges or other sharp features as much as in wiener filtering. This is largely due to the strongly nonlinear
nature of the shrinkage operator, that is optimally adapted to the inherent statistics of natural images.

7.4 Telecommunications
Finally, we mention another emerging application area of great potential: telecommunications. An example ofa real-world communications application where blind separation techniques are useful is the separation of theuser’s own signal from the interfering other users’ signals in CDMA (Code-Division Multiple Access) mobile
communications (Ristaniemi and Joutsensalo, 1999). This problem is semi-blind in the sense that certain additional
prior information is available on the CDMA data model. But the number of parameters to be estimated is often sohigh that suitable blind source separation techniques taking into account the available prior knowledge provide aclear performance improvement over more traditional estimation techniques (Ristaniemi and Joutsensalo, 1999).

8 Conclusion
ICA is a very general-purpose statistical technique in which observed random data are linearly transformed into
components that are maximally independent from each other, and simultaneously have “interesting” distributions.
ICA can be formulated as the estimation of a latent variable model. The intuitive notion of maximum nongaussianity can be used to derive different objective functions whose optimization enables the estimation of the ICA
model. Alternatively, one may use more classical notions like maximum likelihood estimation or minimization ofmutual information to estimate ICA; somewhat surprisingly, these approaches are (approximatively) equivalent. A
computationally very efficient method performing the actual estimation is given by the FastICA algorithm. Applications of ICA can be found in many different areas such as audio processing, biomedical signal processing, image
processing, telecommunications, and econometrics.

References
Amari, S.-I., Cichocki, A., and Yang, H. (1996). A new learning algorithm for blind source separation. In
Advances in Neural Information Processing Systems 8, pages 757–763. MIT Press.
Back, A. D. and Weigend, A. S. (1997). A first application of independent component analysis to extracting
structure from stock returns. Int. J. on Neural Systems, 8(4):473–484.
Bell, A. and Sejnowski, T. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7:1129–1159.

19

Bell, A. and Sejnowski, T. (1997). The ’independent components’ of natural scenes are edge filters. Vision
Research, 37:3327–3338.
Cardoso, J.-F. (1997). Infomax and maximum likelihood for source separation. IEEE Letters on Signal Processing, 4:112–114.
Cardoso, J.-F. and Laheld, B. H. (1996). Equivariant adaptive source separation. IEEE Trans. on Signal Processing, 44(12):3017–3030.
Cichocki, A. and Unbehauen, R. (1996). Robust neural networks with on-line learning for blind identification
and blind separation of sources. IEEE Trans. on Circuits and Systems, 43(11):894–906.
Comon, P. (1994). Independent component analysis—a new concept? Signal Processing, 36:287–314.
Cover, T. M. and Thomas, J. A. (1991). Elements of Information Theory. Wiley.
Delfosse, N. and Loubaton, P. (1995). Adaptive blind separation of independent sources: a deflation approach.
Signal Processing, 45:59–83.
Donoho, D. L., Johnstone, I. M., Kerkyacharian, G., and Picard, D. (1995). Wavelet shrinkage: asymptopia?
Journal of the Royal Statistical Society, Ser. B, 57:301–337.
Friedman, J. (1987). Exploratory projection pursuit. J. of the American Statistical Association, 82(397):249–
266.
Friedman, J. H. and Tukey, J. W. (1974). A projection pursuit algorithm for exploratory data analysis. IEEE
Trans. of Computers, c-23(9):881–890.
Giannakopoulos, X., Karhunen, J., and Oja, E. (1998). Experimental comparison of neural ICA algorithms. In
Proc. Int. Conf. on Artificial Neural Networks (ICANN’98), pages 651–656, Skövde, Sweden.
Gonzales, R. and Wintz, P. (1987). Digital Image Processing. Addison-Wesley.
Huber, P. (1985). Projection pursuit. The Annals of Statistics, 13(2):435–475.
Hyvärinen, A. (1998a). Independent component analysis in the presence of gaussian noise by maximizing joint
likelihood. Neurocomputing, 22:49–67.
Hyvärinen, A. (1998b). New approximations of differential entropy for independent component analysis andprojection pursuit. In Advances in Neural Information Processing Systems, volume 10, pages 273–279. MIT
Press.
Hyvärinen, A. (1999a). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3):626–634.
Hyvärinen, A. (1999b). The fixed-point algorithm and maximum likelihood estimation for independent component analysis. Neural Processing Letters, 10(1):1–5.
Hyvärinen, A. (1999c). Gaussian moments for noisy independent component analysis. IEEE Signal Processing
Letters, 6(6):145–147.
Hyvärinen, A. (1999d). Sparse code shrinkage: Denoising of nongaussian data by maximum likelihood estimation. Neural Computation, 11(7):1739–1768.
Hyvärinen, A. (1999e). Survey on independent component analysis. Neural Computing Surveys, 2:94–128.
Hyvärinen, A. and Oja, E. (1997). A fast fixed-point algorithm for independent component analysis. NeuralComputation, 9(7):1483–1492.
Hyvärinen, A. and Oja, E. (1998). Independent component analysis by general nonlinear Hebbian-like learning
rules. Signal Processing, 64(3):301–313.
Hyvärinen, A., Särelä, J., and Vigário, R. (1999). Spikes and bumps: Artefacts generated by independent component analysis with insufficient sample size. In Proc. Int. Workshop on Independent Component Analysis
and Signal Separation (ICA’99), pages 425–429, Aussois, France.
Jones, M. and Sibson, R. (1987). What is projection pursuit ? J. of the Royal Statistical Society, Ser. A, 150:1–36.

20

Jutten, C. and Hérault, J. (1991). Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture. Signal Processing, 24:1–10.
Karhunen, J., Oja, E., Wang, L., Vigário, R., and Joutsensalo, J. (1997). A class of neural networks for independent component analysis. IEEE Trans. on Neural Networks, 8(3):486–504.
Kiviluoto, K. and Oja, E. (1998). Independent component analysis for parallel financial time series. In Proc. Int.
Conf. on Neural Information Processing (ICONIP’98), volume 2, pages 895–898, Tokyo, Japan.
Lee, T.-W., Girolami, M., and Sejnowski, T. J. (1999). Independent component analysis using an extended
infomax algorithm for mixed sub-gaussian and super-gaussian sources. Neural Computation, 11(2):417–
441.
Luenberger, D. (1969). Optimization by Vector Space Methods. Wiley.
Makeig, S., Bell, A. J., Jung, T.-P., and Sejnowski, T. (1996). Independent component analysis of electroencephalographic data. In Touretzky, D. S., Mozer, M. C., and Hasselmo, M. E., editors, Advances in Neural
Information Processing Systems 8, Cambridge MA. MIT Press.
Mallat, S. G. (1989). A theory for multiresolution signal decomposition: The wavelet representation. IEEE
Trans. on Pattern Analysis and Machine Intelligence, 11:674–693.
Nadal, J.-P. and Parga, N. (1994). Non-linear neurons in the low noise limit: a factorial code maximizes information transfer. Network, 5:565–581.
Olshausen, B. A. and Field, D. J. (1996). Emergence of simple-cell receptive field properties by learning asparse code for natural images. Nature, 381:607–609.
Papoulis, A. (1991). Probability, Random Variables, and Stochastic Processes. McGraw-Hill, 3rd edition.
Pearlmutter, B. A. and Parra, L. C. (1997). Maximum likelihood blind source separation: A context-sensitive
generalization of ICA. In Advances in Neural Information Processing Systems, volume 9, pages 613–619.
Pham, D.-T., Garrat, P., and Jutten, C. (1992). Separation of a mixture of independent sources through a maximum likelihood approach. In Proc. EUSIPCO, pages 771–774.
Ristaniemi, T. and Joutsensalo, J. (1999). On the performance of blind source separation in CDMA downlink. In
Proc. Int. Workshop on Independent Component Analysis and Signal Separation (ICA’99), pages 437–441,
Aussois, France.
Vigário, R. (1997). Extraction of ocular artifacts from EEG using independent component analysis. Electroenceph. Clin. Neurophysiol., 103(3):395–404.
Vigário, R., Jousmäki, V., Hämäläinen, M., Hari, R., and Oja, E. (1998). Independent component analysis
for identification of artifacts in magnetoencephalographic recordings. In Advances in Neural Information
Processing Systems, volume 10, pages 229–235. MIT Press.
Vigário, R., Särelä, J., and Oja, E. (1998). Independent component analysis in wave decomposition of auditory evoked fields. In Proc. Int. Conf. on Artificial Neural Networks (ICANN’98), pages 287–292, Skövde,Sweden.

21

2

0

−2

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

2

0

−2

Figure 1: The original signals.

5

0

−5

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

2

0

−2

Figure 2: The observed mixtures of the source signals in Fig. 1.

22

2

0

−2

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

2

0

−2

Figure 3: The estimates of the original source signals, estimated using only the observed signals in Fig. 2. The
original signals were very accurately estimated, up to multiplicative signs.

23

Figure 4: Basis functions in ICA of natural images. The input window size was 16 × 16 pixels. These basis
functions can be considered as the independent features of images.

24

Figure 5: The joint distribution of the independent components s1 and s2 with uniform distributions. Horizontal
axis: s1 , vertical axis: s2 .

Figure 6: The joint distribution of the observed mixtures x1 and x2 . Horizontal axis: x1 , vertical axis: x2 .

25

Figure 7: The multivariate distribution of two independent gaussian variables.

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
−4

−3

−2

−1

0

1

2

3

4

Figure 8: The density function of the Laplace distribution, which is a typical supergaussian distribution. For
comparison, the gaussian density is given by a dashed line. Both densities are normalized to unit variance.

26

Figure 9: An illustration of projection pursuit and the “interestingness” of nongaussian projections. The data inthis figure is clearly divided into two clusters. However, the principal component, i.e. the direction of maximum
variance, would be vertical, providing no separation between the clusters. In contrast, the strongly nongaussian
projection pursuit direction is horizontal, providing optimal separation of the clusters.

Figure 10: The joint distribution of the whitened mixtures.

27

2

3
4

saccades

5

5

6

6

blinking

1

MEG

1000 fT/cm

EOG

500 µVECG

500 µV
biting

MEG
1
1
2
2
3
3
4
4
5
5
6
6
VEOG
HEOG
ECG

10 s

Figure 11: (From Vigário et al, 1998). Samples of MEG signals, showing artifacts produced by blinking, saccades,
biting and cardiac cycle. For each of the 6 positions shown, the two orthogonal directions of the sensors areplotted.

28

IC1

IC2

IC3

IC4

IC5

IC6

IC7

IC8

IC9
10 s

Figure 12: (From Vigário et al, 1998). Nine independent components found from the MEG data. For each component the left, back and right views of the field patterns generated by these components are shown — full line stands
for magnetic flux coming out from the head, and dotted line the flux inwards.

29

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

Figure 13: (from Kiviluoto and Oja, 1998). Five samples of the original cashflow time series (mean removed,
normalized to unit standard deviation). Horizontal axis: time in weeks.

10

1

0
−10

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

10

2

0
−10

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

3

2
0
−2

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

4

5
0
−5

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

5

5
0
−5

1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49 1 5 9 13 17 21 25 29 33 37 41 45 49

Figure 14: (from Kiviluoto and Oja, 1998). Four independent components or fundamental factors found from the
cashflow data.

30

Figure 15: (from Hyvärinen, 1999d). An experiment in denoising. Upper left: original image. Upper right:
original image corrupted with noise; the noise level is 50 %. Lower left: the recovered image after applying
sparse code shrinkage. Lower right: for comparison, a wiener filtered image.

31