Journal of Machine Learning Research 9 (2008) 2455-2487

Submitted 9/07; Revised 5/08; Published 11/08

Minimal Nonlinear Distortion Principle for Nonlinear Independent
Component Analysis
Kun Zhang
Laiwan ChanKZHANG @ CSE . CUHK . EDU . HK
LWCHAN @ CSE . CUHK . EDU . HK

Department of Computer Science and Engineering
The Chinese University of Hongkong
Hong Kong

Editor: Aapo Hyvärinen

Abstract
It is well known that solutions to the nonlinear independent component analysis (ICA) problem arehighly non-unique. In this paper we propose the “minimal nonlinear distortion” (MND) principle
for tackling the ill-posedness of nonlinear ICA problems. MND prefers the nonlinear ICA solution
with the estimated mixing procedure as close as possible to linear, among all possible solutions.
It also helps to avoid local optima in the solutions. To achieve MND, we exploit a regularization
term to minimize the mean square error between the nonlinear mixing mapping and the best-fitting
linear one. The effect of MND on the inherent trivial and non-trivial indeterminacies in nonlinear
ICA solutions is investigated. Moreover, we show that local MND is closely related to the smoothness regularizer penalizing large curvature, which provides another useful regularization condition
for nonlinear ICA. Experiments on synthetic data show the usefulness of the MND principle forseparating various nonlinear mixtures. Finally, as an application, we use nonlinear ICA with MND
to separate daily returns of a set of stocks in Hong Kong, and the linear causal relations among
them are successfully discovered. The resulting causal relations give some interesting insights into
the stock market. Such a result can not be achieved by linear ICA. Simulation studies also verify
that when doing causality discovery, sometimes one should not ignore the nonlinear distortion in
the data generation procedure, even if it is weak.
Keywords: nonlinear ICA, regularization, minimal nonlinear distortion, mean square error, best
linear reconstruction

1. Introduction
Independent component analysis (ICA) is a popular statistical technique aiming to recover independent sources from their observed mixtures, without knowing the mixing procedure or any specificknowledge of the sources (Hyvärinen et al., 2001; Cardoso, 1998; Cichocki and Amari, 2003). In
the case that the observed mixtures are a linear transformation of the sources, under weak assumptions, ICA can recover the original sources with the trivial permutation and scaling indeterminacies.
Linear ICA is currently a popular method for blind source separation (BSS) of linear mixtures.
However, nonlinear ICA does not necessarily lead to nonlinear BSS. In Hyvärinen and Pajunen (1999), it was shown that solutions to nonlinear ICA always exist and that they are highly
non-unique. Actually, one can easily construct a nonlinear transformation of some non-Gaussian
independent variables to produce independent outputs. Below are a few examples. Let y 1 , ..., yn besome independent variables. Their component-wise nonlinear functions are still mutually indepenc 2008 Kun Zhang and Laiwan Chan.

Z HANG AND C HAN

dent. If we use Gaussianization (Chen and Gopinath, 2001) to transform y i into Gaussian variablesui , any component-wise nonlinear function of U · u, where u = (u 1 , ..., un )T and U is an orthogonal
matrix, still has mutually independent components. Taleb and Jutten (1999) also gave an example inwhich nonlinear mixtures of independent variables are still independent. One can see that nonlinear
BSS is impossible without additional prior knowledge on the mixing model, since the independence
assumption is not strong enough in the general nonlinear mixing case (Jutten and Taleb, 2000; Taleb,
2002).
If we constrain the nonlinear mixing mapping to have some specific forms, the indeterminacies
in the results of nonlinear ICA can be reduced dramatically, and as a consequence, nonlinear ICA
may lead to nonlinear BSS. For example, in Burel (1992), a parametric form of the mixing transformation is assumed known and one just needs to adjust the unknown parameters. The learningalgorithms were improved in Yang et al. (1998). By exploiting the extensions of the DarmoisSkitovich theorem (Kagan et al., 1973) to nonlinear functions, a particular class of nonlinear mixing
mappings, which satisfy an addition theorem in the sense of the theory of functional equations,were considered in Eriksson and Koivunen (2002). In particular, the post-nonlinear (PNL) mixing
model (Taleb and Jutten, 1999), which assumes that the mixing mapping is a linear transformation
followed by a component-wise nonlinear one, has drawn much attention.
In practice, the exact form of the nonlinear mixing procedure is probably unknown. Consequently, in order to model arbitrary nonlinear mappings, one may need to resort to a flexible nonlinear function approximator, such as the multi-layer perceptron (MLP) or the radial basis function(RBF) network, to represent the nonlinear separation system. Almeida (2003) uses the MLP tomodel the separation system and trains the MLP by information-maximization (Infomax). Moreover, the smoothness constraint,1 which is implicitly provided by MLP’s with small initial weights
and with a relatively small number of hidden units, was believed to be a suitable regularization
condition to achieve nonlinear BSS. In Tan et al. (2001), a RBF network is adopted to represent
the separation system, and partial moments of the outputs of the separation system are used forregularization. The matching between the relevant moments of the outputs and those of the original
sources was expected to guarantee a unique solution. But the moments of the original sources may
be unknown. In addition, if the transformation from the original sources to the recovered sources is
non-trivial,2 this regularization could not help to recover the original sources. Variational Bayesiannonlinear ICA (Lappalainen and Honkela, 2000; Valpola, 2000) uses the MLP to model the nonlinear mixing transformation. By resorting to the variational Bayesian inference technique, thismethod can do model selection and avoid overfitting. If we can have some additional knowledge
about the nonlinear mixing transformation and incorporate it efficiently, the results of nonlinear ICA
will be much more meaningful and reliable.
Although we may not know the form of the nonlinearity in the data generation procedure, fortunately, in many cases the nonlinearity for generating natural signals we deal with is not strong.
Hence, provided that the nonlinear ICA outputs are mutually independent, we would prefer the solution with the estimated data generation procedure of minimal nonlinear distortion (MND). This
1. Following Tikhonov and Arsenin (1977), here we use the term “smoothness” in a very general sense. Often it meansthat that the function does not change abruptly and/or that it does not oscillate too much.
2. For the definition of a trivial transformation, one may see Jutten and Taleb (2000). A one-to-one mapping H istrivial if and only if it satisfies Hi (y1 , y2 , ..., yn ) = hi (yσ(i) ), i = 1, 2, ..., n, where hi are arbitrary functions and σ is any
permutation over {1, .2, ..., n}. That is, a trivial mapping of y is a permutation of y i followed by a component-wisetransformation.

2456

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

information can help to reduce the indeterminacies in nonlinear ICA greatly, and moreover, to avoid
local optima in the solutions to nonlinear ICA. The minimal nonlinear distortion of the mixing system is achieved by the technique of regularization. The objective function of nonlinear ICA with
MND is the mutual information between outputs penalized by some terms measuring the level of“closeness to linear” of the mixing system. The mean square error (MSE) between the nonlinear
mixing system and its best-fitting linear one provides such a regularization term. To ensure that
nonlinear ICA results in nonlinear BSS, one may also need to enforce the local MND of the nonlinear mapping averaged at every point, which turns out to be the smoothness regularizer exploiting
second-order partial derivatives.
MND, as well as the smoothness regularizer, can be incorporated in various nonlinear ICA
methods to improve the results. Here we consider two nonlinear ICA methods. The first one is the
MISEP method (Almeida, 2003), where the MLP is used to represent the separation system. As
regularization is powerful for complexity control in neural networks (Bishop, 1995), the structureof the MLP is not optimized during the learning process, that is, it is fixed. The second one is nonlinear ICA based on kernels (Zhang and Chan, 2007a), in which the nonlinear separation system ismodeled using some kernel methods. We then explain why MND helps to alleviate the ill-posedness
in nonlinear ICA solutions, by investigating the effect of MND on trivial and non-trivial indeterminacies in nonlinear ICA solutions systematically. Next, we conduct experiments using synthetic
data to compare the performance of several nonlinear ICA methods. The results confirm the effectiveness of the proposed MND principle to avoid unwanted solutions and to improve the separation
performance. Finally, nonlinear ICA with MND is used to discover linear causal relations in the
Hong Kong stock market and give encouraging results. We also give experimental results on synthetic data, which illustrate that when performing ICA-based causality discovery on the data whose
generation procedure involves nonlinear distortion, one should take into account the nonlinear effectin the ICA separation system, even if it is mild.3

2. Nonlinear ICA with Minimal Nonlinear Distortion
In this section we first briefly review the general nonlinear ICA problem, and then propose theminimal nonlinear distortion (MND) principle for regularization of nonlinear ICA.
2.1 Nonlinear ICA
In the nonlinear ICA model, the observed data x = (x1 , ..., xn )T are assumed to be generated from avector of independent variables s = (s1 , ..., sn )T by a nonlinear transformation:
x = F (s),(1)

where F is an unknown real-valued n-component mixing function. Here for simplicity, we haveassumed that the number of observed variables equals that of the original independent variables.
The general nonlinear ICA problem is to find a mapping G : R n → Rn such thaty = G (x)has statistically independent components. As mentioned in Section 1, the results of nonlinear ICA
are highly non-unique. In order to achieve nonlinear BSS, which aims at recovering the original
sources si , we should resort to additional prior information or suitable regularization constraints.
3. Some preliminary results of this paper were presented at ICML2007 (Zhang and Chan, 2007b).

2457

Z HANG AND C HAN

2.2 With Minimum Nonlinear Distortion
We now propose the MND principle to restrict the space of nonlinear ICA solutions. As a consequence, the ill-posedness of the nonlinear ICA problem is alleviated. Under the condition that the
separation outputs yi are mutually independent, this principle prefers the solution with the estimatedmixing transformation F̂ as close as possible to linear.
Now we need a measure of “closeness to linear” of a mapping. Given a nonlinear mapping F̂ ,its deviation from the affine mapping A∗ , which fits F̂ best among all affine mappings A, is anindicator of its “closeness to linear”, or the level of its nonlinear distortion. The deviation can bemeasured in various ways. The MSE is adopted here, as it greatly facilitates subsequent analysis.
Consequently, the “closeness to linear” of F̂ = G −1 can be measured by the MSE between G −1
and A∗ . We denote this measure by RMSE (θ), where θ denotes the set of unknown parameters in the
nonlinear ICA system. Let x∗ = (x1∗ , · · · , xn∗ )T be the output of the affine transformation from y byA∗ . Let ỹ = [y; 1]. RMSE (θ) can then be written as the MSE between xi and xi∗ :RMSE (θ) = E{(x − x∗ )T (x − x∗ )} , where(2)x∗ = A∗ ỹ, and A∗ = argA min E{(x − Ay)T (x − Ay)}.

Here A∗ is a n × (n + 1) matrix.4 Figure 1 shows the separation system G together with the generation process of RMSE .

y1

x1
.
.
.

xn

(θ)

x1* .
.
xn* .

v1

A*

.
.
.

yn
.
.
.

vn

Figure 1: The separation system G (solid line) and the generation of the regularization term R MSE(dashed line). RMSE = ∑ni=1 v2i , where vi = xi − xi∗ .
With RMSE measuring the level of nonlinear distortion, nonlinear ICA with MND can be formulated as the following constrained optimization problem. It aims to minimize the mutual information
between outputs, that is, I(y), subject to RMSE (θ) ≤ t, where t is a pre-assigned parameter. The Lagrangian for this optimization problem is L(θ, λ) = I(y) + λ[RMSE (θ) − t] with λ ≥ 0. To find θ, weneed to minimize
J = I(y) + λRMSE (θ).
(3)
The non-negative constant λ depends on the pre-assigned parameter t.
Another advantage of the MND principle is that it tends to make the mapping G invertible. In
the general nonlinear ICA problem, it is assumed that both F and G are invertible. But in practice
4. If E(y) = E(x) = 0, x∗ can be obtained as x∗ = A∗ y instead, and here A∗ is a n × n matrix.

2458

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

it is not easy to guarantee the invertibility of the mapping provided by a flexible nonlinear function
approximator, like the MLP. MND pushes G to be close to a linear invertible transformation. Hencewhen nonlinearity in F is not too strong, MND helps to guarantee the invertibility of the nonlinear
ICA separation system G .
2.2.1 S IMPLIFICATION OF RMSE
RMSE , given in Eq. 2, can be further simplified. According to Eq. 2, the derivative of R MSE w.r.t.
∗
T∗TA∗ is ∂R∂AMSE∗ = −2E{(x − A ỹ)ỹ }. Setting the derivative to 0 gives E{(x − A ỹ)ỹ } = 0, which
impliesA∗ = E{xỹT }[E{ỹỹT }]−1 .
(4)
We can see that due to the adoption of the MSE, A∗ can be obtained in closed form. This greatly
simplifies the derivation of learning rules.
Due to Eq. 4, we have E{A∗ ỹỹT A∗T } = E{xỹT }A∗T ,and RMSE then becomes
RMSE = Tr E{(x − A∗ ỹ)(x − A∗ ỹ)T }= Tr E{xxT − A∗ ỹxT − xỹT A∗T + A∗ ỹỹT A∗T }
= Tr E{xxT − A∗ ỹxT − xỹT A∗T + xỹT A∗T }= −Tr E{A∗ ỹxT } + Tr E{xxT }= −Tr E{xỹT }[E{ỹỹT }]−1 E{ỹxT } + const.
(5)
Since yi are independent from each other, they are uncorrelated. We can also easily make y i zeromean. Consequently, E{ỹỹT } = diag{E(y21 ), E(y22 ), ..., E(y2n ), 1}, and RMSE becomes
RMSE = −Tr E{xỹT } · [diag{E(y21 ), ..., E(y2n ), 1}]−1 · E{ỹxT } + constnn

E 2 (x j yi )+ const.
2
j=1 i=1 E(yi )

= −∑∑

(6)

RMSE depends only on the inputs and the outputs of the nonlinear ICA system G (θ). Given a form
for G , the learning rule for nonlinear ICA with MND is derived by minimizing Eq. 3. Note that
RMSE , given in Eq. 2, is inconsistent with certain scaling properties of the observations x. To avoid
this, one needs to normalize the variance of the observations x i through preprocessing, if necessary.
2.2.2 D ETERMINATION OF THE R EGULARIZATION PARAMETER λWe suggest initializing λ with a large value λ0 at the beginning of training and decreasing it to asmall constant λc during the learning process. A large value for λ at the beginning reduces the possibility of getting into unwanted solutions, which may be non-trivial transformations of the original
sources si or local optima. As training goes on, the influence of the regularization term is relaxed,and G gains more freedom. Hopefully, nonlinearity will be introduced, if necessary. The choice of
λc depends on the level of nonlinear distortion in the mixing procedure. If the nonlinear distortionis considerable, we should use a very small value for λc to give the G network enough flexibility. In
our experiments, we found that the separation performance of nonlinear ICA with MND is robust to
the value of λc in a certain range. If the variance of the observations x i is normalized, typical values
used in our experiments are λ0 = 5 and λc = 0.01.
2459

Z HANG AND C HAN

2.3 Relation to Previous Works
The MISEP method has been reported to solve some nonlinear BSS problems successfully, including separating a real-life nonlinear image mixture (Almeida, 2005, 2003). Almeida (2003) claimed
that the MLP itself may provide suitable regularization for nonlinear ICA. Some means were also
used for regularization in the experiments there. For example, first, direct connections between inputs and output units were incorporated in the G network. Direct connections can quickly adapt
the linear part of the mapping G . Second, in Almeida (2005), the G network was initialized with
an identity mapping, and during the first 100 epochs, it was constrained to be linear (by keeping
the output weights of the hidden layer equal to zero). After that, the G network began learning
the nonlinear distortion. G is therefore expected to be not far from linear, and MND is achieved
to some extent. Accordingly, nice experimental results reported there could support the usefulnessof the MND principle. We should mention that the MND principle formulated here, as well as thecorresponding regularizer, provides a way to control the nonlinearity of the mixing mapping. It
can be incorporated by any nonlinear ICA method, including MISEP. Later, we will investigate theeffect of MND on nonlinear ICA solutions theoretically, and compare various related nonlinear ICA
methods empirically.
In the kernel-based nonlinear BSS method (Harmeling et al., 2003), the data are first mapped
to a high-dimensional kernel feature space. Next, a BSS method based on second order temporaldecorrelation is performed. In this way a large number of components are extracted. When the
nonlinearity in data generation is not too strong, the MND principle provides a way to select a
subset of output components corresponding to the original sources. Assume that the outputs y i are
made zero-mean and of unit variance. From Eq. 6 we can see that one can select y i with large
E 2 (x y )∑nj=1 E(yj2 )i = ∑nj=1 E 2 (x j yi ) = ∑nj=1 var(x j ) · corr2 (x j , yi ).
i

It is worth mention that the principle of least mean square error reconstruction has been used fortraining a class of neural networks and gives some interesting results (Xu, 1993). For one-layer networks with linear/nonlinear units, this principle leads to principal component analysis (PCA)/ICA.
We should address that the reconstruction in their work is quite different from that discussed in
Section 2.2 in this paper. In their work, the forward process and the reconstruction process sharethe same weights; in this paper, reconstructed signals are an affine mapping of the outputs, andparameters in the affine mapping are determined by minimizing the reconstruction error.
Smoothness provides a constraint to prevent a neural network from overfitting noisy data. It
is also useful to ensure nonlinear ICA to result in nonlinear BSS (Almeida, 2003). In fact, thesmoothness regularizer exploiting second-order derivatives (Tikhonov and Arsenin, 1977; Poggio
et al., 1985) is also related to the MND principle, as shown below.
2.4 Local Minimal Nonlinear Distortion: Smoothness
RMSE , given in Eq. 2, indicates the deviation of the mapping F̂ from the affine mapping which fitsF̂ globally best. In contrast, one may enforce the local MND of the nonlinear mapping averaged atevery point. We will show that this regularization actually leads to the smoothness regularizer exploiting second-order partial derivatives (Tikhonov and Arsenin, 1977; Poggio et al., 1985; Bishop,1993).
2460

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

For a one-dimensional sufficiently smooth function g(x), we can use the second-order Taylor
expansion to approximate its function value in the vicinity of x in terms of g(x):g(x + ε) ≈ g(x) + ∂g T∂x1· ε + εT Hx ε,2where ε is a small variation of x and Hx denotes the Hessian matrix of g. Let 5i j = ∂x∂i ∂xg j . If we use
the first-order Taylor expansion of g, which is a linear function, to approximate g(x + ε), the squareerror is2 ∂g T
2
1 T
1 n
2
5
ε
ε
·ε ≈
ε Hx ε =g(x + ε) − g(x) −
i
j
i
j
∂x
4
4 i,∑
j=1






n
n
n
n
n
2
1
1
1
ε2i ε2j =52i j
52i j ∑ ε2i = ||ε||4 · ∑ 52i j .
≤∑∑∑4 i, j=1
4 i, j=1
4i, j=1i=1i, j=1
2

The above inequality holds due to the Cauchy’s inequality. Now we can see that in order to achieve
the local MND of g averaged in the domain of x, we just need to minimize the followingZn∑ 52i j dx =Dx i, j=1

Z

Dx


n
2
2
5
5
+
2
∑ ii ∑ i j dx.

 n

i=1

(7)

i, j=1,i< jThis regularizer has been used for achieving the smoothness constraint (see, e.g., Grimson 1982 for
its application in computer vision). When the mapping is vector-valued, we need to apply the above
regularizer to each component of the mapping.
Originally we intended to do regularization on the mixing mapping F̂ , but it is difficult to2xldo since it is hard to evaluate ∂y∂i ∂y. Instead, we do regularization on G , the inverse of F̂ . Thejregularization term in Eq. 3 then becomes
Rlocal (θ) =

Z

n

n

n

 ∂2 y 2∑ ∑ ∑ ∂xi ∂x jD
x l=1 i=1 j=1

l

dx =

Z

n

n

∑ ∑ Pi j dx,D(8)x i=1 j=1

 2 2yl. Nonlinear ICA with a smooth de-mixing mapping can be achieved bywhere Pi j , ∑nl=1 ∂x∂i ∂xjminimizing the mutual information between yi , with Rlocal , given by Eq. 8, as the regularization
 2 22∂ ylterm. There are totally n (n+1)
different
terms
in the integrand of Rlocal . For simplic2∂xi ∂x j
ity and computational reasons, sometimes one may drop the cross derivatives in Eq. 8, that is, 2 2∂ yl
with i 6= j, and consequently obtain the curvature-driven smoothing regularizer proposed∂xi ∂x j
in Bishop (1993), with the number of different terms in the integrand being n 2 .

3. Incorporation of MND in Different Nonlinear ICA Methods
Now we should choose a model for the nonlinear ICA separation system G (θ) and give the learning
rule for nonlinear ICA with MND as well as nonlinear ICA with the smoothness constraint for G .
Two nonlinear ICA methods are considered here. They are MISEP (Almeida, 2003) and nonlinear
ICA based on kernels (Zhang and Chan, 2007a).
2461

Z HANG AND C HAN

3.1 MISEP with MND
Before incorporating MND into the MISEP method (Almeida, 2003) for nonlinear ICA, we give an
overview of this method.
3.1.1 MISEP FOR N ONLINEAR ICAMISEP adopts the MLP to model the separation function G in the nonlinear ICA problem. Figure 2
shows the structure used in this method. This method extends the original Infomax method forlinear ICA (Bell and Sejnowski, 1995) in two aspects. First, the separation system is a nonlinear
transformation, which is modeled by the MLP. Second, the nonlinearities ψ i are not fixed in advance,
but tuned by the Infomax principle, together with G .

x1
xn

y1
.
.
.

.
.
.

ψ1

yn

.
.
.

ψn

u1
un

Figure 2: The network structure used in Infomax and MISEP. G is the separation system, and ψ iare the nonlinearities applied to the separated signals. In MISEP, G is a nonlinear transformation, and both G and ψi are learned by the Infomax principle.
With the Infomax principle, parameters in G and ψi are learned by maximizing the joint entropyof the outputs of the structure in Figure 2, which can be written as H(u) = H(x) + E{log | det J|},
where J = ∂u∂x is the Jacobian of the nonlinear transformation from x to u. As H(x) does not depend on the parameters in G and ψi , it can be considered as a constant. Maximizing H(u) is thus
equivalent to minimizing
J1 (θ) = −E{log | det J|},

(9)where θ denotes the set of unknown parameters. The learning rules for θ were derived by Almeida(2003), in a manner similar to the back-propagation algorithm.
The MLP adopted in this paper has linear output units and a single hidden layer. For the hidden
units, the activation function l(·) may be the logistic sigmoid function, the arctan function, etc.
Direct connections between the inputs and output units are also allowed. Let a = [a 1 , ..., aM ]T be
the inputs to the hidden units, z = [z1 , ..., zM ]T be the output of the hidden units, and W and b
denote the weights and biases, respectively. We use superscripts to distinguish the locations of theseparameters: W(d) denotes the weights from the inputs to output units, W(1) those from the inputs
to the hidden layer, and W(2) those from the hidden layer to the output units. b(1) and b(2) are thebias vectors in the hidden layer and in the output units, respectively. The output of the G network
represented by this MLP takes the form:
y = W(2) · z + W(d) x + b(2) , where
(1)

(1)

zi = l(ai ), and a = W x + b .
2462

(10)

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

3.1.2 MISEP WITH MND
For MISEP with MND, the objective function to be minimized is Eq. 9 regularized by R MSE given
in Eq. 6. The learning rule for θ to minimize Eq. 9 has been considered in Almeida (2003). Hencehere we only give the gradient of RMSE w.r.t. θ.
Using the chain rule, also noting Eq. 10, the gradient of R MSE (θ) w.r.t. W(2) can be obtained:n n h E 2 (x y )E(x j yi ) i ∂yi o
∂RMSE
j i
2
y
−
x ·
=
E
= E{K · zT },(11)i∑ E 2 (y2 )
2) i
(2)
E(y
∂W(2)
∂W
i
i
i=1
i
h 2
E(x y )E (x y )where K , [K1 , ..., Kn ]T with its i-th element being Ki = 2 ∑ j E 2 (yj2 )i yi − E(yj2 )i x j , and z = [z1 , z2 ,ii..., zM ]T is the output of the hidden layer of the MLP. For the gradient of R MSE w.r.t. W(1) , W(d) ,b(2) , and b(1) , see Appendix A.
3.1.3 MISEP WITH S MOOTHNESS C ONSTRAINT ON G
The mapping provided by a MLP may not be smooth enough to make nonlinear ICA result in
nonlinear BSS. So here we also implement MISEP with the smoothness constraint on G . The
objective function to be minimized becomes Eq. 9 regularized by R local given in Eq. 8. Pi j appears
in the expression of Rlocal . We first derive its gradient w.r.t. θ in a way analogous to that in Bishop(1993); see Appendix B.
local, the integral in Eq. 8 is difficult to evaluate. Below are two ways to tackle
In calculation of ∂R∂θthis problem. A very simple way to approximate Eq. 7 is to use the average of the integrand over allobservations instead of the integral (ignoring a constant scaling factor), just as Bishop (1993) does:
n n n
o
(1)
Rlocal (θ) = E ∑ ∑ Pi j .
(12)
i=1 j=1

This approximation actually assumes that the distribution of x is close to uniform, as seen frombelow. Eq. 8 can be rewritten asZon 1 n n
1 n n
(13)
Rlocal (θ) =p(x) ·Pdx=
EP∑ ∑ ij∑ ∑ ij .
p(x) i=1p(x) i=1
Dx
j=1
j=1
If p(x) is a constant in the domain Dx , Eq. 12 is equivalent to Eq. 13; otherwise, the approximation
using Eq. 12 may result in large error, and we may need another way to approximate the integral in
Eq. 8.
When the nonlinear ICA algorithm has run for a certain number of epochs, u, the output of
the system in Figure 2, has approximately independent components and is approximately uniformly distributed in [0, 1]n . This means that p(u) is approximately 1. As p(x) = p(u) · | det J|,
one
is approximately equal to | det J|. Consequently Eq. 13 becomes R local (θ) ≈n can see that p(x) onn1
E | det J| ∑i=1 ∑ j=1 Pi j . The gradient of Rlocal (θ) is
n 1
n n∂Pi j o∂Rlocal (θ)
≈E
.
∑
∑
∂θ
| det J| i=1 j=1 ∂θ

(14)As J = ∂u∂x can be easily calculated according to the network structure in Figure 2, Eq. 14 is alsoeasy to evaluate, using Eq. 20 of Appendix B.
2463

Z HANG AND C HAN

3.2 MND for Nonlinear ICA Based on Kernels
Nonlinear ICA based on kernels (Zhang and Chan, 2007a) exploits kernel methods to construct the
separation system G , and unknown parameters are adjusted by minimizing the mutual information
between outputs yi .5 We have applied the MND principle and the smoothness regularizer to nonlinear ICA based on kernels; for details, see Zhang and Chan (2007a). Note that unlike the mapping
provided by a MLP, which is comparatively smooth, the mapping constructed by kernel methods
may not be smooth. So it is quite necessary to explicitly enforce the smoothness constraint for
nonlinear ICA based on kernels.

4. Investigation of the Effect of MND
In this section we intend to explain why the MND principle, including the smoothness regularization, helps to alleviate the ill-posedness of nonlinear ICA from a mathematical viewpoint. There
are two types of indeterminacies in solutions to nonlinear ICA, namely trivial indeterminacies andnon-trivial indeterminacies. Trivial indeterminacies mean that the estimate of s j produced by nonlinear ICA may be any nonlinear function of s j ; non-trivial indeterminacies mean that the outputs
of nonlinear ICA, although mutually independent, are still a mixing of the original sources. Let us
begin with the effect of MND on trivial indeterminacies.
4.1 For Trivial Indeterminacies
Let us assume in this section that, in the solutions of nonlinear ICA, each component depends only
on one of the sources. Before presenting the main result, let us first give the following lemma.
Lemma 1 Suppose that we are given the random vector d = (d 1 , d2 , · · · , dn )T . Let Ry be the meansquare error of reconstructing d from the variable y with the best-fitting linear transformation, that
is, Ry = mina E{||d − a · y||2 }, where a = (a1 , a2 , · · · , an )T . The variable y which gives the minimum
Ry is the first non-centered principal component of d multiplied by a constant, and if y is constrained
to be zero-mean, it is the first principal component of d multiplied by a constant.
See Appendix C for a proof. Now let us consider a particular kind of nonlinear mixtures, inwhich each observed nonlinear mixture xi is assumed to be generated by
xi = fi1 (s1 ) + fi2 (s2 ) + · · · + fin (sn ),(15)

where fi j are invertible functions. We call such nonlinear mixtures distorted source (DS) mixtures,
since each observation is a linear mixture of nonlinearly distorted sources. For this nonlinear mixing model, we have the following theorem on the effect of MND on trivial indeterminacies in its
nonlinear ICA solutions. Here the following assumptions are made:A1: In the output of nonlinear ICA, each component depends only on one of the sources and iszero-mean.
A2: The nonlinear ICA system has enough flexibility to reach the minimum of the MND regularization term RMSE defined by Eq. 2.
5. The difference between nonlinear ICA based on kernels discussed here and the kernel-based nonlinear BSS method
by Harmeling et al. (2003) should be made clear. Both of them use kernels. However, the former produces statisticallyindependent outputs, while the latter exploits the temporal structure of the sources for separation.

2464

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

Theorem 1 Suppose that each observed nonlinear mixture x i is generated according to Eq. 15.
Under assumptions A1 & A2, the estimate of s j produced by nonlinear ICA with MND is the firstprincipal component of f∗ j (s j ) = [ f1 j (s j ), · · · , fn j (s j )]T , multiplied by a constant.
See Appendix D for a proof. The DS mixing model Eq. 15 may be restrictive. Now let us
consider the case where nonlinearity in F is mild such that F can be well approximated by its
Maclaurin expansion of degree 3. Let
Oi, j =

∂2 x i∂xi∂3 x i, Oi, jk =, and Oi, jkl =.
∂s j s j =0
∂s j ∂sk s j ,sk =0∂s j ∂sk ∂sl s j ,sk ,sl =0

The following theorem discusses the effect of MND on trivial indeterminacies of nonlinear ICA
solutions in this case. In particular, it states that by incorporating MND into the nonlinear ICA
system, trivial indeterminacies in nonlinear ICA solutions are overcome; it shows how the outputsof the nonlinear ICA system, as the estimate of the sources, are related to the original sources s i andthe mixing system F .
Theorem 2 Suppose that each component of the mixing mapping F = ( f 1 , · · · , fn )T in Eq. 1 isgenerated by the following Maclaurin series of degree 3:
xi = fi (s) = fi (0) + ∑ Oi, j s j +j1
1
Oi, jk · s j sk + ∑ Oi, jkl · s j sk sl ,∑2 j,k
6 j,k,lwhere E{s j } = 0 and E{s2j } = 1, for j = 1, · · · , n. Let


1
1
1
Di j (s j ) , Oi, j + ∑ Oi, jkk · s j + Oi, j j · s2j + Oi, j j j · s3j .
2 k6= j2
6
And let D̃i j (s j ) be the centered version of Di (s j ), that is, D̃i j (s j ) = Di j (s j ) − Ei {Di j (s j )}. Underassumptions A1 & A2, the estimate of s j produced by nonlinear ICA with MND is the first principal
component of D̃∗ j (s j ) = [D̃1 j (s j ), · · · , D̃n j (s j ), ]T , multiplied by a constant.
See Appendix E for a proof. Under the condition that nonlinear distortion in the mixing mapping

F is not strong, D̃i j (s j ) would not be far from linear. Moreover, if the nonlinear part of D̃i j (s j )varies for different i, the estimate of s j is expected to be closer to linear than D̃i j (s j ), because it is
the first principal component (PC) of D̃∗ j (s j ). To summarize, Theorems 1 and 2 show that trivial
indeterminacies in nonlinear ICA solutions can be overcome by the MND principle; and whenthe mixing mapping is not strong, the nonlinear distortion in the nonlinear ICA outputs w.r.t. theoriginal sources is weak.
4.1.1 R EMARKIn the proof of Theorems 1 and 2, we have made use of the fact that mutual information is invariant
to any component-wise strictly monotonic nonlinear transformation of the variables. Consequently,
trivial transformations do not affect the first term in Eq. 3, and they can be determined by minimizing
RMSE only, as claimed in the theorems. However, in practical implementations of nonlinear ICA
algorithms, one needs to estimate the densities of yi or their variations. Due to estimation error, the2465

Z HANG AND C HAN

gradient of the mutual information I(y1 , · · · , yn ) may be sensitive to the distribution of yi , or it may
be slightly affected by trivial transformations. This may cause the results of Theorems 1 and 2 to beviolated slightly.
Fortunately, this phenomenon can be avoided easily. To model the trivial transformations, weapply a separate nonlinear function approximator (such as a MLP) to each output of nonlinear
ICA to generate the final nonlinear ICA result. These nonlinear function approximators are thenlearned by minimizing RMSE (Eq. 6). This provides a way to tackle the trivial indeterminacies; after
performing nonlinear ICA with any nonlinear ICA method, if we know that there only exist trivial
indeterminacies, we can adopt the above technique to determine the trivial transformations.
4.2 For Non-Trivial Indeterminacies
Now let us investigate the effect of MND on non-trivial indeterminacies in nonlinear ICA solutions.
Generally speaking, there exist an infinite number of ways in which non-trivial indeterminacies
occur, and it is impossible to formulate all of them. Hyvärinen and Pajunen (1999) gave somefamilies of non-trivial transformations preserving mutual independence.
4.2.1 A PARTICULAR C LASS OF N ON -T RIVIAL I NDETERMINACIESFor the convenience of analysis, here we consider the following manner to construct non-trivial
transformations preserving mutual independence. First, using the Gaussianization technique (Chenand Gopinath, 2001), we transform each of the independent variables s i to a standard Gaussianvariable ui with an strictly increasing function qi , that is, ui = qi (si ). Clearly ui are mutually independent. Second, we can apply an orthogonal transformation U to u = (u 1 , · · · , un )T . The components of e = Uu are still jointly Gaussian and mutually independent. 6 Finally, let y = r(e), wherer = (r1 , · · · , rn )T is a component-wise function with each ri strictly increasing. Components of y arestill mutually independent. That is, y is always a solution to nonlinear ICA of the nonlinear mixture
x = F (s). The procedure transforming s to y can be described as r ◦ U ◦ q, as shown in Figure 3.
When U is a permutation matrix, this transformation is trivial; otherwise it is not.
.
.
.

.
.
.

Figure 3: A non-trivial transformation from s to y preserving independence, that is, r ◦ U ◦ q.
4.2.2 E FFECT OF MND
To see the effect of MND on y in Figure 3 (recall that y is a solution to nonlinear ICA of x =F (s)), we need to find how MND affects U, as well as ri . First, let us consider the case where
the outputs yi are Gaussian, meaning that each component of r̃ is a linear mapping. Without lossof generality, we further assume that yi are zero-mean and of unit variance, that is, E{yyT } = I.
Consequently, ri are identity mappings and y =e = Uu. Assuming xi are zero-mean, according
to Eq. 5, We have RMSE
= −Tr E{xyT }E{yxT } + const = −Tr E{xuT }UT UE{uxT } + const =−Tr E{xuT }E{uxT } + const. In this case RMSE is not affected by U, and the MND principle
6. U may depend on ||u||. In other words, U may be different for u of different norms.

2466

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

could not help to avoid such non-trivial indeterminacies. We have empirically found that in general,
when yi are close to Gaussian, the separation performance tends to be bad. To make sure thatthe separation result is reliable, one should check the non-Gaussianity of y i after the algorithm
converges.
Next, suppose that that both sl and y j are non-Gaussian. ri are then nonlinear. Consider theextreme case that the mixing mapping F is linear; in order to minimize R MSE (Eq. 2), U in Figure 3
must be a permutation matrix. One can then image that if nonlinearity in F is weak enough, U in
Figure 3 should be approximately a permutation matrix, meaning that the original sources s could
be recovered.
However, if nonlinearity in F is strong, U may not be a permutation matrix, and non-trivial
transformations from s to y may occur. This is actually quite natural. Consider the mixing mappingx = F (s) which can be decomposed as a non-trivial transformation of s shown in Figure 3 (denoteby z its output), followed by a nonlinear transformation x = F L (z) which is close enough to linear. In
this situation, the output of nonlinear ICA with MND would be an estimate of z, and if no additional
knowledge of the mixing mapping is given, it is impossible to recover the original sources s i .
Below we give an two-channel example to illustrate the relationship between R MSE and theorthogonal matrix U when nonlinearity in F is strong. The two independent sources are a uniformly distributed signal and a super-Gaussian signal, and their scatter plot is given in Figure 9(a).
The observations xi , whose scatter plot is shown in Figure 4(a), are generated by applying a 23-2 MLP to the source signals. From this figure we can see that nonlinearity in the mixing procedure is comparatively strong. The orthogonal matrix U in Figure 3 is parameterized as U =[cos(α), − sin(α); sin(α), cos(α)]. From Eq. 6 and Figure 3, one can see that R MSE depends on α
and ri . For each value of α, ri (i = 1, 2) are modelled by a 1-6-1 MLP and they are learned byminimizing RMSE . Finally, minri RMSE is a function of α, with a period of 90 degrees, as plottedin Figure 4(b). In this example, α determined by the MND principle is about 11 degrees. It is notthat close to zero, but it is still comparatively small and consequently the sources s i are recovered
approximately.
(a)

(b)

5

−1.58

Best−fitting linear mixture
Nonlinear mixture

−1.6
−1.62

min R

0

−1.66

ri

x2

MSE

−1.64

−1.68
−1.7
−1.72
−1.74

−5
−5

0

x1

−1.76
−90 −70 −50 −30 −10 0 10

5

α (degree)

30

50

70

90

Figure 4: (a) Nonlinear mixtures of a sinusoid source signal and a super-Gaussian source signal(whose scatter plot is given in Figure 9.a) generated by a 2-3-2 MLP. x-mark points show
linear mixtures of the sources which fit the nonlinear mixtures best. (b) min ri RMSE as afunction of α, whose minimum is achieved at α ≈ 11 degrees.
2467

Z HANG AND C HAN

5. Simulations
In this section we investigate the performance of the proposed principle for solving nonlinear ICA
using synthetic data. The experiments in Zhang and Chan (2007a) have empirically shown thatboth MND and the smoothness constraint are useful to ensure nonlinear ICA based on kernels toresult in nonlinear BSS, when nonlinear distortion in the mixing procedure is not very strong. As
its performance depends somewhat crucially on the choice of the kernel function, nonlinear ICA
based on kernels is not used for comparison here. The following six methods (schemes) were used
to separate various nonlinear mixtures:
1. MISEP: The MISEP method (Almeida, 2003) with parameters θ randomly initialized. 7 Notethat in this method, the smoothness constraint has been implicitly incorporated to some extent,
due to the property of the adopted MLP.
2. Linear init.: The MISEP method with G initialized as a linear mapping. This was achieved byadopting the regularization term Eq. 2 with λ = 5 (which is very large) in the first 50 epochs.
3. MND: The MISEP method incorporating MND, with RMSE , the mean square error of the best
linear reconstruction, as the regularization term (Section 2.2). The regularization parameter λdecayed from λ0 = 5 to λc = 0.01 in the first 350 epochs. After that λ was fixed as λc .
4. Smooth (I): The MISEP method with the smoothness regularizer (Section 2.4) explicitly incorporated. λ decayed from 1 to 0.004 in the first 350 epochs.
5. Smooth (II): Same as Smooth (I), but λ was fixed to 0.007.
6. VB-NICA: Bayesian variational nonlinear ICA (Lappalainen and Honkela, 2000; Valpola,
2000).8 PCA was used for initialization. After obtaining nonlinear factor analysis solutions
using the package, we applied linear ICA (FastICA by Hyvärinen 1999 was used) to achieve
nonlinear BSS.
In addition, in order to show the necessity of nonlinear ICA methods for separating nonlinear mixtures, linear ICA (FastICA was adopted) was also used to separate the nonlinear mixtures.
It was addressed in Section 2.3 that the incorporation of direct connections between inputs andoutput units in the MLP representing G implicitly and roughly implements the MND principle. To
check that, in our experiments, the MLP without direct connections and that with direct connections
were both adopted to represent G , for comparison reasons. Like in Almeida (2003), the MLP has
20 arctan hidden units, 10 of which are connected to each of the output units of G .
We use the signal to noise ratio (SNR) of yi relative to si , denoted by SNR(yi ), to assess theseparation performance of si . Besides, we apply a flexible nonlinear transformation h to y i to minimize the MSE between h(yi ) and si , and use the SNR of h(yi ) relative to si as another performancemeasure. In this way possible trivial transformations between s i and yi are eliminated. In our experiments h was implemented by a two-layer MLP with eight hidden units with the hyperbolic tangent
activation function and a linear output unit. This MLP was trained using the MATLAB neuralnetwork toolbox.
7. Source code is available at http : //www.lx.it.pt/ ∼ lbalmeida/ica/mitoolbox.html.
8. Source code is available at http : //www.cis.hut.fi/projects/bayes/. The following MATLAT commands were used
to produce the ouput y: [nlfa sources, net, params, status, fs] = nlfa(x, ’searchsources’, 2, ’hidneurons’, 15, ’iters’,2000); y = fastica(nlfa sources.e, ’approach’, ’symm’, ’g’, ’tanh’);2468

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

Three kinds of nonlinear mixtures were investigated. They are distorted source (DS) mixtures,
post-nonlinear (PNL) mixtures, and generic nonlinear (GN) mixtures which are generated by a MLP.
Both super-Gaussian and sub-Gaussian sources were used.
5.1 For Distorted Source Mixtures
We first considered the DS mixtures defined in Eq. 15. Specifically, in the experiments the twochannel mixtures xi were generated according to x1 = a11 s1 + f12 (s2 ), x2 = f21 (s1 ) + a22 s2 , where
a11 = a22 = 1, and f 12 (si ) = f21 (si ) = 3 tanh(si /4) + 0.1si . We used two super-Gaussian source
signals, which are generated by si = 35 ni + 25 n3i , where ni are independent Gaussian signals. Eachsignal has 1000 samples. Figure 5 shows the scatter plot of the sources s i and that of the observations
xi . To see the level of nonlinear distortion in the mixing transformation, we also give the scatter plotof the affine transformation of si which fits xi the best.
(a)

(b)
8

15

6

10

4

5

2

Best−fitting linear mixture
Nonlinear mixture

x2

s2

20

0

0

−5

−2

−10

−4

−15
−20

−10

0

10

−6
−6

20

s1

−4

−2

0

2

4

6

x1

Figure 5: (a) Scatter plot of the sources si generating the DS mixtures. (b) Scatter plot of the DS
mixtures xi . x-mark points are linear mixtures of si which fit xi best.

To reduce the random effect, all methods were repeated for 40 runs, and in each run the MLP
was randomly initialized. We found that the separated results in the two channels have a similar
SNR, so for saving space, here we just give the SNR in the first channel. Figure 6 compares theboxplot of SNR(y1 ) and SNR(h(y1 )) for different methods. In Figure 6 (a, b), the MLP has nodirect connections between inputs and output units, while in (c, d) the MLP has direct connections.
We can see that in this case the methods MND, Smooth(I), and Smooth(II) give very high SNR,
and at the same time, produce fewest unwanted results. Moreover, the MLP with direct connections
behaves better than that without direct connections. The performance of VB-NICA is not very
good. The reason may be that this method does not take into account the very useful informationthat nonlinearity in the mixing mapping is not very strong. It should be noted that VB-NICA may
not exhibit its potential powerfulness in the experiments, since the source number is given and nonoise is considered.
2469

Z HANG AND C HAN

(a)

(b)
(dB )

(dB )

18
16

20

14

SNR(h(y1))

1

SNR(y )

12
10
8

15

10

FastICA

FastICA

6
4

5

2
0

0

MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICAMISEP Linear init.

MND

Method

(c)

(d)
30

(dB )

(dB )

20
18
16
14

25

20

12

SNR(h(y1))

SNR(y1)

Smooth(I) Smooth(II) VB−NICA
Method10
8

FastICA

15

10

6
4

FastICA

5

2
0

MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA
Method0

MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA
MethodFigure 6: Boxplot of the SNR of separating the DS mixtures by the MLP without or with directconnections between inputs and output units. Top: Without direct connections. Bottom:
With direct connections. (a, c) SNR(y1 ). (b, d) SNR(h(y1 )).

5.2 For Post-Nonlinear Mixtures
The second experiment is to separate PNL mixtures. We used two sub-Gaussian source signals,
which are a uniformly distributed white signal and a sinusoid waveform. The sources were firstmixed with the mixing matrix A = [-0.2261, -0.1189; -0.1706, -0.2836], producing linear mixturesz. The observations were then generated as x1 = z1 /2.5 + tanh(3z1 ) and x2 = z2 + z32 /1.5. Figure 7
shows the scatter plot of the sources and that of the PNL mixtures (after standardization). Figure 8
gives the separation performance of s1 by various methods.9 In this case, the proposed nonlinear
ICA with MND (labelled by MND) also gives almost the best results; especially for the MLP without direct connections, the result of nonlinear ICA with MND is clearly the best. Again, the MLP
with direct connections produces better results. Moreover, one can see that compared to the DS
9. If we use the PNL mixing model (Taleb and Jutten, 1999) to separate such mixtures, theoretically the sources could
be well recovered. But in this paper we assume that the form of the mixing procedure is unknown, and treat it as a
general nonlinear ICA problem.

2470

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

mixtures in Section 5.1, the PNL mixtures considered here are comparatively hard to be separated
by the MLP structure.
(a)

(b)

1.5

Best−fitting linear mixture
Nonlinear mixture

2.5
2

1

1.5
0.5

1

x2

s2

0.5
0

0
−0.5

−0.5

−1
−1.5

−1

−2
−1.5
−2

−1

0

1

−2.5

2

s1

−2

−1

0

1

2

x1

Figure 7: (a) Scatter plot of the sources si generating PNL mixtures. (b) Scatter plot of the PNL
mixtures xi .

5.3 For Generic Nonlinear Mixtures
We used a 2-2-2 MLP to generate nonlinear mixtures from sources. Hidden units have the arctan
activation function. The weights between the input layer and the hidden layer are random numbers
between -1 and 1. They are not large such that the mixing mapping is invertible and the nonlinear
distortion produced by the MLP would not be very strong. The sources used here were the firstsource in Experiment 1 (super-Gaussian) and the second one in Experiment 2 (sub-Gaussian). Figure 9 shows the scatter plot of the sources and that of the GN mixtures. The performance of various
methods for separating such mixtures is given in Figures 10. Apparently nonlinear ICA with MND
gives the best separation results in this case.
Summed over all the three cases discussed above, we can see that MISEP with MND produces
promising results for the general nonlinear ICA problem, provided that nonlinearity in the mixing
mapping is not very strong. Specifically, it gives the fewest unwanted solutions, and its separation
performance is very good. Moreover, the MLP with direct connections usually performs better thanthat without direct connections, but we also found that in some cases it got stuck into unwanted
solutions more easily.
5.4 On Trivial Indeterminacies
In Section 4.1 we have discussed the effect of the MND principle on trivial indeterminacies ofnonlinear ICA solutions. In particular, Theorem 1 states that for DS mixtures, if there are only
trivial indeterminacies, each output of nonlinear ICA with MND is the PC of the contributions of
the corresponding source to all mixtures. Now let us illustrate this with the help of the DS mixtures
used in Section 5.1.
2471

Z HANG AND C HAN

(a)

(b)
(dB )

(dB )

22
20

20

18
16

SNR(h(y1))

SNR(y1)

14
12
10

10

FastICA

8

15

FastICA

6
5

4
2
MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA
MethodMISEP Linear init.

(dB )

20
18

22
20

16

18

14

16

12
10

FastICA

8

Smooth(I) Smooth(II) VB−NICA
Method(d)

SNR(h(y1))

SNR(y1)

(dB )

(c)

MND

14
12
10

FastICA

8
6
6
4

4

2

2
MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA

MethodMISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA
MethodFigure 8: Boxplot of the SNR of separating the PNL mixtures by the MLP without or with directconnections between inputs and output units. Top: Without direct connections. Bottom:
With direct connections. (a, c) SNR(y1 ). (b, d) SNR(h(y1 )).

Figure 11 shows the relationship between yi obtained by MISEP with MND in one run and the
PC of f∗i (si ) = [ f1i (si ), f2i (si )]T . We can see that each yi is actually not very close to the corresponding PC, which may be caused by two reasons. First, there may exist some weak non-trivial
transformation in the solution, as seen from the points close to the origin in Figure 11(b); y 2 is not
solely dependent on s2 , but also slightly affected by s1 . The other reason is the error in estimating
the density of yi or its variation involved in the MISEP method, as explained in Section 4.1.1. We
use the method proposed there to avoid the effect of the estimation error: a 1-8-1 MLP, denoted by
τi , is applied to each yi , and τi (yi ) is taken as the final nonlinear ICA output. Each τi is learned byminimizing RMSE (Eq. 6). The resulting τi (yi ) is almost identical to the corresponding PC of f∗i (si ),as seen from Figure 12. This has confirmed Theorem 1 and the validity of the method for tackling
trivial indeterminacies proposed in Section 4.1.1.
2472

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

(a)

(b)
6

8

Best−fitting linear mixture
Generic nonlinear mixture

6

4
4

2

x2

s2

2
0

0

−2

−2
−4

−4
−6
−8
−2

−1

0
s1

1

−6
−5

2

0

x

5

1

Figure 9: (a) Scatter plot of the sources si . (b) Scatter plot of the GN mixtures xi .

6. Application to Causality Discovery in the Hong Kong Stock Market
In this section we give a real-life application of nonlinear ICA with MND. Specifically, we use thismethod to discover linear causal relations among the daily returns of a set of stocks. The empiricalresults were ever reported in Zhang and Chan (2006), without much detail of the method.
6.1 Introduction
It is well known that financial assets are not independent of each other, and that there may be somerelations among them. Such relations can be described in different ways. In risk management,
correlations are used to describe them and help to construct portfolios. The business group, which
is a collection of firms bound together in some formal and/or informal ways, focuses on ties between
financial assets and has attracted a lot of interest (Khanna and Rivkin, 2006). But these descriptions
cannot tell us the causal relations among the financial assets.
The return of a particular stock may be influenced by those of other stocks, for many reasons,such as the ownership relations and financial interlinkages (Khanna and Rivkin, 2006). According to
the efficient market hypothesis, such influence should be reflected in the stock returns immediately.
In this part we aim to discover the causal relations among selected stocks by analyzing their daily
returns.10
Traditionally, causality discovery algorithms for continuous variables usually assume that thedependencies are of a linear form and that the variables are Gaussian distributed (Pearl, 2000).
Under the Gaussianity assumption, only the correlation structure of variables is considered and all
higher-order information is neglected. As a consequence, one obtains some possible causal dia10. In other words, here we aim to find the “instantaneous” causality in the stock market. In contrast, Granger causality(Granger, 1980) analysis has become an important tool to find the “lagged” causality between time series. A timeseries x1 “Granger causes” another series x2 if by incorporating the past history of x1 can improve a prediction of x2
over a prediction based only on the history of x2 alone. The efficient market hypothesis implies no significant Granger
causality between stock returns. In fact, we have applied the approach by Reale and Tunnicliffe Wilson (2001) andpartial directed coherence (Baccala and Sameshima, 2001) to find the Granger causality among the selected stocks,
and very few Granger causal relations were found.

2473

Z HANG AND C HAN

(a)

(b)
(dB )

20

16

18

14

16

12

14

SNR(h(y1))

SNR(y1)

(dB )

18

10
8

FastICA

6

12
10

FastICA

8
6

4
4
2

2

0

0
MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICAMISEP Linear init.

MND

Method(c)

Smooth(I) Smooth(II) VB−NICA
Method

(d)
(dB )

(dB )

18
16

18
16

14
14

SNR(h(y1))

SNR(y1)

12
10
8

12
10

FastICA

8

FastICA
6

6

4

4

2

2
MISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA

MethodMISEP Linear init.

MND

Smooth(I) Smooth(II) VB−NICA
MethodFigure 10: Boxplot of the SNR of separating the GN mixtures by the MLP without or with directconnections between inputs and output units. Top: Without direct connections. Bottom:
With direct connections. (a, c) SNR(y1 ). (b, d) SNR(h(y1 )).

grams which are equivalent in their correlation structure, and cannot find the true causal directions.
Recently, it has been shown that the non-Gaussianity distribution of the variables allows us to distinguish the explanatory variable from the response variable, and consequently, to identify the fullcausal model (Dodge and Rousson, 2001; Shimizu et al., 2006).
In particular, in Shimizu et al. (2006) an elegant and efficient method was proposed for identifying the linear, non-Gaussian, acyclic causal model (abbreviated LiNGAM) by exploiting ICA. If
the data are generated according to the LiNGAM model, theoretically, the ICA de-mixing matrix
W can be permuted to lower triangularity. However, in practice, this may not hold, due to the finite
sample effect, the existence of unobserved confounder variables (Pearl, 2000), or mild nonlinearity and noise that are often encountered in the data generation procedure. To tackle possible mild
nonlinearity in the data generation procedure, we use nonlinear ICA with MND, instead of linear
ICA, to separate the observed data. As the nonlinear distortion is mild, it can be neglected andconsequently, linear causal relations among the observed data can be discovered.
2474

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

(a)

(b)

15

10
8

10

6
4

5
y

2

y1

2

0

0
−2
−4

−5

−6

−10
−10

−10
−10

−8

−5
0
5
Principal component of f (s )
i1

10

1

−5
0
5
Principal component of f (s )
i2

10

2

Figure 11: (a) y1 recovered by MISEP with MND versus the PC of the contributions of s 1 to the DS
mixtures used in Section 5.1. The SNR of y1 w.r.t. the PC of the contributions of s1 is
13.48dB. The dashed line is the linear function fitting the points best. (b) y 2 versus thePC of the contributions of s2 to the DS mixtures. The SNR is 9.12dB.
(a)

(b)

15

10
8
6

10

2
2

1

2

τ (y )

51

τ (y )

4

0
−2

0

−4
−6

−5

−8
−5

0
5
Principal component of f (s )
i1

10

1

−10
−10

−5
0
5
Principal component of f (s )
i2

10

2

Figure 12: (a) τ1 (y1 ) versus the PC of the contributions of s1 to xi . τ1 is modelled by a 1-8-1 MLP
and is learned by minimizing RMSE (Eq. 6). The SNR is 20.99dB. (b) τ2 (y2 ) versus thePC of the contributions of s2 to xi . The SNR is 18.64dB.
6.2 Causality Discovery by ICA: Basic Idea
The LiNGAM model assumes that the generation procedure of the observed data follows the following properties (Shimizu et al., 2006). 1. It is recursive. This is, the observed variables x i , i = 1, ..., n,can be arranged in a causal order, such that no later variable causes any earlier variable. This causal
order is denoted by k(i). 2. The value of xi is a linear function of the values assigned to the earlier
variables, plus a disturbance term ei and an optional constant ci : xi = ∑k( j)<k(i) bi j x j + ei + ci . 3.
2475

Z HANG AND C HAN

ei are independent continuous-valued variables with non-Gaussian distributions (or at most one isGaussian).
After centering of the variables, the causal relations among x i can be written in the matrix
form: x = Bx + e, where x = (x1 , ..., xn )T , e = (e1 , ..., en )T , and the matrix B can be permuted (bysimultaneous equal row and column permutations) to strict lower triangularity if one knows thecausal order k(i) of xi . We then have e = Wx, where W = I − B. This is exactly the ICA separation
procedure (Hyvärinen et al., 2001). Therefore, the LiNGAM model can be estimated by ICA. We
can permute the rows of the ICA de-mixing matrix W such that it produces a matrix f
W withoutany zero on its diagonal (or in practice, ∑i |f
Wii | is maximized). Dividing each row of f
W by the0fcorresponding diagonal entry gives a new matrix W with all entries on its diagonal equal to 1.
eFinally, by applying equal row and column permutations on B = I − f
W0 , we can find the matrix Bewhich is as close as possible to strictly lower triangularity. B contains the causal relations of xi . For
details, see Shimizu et al. (2006).
6.3 With Nonlinear ICA with Minimal Nonlinear Distortion
We now consider a general case of the nonlinear distortion often encountered in the data generation
procedure, provided that the nonlinear distortion is smooth and mild. We use the MLP structure described in Section 3.1.1, which is a linear transformation coupled with an ordinary MLP, as shownin Figure 13, to model the nonlinear transformation from the the observed variables x i to the disturbance variables ei .
According to Figure 13, we have e = Wx + h(x), and consequently x = (I − W)x − h(x) +e, where h(x) denotes the output of the MLP. As it is difficult to analyze the relations among x i
implied by the nonlinear transformation h(x), we expect that h(x) is weak such that its effect canbe neglected. The linear causal relations among xi can then be discovered by analyzing W.
h(x)

MLP

x

e

W

Figure 13: Structure used to model the transformation from the observed data x i to independent
disturbances ei . h(x) accounts for nonlinear distortion if necessary.
In order to do causality discovery, the separation system in Figure 13 is expected to exhibit thefollowing properties. 1. The outputs ei are mutually independent, since independence of e i is a
crucial assumption in LiNGAM. This can be achieved since nonlinear ICA always has solutions.
2. The matrix W is sparse enough such that it can be permuted to lower triangularity. This can be
enforced by incorporating the L1 (Hyvärinen and Karthikesh, 2000) or smoothly clipped absolutedeviation (SCAD) penalty (Fan and Li, 2001) on the entries of W. 3. The nonlinear mapping modeled by the MLP is weak enough such that we just care about the linear causal relations indicatedby W. To achieve that, we use MISEP with MND given in Section 3.1. In addition, we initialize thesystem with linear ICA results. That is, W is initialized by the linear ICA de-mixing matrix, andthe initial values for weights in the MLP h(x) are very close to 0. The training process is terminatedi (x))
once the LiNGAM property holds for W. After the algorithm terminates, var(hvar(ei ) can be used tomeasure the level of nonlinear distortion in each channel, if needed.
2476

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

6.4 Simulation Study
We examined the performance of the scheme discussed in Section 6.3 for identifying linear causal
relations using simulated data. To make the nonlinear distortion in the data generation procedure
weak, we used the structure in Figure 13 to generate the 8-dimensional observed data x i from someindependent and non-Gaussian variables ei , that is, xi are generated by a linear transformation coupled with a MLP.
The linear transformation in the data generation procedure was generated by A = (I − B) −1 . It
satisfies the LiNGAM property since B was made strict lower triangular. The magnitude of nonzero entries of B is uniformly distributed between 0.05 and 0.5, and the sign is random. To examine
if spurious causal relations would be caused, we also randomly selected 9 entries in the strict lower
triangular part of B and set them to zero. The disturbance variables were obtained by passing
independent Gaussian variables through power non-linearities with the exponent between 1.5 and2. The variances of ei were randomly chosen between 0.2 and 1. These settings are similar to thosein the simulation studies by Shimizu et al. (2006). The sample size is 1000. The nonlinear part isa 8-10-8 MLP with the arctan activation function in the hidden layer. The weights from the inputs
to the hidden layer are between -3 and 3, that is, they are comparatively large, while those fromthe hidden layer to the outputs are small, such that the nonlinear distortion is weak. The nonlinear
distortion level in the generation procedure is measured by the ratio of the variance of the MLP
output to that of the linear output. We considered two cases where the nonlinear distortion level is0.01 and 0.03, respectively.
We used the scheme detailed in Section 6.3 to identify the linear causal relations among x i . The
SCAD penalty was used, and there are 10 arctan hidden units connected to each output of the MLP.
We repeated the simulation for 100 trials. In each trial the maximum iteration number was set to 800.
The results are given in Table 1 (numbers in parentheses are corresponding standard errors). The
failure rate (the chance that LiNGAM does not hold for W within 800 iterations), the percentages ofcorrectly identified non-zero edges, correctly identified large edges (with the magnitude larger thani (x))
0.2), and spurious edges in the successful cases, and the resulting nonlinear distortion level var(hvar(ei )in the separation system are reported. We can see that W almost always satisfies the LiNGAM
property, and that most causal relations (especially large ones) are successfully identified. The
setting λ = 0.12, meaning that MND is explicitly incorporated, gives better results than λ = 0 does,
although the difference is not large. This is not surprising because even with λ = 0, nonlinear ICA
with the separation structure of Figure 13 and with W initialized by linear ICA could achieve MND
to some extent. However, when λ = 0.12, the nonlinear distortion in the separation system is much
weaker, and we found that estimated values of the entries of B are closer to the true ones. The
penalization parameter for SCAD, λSCAD , plays an important role. A larger λSCAD would make W
satisfy the LiNGAM property more easily, but as a price, in the result more causal relations tend todisappear or be weaker.
For comparison, we also used linear ICA with the de-mixing matrix penalized by SCAD 11 forcausality discovery. The result is reported in Table 2. Even when λ SCAD is very large, which causes
many causal relations to disappear, as seen from the table, there is still a high probability thatthe resulting de-mixing matrix fails to satisfy the LiNGAM property. These results show that for
11. The algorithm can be derived by maximizing the ICA likelihood penalized by the SCAD penalty on each entry of thede-mixing matrix. We used the natural gradient learning rule, with the score function adaptively estimated from the
data.

2477

Z HANG AND C HAN

Nonl.
level in F
0.01
0.03

Settings(λ, λSCAD )(0.12, 0.06)(0, 0.06)(0.12, 0.10)(0, 0.10)

Fail. in
800 iter.
3%
3%
1%
1%

Edges identified
88% (11%)87% (12%)79% (14%)76% (12%)Large edges
identified
99% (3%)97% (4%)92% (7%)89% (8%)Spuriousedges
7% (10%)7% (9%)
9% (11%)8% (10%)Nonl. levelvar(hi (x))
var(ei )∼ 0.03∼ 0.06∼ 0.08∼ 0.13

Table 1: Simulation results of identifying linear causal relations among x i with the nonlinear ICA
structure Figure 13 and the SCAD penalty (100 trials). Numbers in parentheses are corresponding standard errors.
Nonl.
level in F
0.01
0.03

Settings

Fail. rate

λSCAD = 0.2
λSCAD = 0.25

41%
54%

Edges identified
67% (13%)52% (12%)Large edges
identified
79% (15%)58% (17%)

Spuriousedges
4% (5%)
4% (7%)

Table 2: Simulation results of identifying linear causal relations among x i by linear ICA with SCAD
penalized de-mixing matrix (100 trials).

the data whose generation procedure has weak nonlinear distortion and approximately satisfies the
LiNGAM property, nonlinear ICA with MND, together with the SCAD penalty, is useful to identify
their linear causal relations.
6.5 Empirical Results
The Hong Kong stock market has some structural features different from the US and UK markets(Ho et al., 2004). One typical feature is the concentration of market activities and equity ownershipin relatively small group of stocks, which probably makes causal relations in the Hong Kong stock
market more obvious.
6.5.1 DATA
Here we aim at discovering the causality network among 14 stocks selected from the Hong Kongstock market.12 The selected 14 stocks are constituents of Hang Seng Index (HSI). 13 They arealmost the largest companies of the Hong Kong stock market. We used the daily dividend/split
adjusted closing prices from Jan. 4, 2000 to Jun. 17, 2005, obtained from the Yahoo financedatabase. For the few days when the stock price is not available, we used simple linear interpolation
to estimate the price. Denoting the closing price of the ith stock on day t by Pit , the corresponding
P −Pi,t−1
return is calculated by xit = itPi,t−1. The observed data are xt = (x1t , ..., x14,t )T . Each return series
contains 1331 samples.
Recently ICA has been exploited as a possible way to explain the driving forces for financial
returns (Back, 1997; Kiviluoto and Oja, 1998; Chan and Cha, 2001). We conjecture that nonlinear
ICA would be more suitable than linear ICA to serve this task, since it seems reasonable that the
12. For saving space, they are not listed here; see the legend in Figure 15.
13. The only exception is Hang Lung Development Co. Ltd (0010.hk), which was removed from HSI on Dec. 2, 2002.

2478

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

ICA mixing model varies slightly for returns at different levels. So we use nonlinear ICA with MND
to analyze the stock returns and to do causality discovery. However, we should be aware that it isprobably very hard to discover causal relations among the selected stocks, since the financial data
are somewhat non-stationary, the data generation mechanism is not clear, and there may be many
confounder variables.
6.5.2 R ESULTS
We first applied a standard ICA algorithm to perform ICA on the data xt . The natural gradientalgorithm (Amari et al., 1996) with the score function adaptively estimated from data was adopted.
We used the LiNGAM software14 to permute W and obtain the matrix B = I− f
W0 . B seems unlikely
to be lower-triangular; in fact, the ratio of the sum of squares of its upper-triangular entries tothat of all entries is 0.24, which is very large. We also exploited linear ICA with the de-mixing
matrix penalized by SCAD to do causality discovery. It was found that the learned de-mixing
matrix W does not follow LiNGAM for λSCAD ≤ 0.25. The value 0.25 for λSCAD is so large thatstatistical independence between outputs is affected. (In fact, most correlations between outputs
have a magnitude larger than 0.1 when λSCAD = 0.25.) We may conclude that the data do not satisfy
the LiNGAM model.
We then adopted the method proposed in Section 6.3. The SCAD penalty was applied to entriesof W with λSCAD = 0.04. The regularization parameter for nonlinear ICA with MND (Eqs. 11
and 16–19) was λ = 0.14. After 195 epochs, W satisfies the LiNGAM assumption and the training
process is terminated. Figure 14 shows the scatter plot of each output e i and its linear part, fromwhich we can see that the nonlinear distortion is weak. Based on the learned W, we found thelinear causal relations among these stocks, as shown in Figure 15. This figure was plotted using the
LiNGAM software.
10

20

10

10

10

0

0

0

0

0

e

−10
−10
10

0

0

−10
−10
10

0

0

0

0

−10
−10

−20
10 −2020−20
10 −2010−10
10 −10

0

0

0

0

0

−10
20 −1010

0−10
10 −1010

0−10
20 −10
10

0

0

0

−10
10 −1010−10
10 −10
10

0

10

0

10

0

0

−10
10 −10

0

10

0

−10
10 −100−10
10 −10Wx

Figure 14: Scatter plot of each output of the system in Figure 13 and its linear part. The nonlinear
i (x))
distortion level var(hvar(ei ) is 0.0485, 0.0145, 0.0287, 0.2075, 0.0180, 0.0753, 0, 0.0001,
0.0193, 0.0652, 0.0146, 0.0419, 0.0544, and 0.0492, respectively, for the 14 outputs e i .
14. It is available at http : //www.cs.helsinki.fi/group/neuroinf/lingam/.

2479

Z HANG AND C HAN

x1: Cheung Kong (0001.hk)
x2: CLP Hldgs (0002.hk)
x3: HK & China Gas (0003.hk)
x4: Wharf (Hldgs) (0004.hk)
x5: HSBC Hldg (0005.hk),
x6: HK Electric (0006.hk)
x7: Hang Lung Dev (0010.hk)
x8: Hang Seng Bank (0011.hk)
x9: Henderson Land (0012.hk)
x10: Hutchison (0013.hk)
x11: Sun Hung Kai Prop (0016.hk)
x12: Swire Pacific ’A’ (0019.hk)
x13: Bank of East Asia (0023.hk)
x14: Cathay Pacific Air (0293.hk)

Figure 15: Causal diagram of the 14 stocks.
Figure 15 gives some interesting findings. 1. Ownership relations tend to cause causal relations.
If A is a holding company of B, there tends to be a causal relation from B to A. There are two
significant relations x8 → x5 and x10 → x1 . In fact, x5 owns some 60% of x8 , and x1 holds about
50% of x10 . 2. Stocks belonging to the same subindex tend to be connected together. For example,
x2 , x3 , and x6 , which are linked together, are the only three constituents of Hang Seng Utilities Index.
x1 , x9 , and x11 are constituents of Hang Seng Property Index. 3. Large bank companies are the causeof many stocks. Here x5 and x8 are the two largest banks in Hong Kong. 4. Returns of stocks in
Hang Seng Property Index tend to depend on many other stocks, while they hardly influence otherstocks. Note that Here x1 , x9 , and x11 are in Hang Seng Property Index.

7. Conclusion
We have proposed the “minimal nonlinear distortion” principle to overcome the ill-posedness of thenonlinear ICA problem. With this principle, the nonlinear ICA solution whose estimated mixing
system is close to linear would be preferred. This principle was implemented by a regularization
technique that minimizes the mean square error of the best linear reconstruction of the observed
mixtures. We explained how the proposed principle overcomes trivial and non-trivial indeterminacies in nonlinear ICA solutions. Experimental results on synthetic data in various situations showedthat nonlinear ICA with minimal nonlinear distortion behaves very well and confirmed our theoretical claims. Since nonlinearity is usually encountered in practice and is not very strong in manycases, nonlinear ICA with minimal nonlinear distortion is expected to be capable of solving somereal-life problems. Its successful application to causality discovery in the Hong Kong stock market illustrated the applicability of the method and the validity of the “minimal nonlinear distortion”
2480

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

principle for some real problems. The result also supports the independent factor model in finance tosome extent. Finally, it should be noted that solutions to nonlinear ICA or nonlinear BSS rely heavily on the prior information on the sources or the mixing mappings. “Minimal nonlinear distortion”
is one type of such information for some problems. If more precise prior information, such as theform of the mixing mapping, the temporal structure of the sources, etc., is available, the separation
result may be more meaningful.

Acknowledgments

This work was partially supported by a grant from the Research Grants Council of the Hong Kong
Special Administration Region, China. We are very grateful to the action editor and the anonymousreferees for their valuable comments and suggestions. The first author would like to thank Haixuan
Yang, Gang Li, and Wan Zhang for helpful discussions.

Appendix A. Gradient of RMSE

(2)

Let H = diag{h0 (a1 ), h0 (a2 ), ..., h0 (aM )}, and W j denote the j-th column of W(2) . We have

∂RMSE (θ)∂W(1)∂RMSE (θ)∂W(d)∂RMSE (θ)∂b(2)∂RMSE (θ)∂b(1)

h M ∂yn n∂a j io
∂yi o
i
K
·
=
E
·
∑ i ∑ ∂a j ∂W(1)
∂W(1)
j=1
i=1
i=1
n M h ∂y T i ∂a on M  ∂a j o
j
(2)T
0
h
(a
)
·
W
·
K
·
= E ∑
K ·
=
E
j
∑
j
∂a j
∂W(1)
∂W(1)
j=1
j=1
= E

n n

∑ Ki ·= E{H · W(2)T · K · xT },
n n
∂yi o
= E ∑ Ki ·
∂W(d)
i=1

(16)

= E{KxT },(17)= E{K},(18)= E{H · W(2)T · K}.

(19)

2481

Z HANG AND C HAN

Appendix B. Gradient of Pi j in Eq. 8
 
 2 
∂yl
∂ yl
∂2
∂
=
Noting that ∂θsince θ is independent from xi , we can obtain the following∂xi ∂x j∂xi ∂x j ∂θrule after tedious derivation:
Pi j
(2)

∂wlm
∂Pi j

(1)

∂wmk
∂Pi j

(1)
∂bm

∂Pi j
∂W(d)
∂Pi j
b(2)

=

∂2 y l∂2 z m·,∂xi ∂x j ∂xi ∂x j

(20)


(1)
(1)
(1)
(1)
= ∆i jm · h00 (am )[wmi · δk j + wm j · δik ] + h000 (am ) · wm j · wmi · xk ,
(1)

(1)

= ∆i jm · h000 (am ) · wm j · wmi ,
= 0,
= 0,

yl
yl
∂ zm
∂ zm
00
where ∆i jm = ∑nl=1 wlm · ∂x∂i ∂x, ∂x∂i ∂x= ∑M
m=1 wlm · ∂xi ∂x j , ∂xi ∂x j = h (am ) · wmi · wm j , and δik is thejj
Kronecker delta function.
(2)

2

2

(2)

2

2

(1)

(1)

Appendix C. Proof of Lemma 1
Proof. The mean square error of reconstructing d from y with the linear transformation a isE{||d − a · y||2 } = E{(d − a · y)T (d − a · y)}= E{dT · d − 2aT d · y + aT a · y2 }noaT d 2 (aT d)2= E aT a · y − T+ dT d− Ta aa an (aT d)2 onT d 2 oa−E+ E{dT d}.
= aT a · E y − T
a a
aT a

(21)

The first term of Eq. 21 is always non-negative. No matter what value a takes, in order to minimize
Eq. 21, we should choose
y = aT d · (aT a)−1(22)
to make this term vanish, meaning that y is the linear combination of d i with the coefficients a ·(aT a)−1 .
Next, when the first term of Eq. 21 vanishes, minimizing this function√ w.r.t. a is reduced to
T
T
2
T
−1
T
T
−1
0
maximizing E{(a d) · (a a) } = E{a dd a · (a a) }. Letting a = a/ aT a, this is equivalent
to the constrained optimization problem: max a0T · E{ddT } · a0 , s.t. a0T a0 = 1. Clearly this is the
PCA problem. So a0 is the eigenvector of E{ddT } associated with the largest eigenvalue, andaccording to Eq. 22, y is the principal component of d multiplied by a constant.
Now let us consider the case where y is constrained to be zero-mean. Let d = E{d}, and d̃ =T
d−d. We have E{||d−a·y||2 } = E{(d̃−a·y+d)T (d̃−a·y+d)} = E{(d̃−a·y)T (d̃−a·y)}+d d.
T
d d can be considered as a constant. Using the result above, we can see that when R y is minimized,
y is the principal component of d̃ multiplied by a constant. (Q.E.D)
2482

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

Appendix D. Proof of Theorem 1
Proof: As it has been assumed here that each output of nonlinear ICA depends only on one of thesources, we can denote by h j (s j ) the estimate of s j produced by nonlinear ICA. For the sake ofsimplicity, we make both xi and h j (s j ) zero-mean, that is, E{xi } = E{h j (s j )} = 0. So the matrixA∗ in Eq. 2 is n × n. Denote by a∗i j the (i, j)th entry of A∗ . RMSE defined by Eq. 2 is
RMSE

=

i

=

n

∑ E xi − ∑ a∗i j h j (s j )
n

j

= ∑E
i

n

h

∑ fi j (s j ) − a∗i j h j (s j )
j

io2



∑ ∑E
i

o2

2
fi j (s j ) − a∗i j h j (s j )

j

+∑E
k6=l




o
.
fik (sk ) − a∗ik hk (sk ) · fil (sl ) − a∗il hl (sl )As E{hk (sk )hl (sl )} = E{hk (sk ) fil (sl )} = 0 for k 6= l, the above equation becomes
RMSE

=

n

i

=

2

∑ ∑ E fi j (s j ) − a∗i j h j (s j ) + ∑ E fik (sk ) fil (sl )nj∑ ∑ E fi j (s j ) − a∗i j h j (s j )
j

i

2 o

k6=l

o

+ const.

One can see that minimization of the above function can be achieved by minimizing ∑i E fi j (s j ) −2a∗i j h j (s j ) independently for each j. That is, h j (s j ) and a∗i j are adjusted to minimize ∑i E fi j (s j ) −2a∗i j h j (s j ) . According to Lemma 1, h j (s j ) produced by nonlinear ICA with MND is the firstprincipal component of f∗ j (s j ) = [ f1 j (s j ), · · · , fn j (s j )]T , multiplied by a constant. (Q.E.D)

Appendix E. Proof of Theorem 2
Proof. Denote by h j (s j ) the estimate of s j produced by nonlinear ICA, and assume that both x i andh j (s j ) zero-mean. Denote by a∗i j the (i, j)th entry of A∗ . Note that ∑ j,k,l Oi, jkl · s j sk sl = ∑ j Oi, j j j ·s j s j s j + 3 · ∑ j ∑k6= j Oi, jkk · s j s2k + ∑ j ∑k6= j ∑ l6= j Oi, jkl · s j sk sl , and that E{s j } = 0 and E{s2j } = 1. RMSE
l6=kdefined by Eq. 2 becomesno2
RMSE = ∑ E xi − ∑ a∗i j h j (s j )
=

n

i

j

1

1

∑ E fi (0) + ∑ Oi, j · s j + 2 ∑ Oi, jk · s j sk + 6 ∑ Oi, jkl · s j sk sl − ∑ a∗i j h j (s j )
i
n

j

j,k

j,k,l

j

o2

n

h
1
3
1
= ∑ E fi (0) + ∑ Oi, j · s j + Oi, j j · s2j + Oi, j j j · s3j + ∑ Oi, jkk · s j s2k
2
6
6
j
i=1
k6= j
i 1
o2
1
−a∗i j h j (s j ) + ∑ ∑ Oi, jk · s j sk + ∑ ∑ ∑ Oi, jkl · s j sk sl .
2 j k6= j
6 j k6= j l6= j
l6=k

2483

(23)

Z HANG AND C HAN

Bearing in mind that s j are mutually independent, and also taking all the terms independent ofh j (s j ) and a∗i j as constants, we can re-write Eq. 23 as
RMSE
io2
n h
1
3
1
+ const= ∑ E ∑ Oi, j · s j + Oi, j j · s2j + Oi, j j j · s3j + ∑ Oi, jkk · s j s2k − a∗i j h j (s j )2
6
6 k6= j
i
j
n h
i 1
o2
1
1
= ∑ E ∑ Oi, j · s j + Oi, j j · s2j + Oi, j j j · s3j − a∗i j h j (s j ) + ∑ ∑ Oi, jkk · s j s2k + const
2
6
2 j k6= j
j
i
nh 
i
2
1
1
= ∑ E ∑ Oi, j · s j + Oi, j j · s2j + Oi, j j j · s3j − a∗i j h j (s j )
2
6
i
j
o

− ∑ a∗i j h j (s j ) · ∑ Oi, jkk · s j s2k + const
j

=

n



k6= j

2
1
1
2
3
∗
O
·
s
+
O
·
s
+
O
·
s
−
a
h
(s
)
∑ i, j j 2 i, j j j 6 i, j j j j i j j jijo− ∑ a∗i j h j (s j ) · ∑ Oi, jkk · s j + const

∑E

k6= j

j

n

o2

1
1
1
Oi, jkk · s j + Oi, j j · s2j + Oi, j j j · s3j − a∗i j h j (s j ) + const∑
2 k6= j
2
6
i j
h
i

2
= ∑ ∑ E Di j (s j ) − a∗i j h j (s j ) + const.

=

∑∑E
j

Oi, j +

i

Note that there is no dependence relationship between h j (·), as well as a∗i j , with different j. To
minimize the above function, we just need to adjust h j (s j ) and a∗i j to minimize ∑i E Di (s j ) −2a∗i j h j (s j ) , independently for each j. According to Lemma 1, h j (s j ) is the first principal component of D̃∗ j (s j ) = [D̃1 j (s j ), · · · , D̃n j (s j ), ]T , multiplied by a constant. (Q.E.D)

References
L.B. Almeida. MISEP - linear and nonlinear ICA based on mutual information. Journal of Machine
Learning Research, 4:1297–1318, 2003.
L.B. Almeida. Separating a real-life nonlinear image mixture. Journal of Machine Learning Research, 6:1199–1229, 2005.
S. Amari, A. Cichocki, and H.H. Yang. A new learning algorithm for blind signal separation. In
David S. Touretzky, Michael C. Mozer, and Michael E. Hasselmo, editors, Advances in Neural
Information Processing Systems, volume 8, pages 757–763, Cambridge, MA, 1996. MIT Press.
L.A. Baccala and K. Sameshima. Partial directed coherence: a new concept in neural structure
determination. Biol. Cybern., 84:463–474, 2001.
A.D. Back. A first application of independent component analysis to extracting structure from stock
returns. International Journal of Neural Systems, 8(4):473–484, August 1997.
A.J. Bell and T.J. Sejnowski. An information-maximization approach to blind separation and blind
deconvolution. Neural Computation, 7(6):1129–1159, 1995.
2484

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

C.M. Bishop. Curvature-driven smoothing: a learning algorithm for feedforward networks. IEEE
Transactions on Neural Networks, 4(5):882–884, 1993.
C.M. Bishop. Regularization and complexity control in feed-forward networks. In Proc. International Conference on Artificial Neural Networks (ICANN’95), volume 1, pages 141–148, 1995.
G. Burel. Blind separation of sources: a nonlinear neural algorithm. Neural Networks, 5(6):937–
947, 1992.
J.F. Cardoso. Blind signal separation: Statistical principles. Proceeding Of The IEEE, special issue
on blind identification and estimation, 9(10):2009–2025, 1998.
L. Chan and S.M. Cha. Selection of independent factor model in finance. In proceedings of 3rdInternational Conference on Independent Component Analysis and blind Signal Separation, San
Diego, California, USA, December 2001.
S.S. Chen and R.A. Gopinath. Gaussianization. In Todd K. Leen, Thomas G. Dietterich, and Volker
Tresp, editors, Advances in Neural Information Processing Systems 13, pages 423–429. MIT
Press, 2001.
A. Cichocki and S. Amari. Adaptive Blind Signal and Image Processing: Learning Algorithms andApplications. John Wiley & Sons, UK, corrected and revisited edition, 2003.
Y. Dodge and V. Rousson. On asymmetric properties of the correlation coefficient in the regression
setting. The American Statistician, 55(1):51–54, 2001.
J. Eriksson and V. Koivunen. Blind identification of class of nonlinear instantaneous ICA models. In
Proc. of the XI European SIgnal Proc. Conf. (EUSIPCO 2002), volume 2, pages 7–10, Toulouse,
France, Sept. 2002.
J. Fan and R. Li. Variable selection via nonconcave penalized likelihood and its oracle properties.
J. Amer. Statist. Assoc., 96:1348–1360, 2001.
C. Granger. Testing for causality: A personal viewpoint. Journal of Economic Dynamics andControl, 2:329-352, 1980.
W.E.L. Grimson. A computational theory of visual surface interpolation. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 298:395–427, 1982.
S. Harmeling, A. Ziehe, M. Kawanabe, and K.R. Müller. Kernel-based nonlinear blind source
separation. Neural Computation, 15:1089–1124, 2003.
R.Y. Ho, R. Strange, and J. Piesse. The structural and institutional features of the Hong Kong stock
market: Implications for asset pricing. Research Paper 027, The Management Centre Research
Papers, King’s College London, 2004.
A. Hyvärinen. Fast and robust fixed-point algorithms for independent component analysis. IEEE
Transactions on Neural Networks, 10(3):626–634, 1999.
2485

Z HANG AND C HAN

A. Hyvärinen and R. Karthikesh. Sparse priors on the mixing matrix in independent component
analysis. In Proc. Int. Workshop on Independent Component Analysis and Blind Signal Separation (ICA2000), pages 477–452, Helsinki, Finland, 2000.
A. Hyvärinen and P. Pajunen. Nonlinear independent component analysis: Existence and uniqueness results. Neural Networks, 12(3):429–439, 1999.
A. Hyvärinen, J. Karhunen, and E. Oja. Independent Component Analysis. John Wiley & Sons, Inc,
2001.
C. Jutten and A. Taleb. Source separation: From dusk till dawn. In 2nd International Workshop
on Independent Component Analysis and Blind Signal Separation (ICA 2000), pages 15–26,
Helsinki, Finland, 2000.
A.M. Kagan, Y.V. Linnik, and C.R. Rao. Characterization Problems in Mathematical Statistics.
Wiley, New York, 1973.
T. Khanna and J.W. Rivkin. Interorganizational ties and business group boundaries: Evidence froman emerging economy. Organization Science, 17(3):333-352, 2006.
K. Kiviluoto and E. Oja. Independent component analysis for parallel financial time series. In Proc.
ICONIP’98, volume 2, pages 895–898, Tokyo, Japan, 1998.
H. Lappalainen and A. Honkela. Bayesian nonlinear independent component analysis by multilayer perceptron. In M.Girolami, editor, Advances in Independent Component Analysis, pages
93–121. Spring-Verlag, 2000.
J. Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, Cambridge,
2000.
T. Poggio, V. Torre, and C. Koch. Computational vision and regularization theory. Nature, 317:
314–319, 1985.
M. Reale and G. Tunnicliffe Wilson. Identification of vector ar models with recursive structural
errors using conditional independence graphs. Statistical Methods and Applications, 10(1-3):
49–65, 2001.
S. Shimizu, P.O. Hoyer, A. Hyvärinen, and A.J. Kerminen. A linear non-Gaussian acyclic model
for causal discovery. Journal of Machine Learning Research, 7:2003–2030, 2006.
A. Taleb. A generic framework for blind source separation in structured nonlinear models. IEEE
Transactions on Signal Processing, 50(8):1819–1830, 2002.
A. Taleb and C. Jutten. Source separation in post-nonlinear mixtures. IEEE Trans. on Signal
Processing, 47(10):2807–2820, 1999.
Y. Tan, J. Wang, and J. M. Zurada. Nonlinear blind source separation using a radial basis function
network. IEEE Trans. on Neural Networks, 12(1):124–134, 2001.
A.N. Tikhonov and V.A. Arsenin. Solutions of Ill-posed Problems. Winston & Sons, Washington,
1977.
2486

M INIMAL N ONLINEAR D ISTORTION FOR N ONLINEAR ICA

H. Valpola. Nonlinear independent component analysis using ensemble learning: Theory. In Proc.
Int. Workshop on Independent Component Analysis and Blind Signal Separation (ICA2000),
pages 251–256, Helsinki, Finland, 2000.
L. Xu. Least mean square error reconstruction principle for self-organizing neural-nets. Neural
Networks, 6:627–648, 1993.
H.H. Yang, S. Amari, and A. Cichocki. Information-theoretic approach to blind separation ofsources in nonlinear mixture. Signal Processing, 64(3):291–300, 1998.
K. Zhang and L. Chan. Kernel-based nonlinear independent component analysis. In Proc. Int.
Workshop on Independent Component Analysis and Signal Separation (ICA2007), pages 301–
308, London, UK, Sept. 2007a.
K. Zhang and L. Chan. Nonlinear independent component analysis with minimum nonlinear distortion. In the 24th Annual International Conference on Machine Learning (ICML 2007), pages
1127–1134, Corvallis, OR, US, Jun. 2007b.
K. Zhang and L. Chan. Extensions of ICA for causality discovery in the Hong Kong stock market.
In Proc. 13th International Conference on Neural Information Processing (ICONIP 2006), pages
400–409, Hong Kong, 2006.

2487