ARTICLE

Communicated by Shun-ichi Amari

On the Achievability of Blind Source Separationfor High-Dimensional Nonlinear Source Mixtures
Takuya Isomura

Taro Toyoizumi
taro.toyoizumi@riken.jp
Laboratory for Neural Computation and Adaptation, RIKEN Center for Brain
Science, Wako, Saitama 351-0198, Japan, and Department of Mathematical
Informatics, Graduate School of Information Science and Technology,
University of Tokyo, Bunkyo-ku, Tokyo 113-8656, JapanFor many years, a combination of principal component analysis (PCA)
and independent component analysis (ICA) has been used for blindsource separation (BSS). However, it remains unclear why these linear methods work well with real-world data that involve nonlinear
source mixtures. This work theoretically validates that a cascade of linear PCA and ICA can solve a nonlinear BSS problem accurately—when
the sensory inputs are generated from hidden sources via nonlinear mappings with sufficient dimensionality. Our proposed theorem, termed theasymptotic linearization theorem, theoretically guarantees that applying
linear PCA to the inputs can reliably extract a subspace spanned by thelinear projections from every hidden source as the major components—
and thus projecting the inputs onto their major eigenspace can effectively
recover a linear transformation of the hidden sources. Then subsequent
application of linear ICA can separate all the true independent hidden
sources accurately. Zero-element-wise-error nonlinear BSS is asymptotically attained when the source dimensionality is large and the input dimensionality is sufficiently larger than the source dimensionality. Our
proposed theorem is validated analytically and numerically. Moreover,the same computation can be performed by using Hebbian-like plasticity rules, implying the biological plausibility of this nonlinear BSS strategy. Our results highlight the utility of linear PCA and ICA for accurately
and reliably recovering nonlinearly mixed sources and suggest the importance of employing sensors with sufficient dimensionality to identify
true hidden sources of real-world data.

Neural Computation 33, 1433–1468 (2021)
https://doi.org/10.1162/neco_a_01378© 2021 Massachusetts Institute of Technology.
Published under a Creative Commons
Attribution 4.0 International (CC BY 4.0) license.

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

takuya.isomura@riken.jp
Laboratory for Neural Computation and Adaptation and Brain Intelligence Theory
Unit, RIKEN Center for Brain Science, Wako, Saitama 351-0198, Japan

1434

T. Isomura and T. Toyoizumi

1 Introduction

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Blind source separation (BSS) involves the separation of mixed sensory inputs into their hidden sources without knowledge of the manner in which
they were mixed (Cichocki, Zdunek, Phan, & Amari, 2009; Comon & Jutten, 2010). Among the numerous BSS methods, a combination of principal
component analysis (PCA) (Pearson, 1901; Oja, 1982, 1989; Sanger, 1989; Xu,
1993; Jolliffe, 2002) and independent component analysis (ICA) (Comon,
1994; Bell & Sejnowski, 1995, 1997; Amari, Cichocki, & Yang, 1996; Hyvarinen & Oja, 1997) is one of the most widely used approaches. In thiscombined PCA–ICA approach, PCA yields a low-dimensional concise representation (i.e., the major principal components) of sensory inputs thatmost suitably describes the original high-dimensional redundant inputs.
Whereas, ICA provides a representation (i.e., encoders) that separates thecompressed sensory inputs into independent hidden sources. A classicalsetup for BSS assumes a linear generative process (Bell & Sejnowski, 1995),
in which sensory inputs are generated as a linear superposition of independent hidden sources. The linear BSS problem has been extensively studied
both analytically and numerically (Amari, Chen, & Cichocki, 1997; Oja &Yuan, 2006; Erdogan, 2009), where the cascade of PCA and ICA is guaranteed to provide the optimal linear encoder that can separate sensory inputs
into their true hidden sources, up to permutations and sign-flips (Baldi &Hornik, 1989; Chen, Hua, & Yan, 1998; Papadias, 2000; Erdogan, 2007).
Another crucial perspective is the applicability of BSS methods to realworld data generated from a nonlinear generative process. In particular,
the aim of nonlinear BSS is to identify the inverse of a nonlinear generative
process that generates sensory inputs and thereby infer their true independent hidden sources based exclusively on the sensory inputs. Although thecascade of linear PCA and ICA has been applied empirically to real-world
BSS problems (Calhoun, Liu, & Adali, 2009), no one has yet theoretically
proven that this linear BSS approach can solve a nonlinear BSS problem. To
address this gap, this work demonstrates mathematically that the cascadeof PCA and ICA can solve a nonlinear BSS problem accurately when thesource dimensionality is large and the input dimensionality is sufficiently
larger than the source dimensionality so that various nonlinear mappingsfrom sources to inputs can be determined from the inputs themselves.
In general, there are five requirements for solving the nonlinear BSS
problem. The first two requirements are related to the representation capacity of the encoder: (1) the encoder’s parameter space must be sufficiently
large to accommodate the actual solution that can express the inverse of
the true generative process. (2) However, this parameter space should not
be too large; otherwise, a nonlinear BSS problem can have infinitely many
spurious solutions wherein all encoders are independent but dissimilar to
the true hidden sources (Hyvarinen & Pajunen, 1999; Jutten & Karhunen,2004). Hence, it is important to constrain the representation capacity of theAchievability of Nonlinear Blind Source Separation

1435

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

encoder in order to satisfy these opposing requirements. A typical approachfor solving the nonlinear BSS problem involves using a multilayer neural network—with nonlinear activation functions—that learns the inverseof the generative process (Lappalainen & Honkela, 2000; Karhunen, 2001;
Hinton & Salakhutdinov, 2006; Kingma & Welling, 2013; Dinh, Krueger, &Bengio, 2014). The remaining three requirements are related to the unsupervised learning algorithms used to identify the optimal parameters forthe encoder: (3) the learning algorithm must have a fixed point at which
the network expresses the inverse of the generative process; (4) the fixed
point must be linearly stable so that the learning process converges to thesolution; and (5) the probability of not converging to this solution should
be small, that is, most realistic initial conditions must be within the basin of
attraction of the true solution.
Approaches using a nonlinear multilayer neural network satisfy requirement 1 when the number of neurons in each layer is sufficient (universality: Cybenko, 1989; Hornik, Stinchcombe, & White, 1989; Barron, 1993);
moreover, learning algorithms that satisfy requirements 3 and 4 are also
known (Dayan, Hinton, Neal, & Zemel, 1995; Friston, 2008; Friston, TrujilloBarreto, & Daunizeau, 2008). However, reliable identification of the true
hidden sources is still necessary because the encoder can have infinitely
many spurious solutions if its representation capacity is too large (i.e., if requirement 2 is violated). As previously indicated, to the best of our knowledge, there is no theoretical proof that confirms a solution for a nonlinear
BSS problem (Hyvarinen & Pajunen, 1999; Jutten & Karhunen, 2004), except for some cases wherein temporal information—such that each independent source has its own dynamics—is available (Hyvarinen & Morioka,
2016, 2017; Khemakhem, Kingma, Monti, & Hyvarinen, 2020). Moreover,even when requirement 2 is satisfied, there is no guarantee that a learningalgorithm will converge to the true hidden source representation because
it might be trapped in a local minimum wherein outputs are still not independent of each other. Thus, for a nonlinear BSS problem, it is of paramountimportance to simplify the parameter space of the inverse model in orderto remove spurious solutions and prevent the learning algorithm from attaining a local minimum (to satisfy requirements 2 and 5) while retaining
its capacity to represent the actual solution (requirement 1). Hence, in thiswork, we apply a linear approach to solving a nonlinear BSS problem inorder to ensure requirements 2 and 5 are satisfied. We demonstrate that thecascade of PCA and ICA can reliably identify a good approximation of theinverse of nonlinear generative processes asymptotically under the condition where the source dimensionality is large and the input dimensionality
is sufficiently larger than the source dimensionality (thus satisfying requirements 1–5). Although such a condition is different from the case typicallyconsidered by earlier work—where the sources and inputs have the same
dimensionality—the condition we consider turns out to be apt for the mathematical justification of the achievability of the nonlinear BSS.

1436

T. Isomura and T. Toyoizumi

2 Results
2.1 Overview. Our proposed theorem, referred to as asymptotic linearization theorem, is based on an intuition that when the dimensionalityof sensory inputs is significantly larger than that of hidden sources, these
inputs must involve various linear and nonlinear mappings of all hidden
sources—thus providing sufficient information to identify the true hidden
sources without ambiguity using an unsupervised learning approach. We
consider that Ns -dimensional hidden sources s ≡ (s1 , . . . , sNs )T generate Nx dimensional sensory inputs x ≡ (x1 , . . . , xNx )T through a generative processcharacterized by an arbitrary nonlinear mapping x ≡ F(s) (see Figure 1).
Here, the hidden sources are supposed to follow independently an identicalprobability distribution with zero mean and unit variance ps (s) ≡ i pi (si ).
A class of nonlinear mappings F(s) can be universally approximated by aspecific but generic form of two-layer network. Suppose A ∈ RN f ×Ns andB ∈ RNx ×N f as higher- and lower-layer mixing matrices, respectively; a ∈RN f as a constant vector of offsets; y ≡ As as a linear mixture of sources,f (•) : R → R as a nonlinear function; and f ≡ ( f1 , . . . , fN f )T ≡ f (As +a) ≡ f (y + a) as N f -dimensional nonlinear bases. The sensory inputs are
given asx = B f (As + a),(2.1)

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Figure 1: Structures of nonlinear generative process (top) and linear neuralnetwork (bottom). Hidden sources s = (s1 , . . . , sNs )T generate sensory inputs
x = (x1 , . . . , xNx )T via nonlinear mappings characterized by nonlinear bases f =( f1 , . . . , fN f )T , where A = {Akl } and a = (a1 , . . . , aN f )T are gaussian-distributed
higher-layer mixing weights and offsets, respectively, and B = {B jk } is a lowerlayer mixing weight matrix. Encoders comprise a single-layer linear neural network, where u = (u1 , . . . , uNs )T are neural outputs and W = {Wi j } is a synaptic
weight matrix. Equations on the right-hand-side panels summarize the generation of sensory inputs from hidden sources through nonlinear mixtures and theinversion of this process using a linear neural network.

Achievability of Nonlinear Blind Source Separation

1437

x − E[x] = BHs
Bφ .
 + signal(2.2)

residual

The first term on the right-hand side represents the signal comprising thelinear components of the hidden sources, whereas the second term represents the residual introduced by the nonlinearity in the generative process.
Further, because φ is uncorrelated with s, the covariance of the bases isdecomposed into their linear and nonlinear counterparts Cov[ f ] ≡ E[( f −E[ f ])( f − E[ f ])T ] = Cov[Hs] + Cov[φ] = HH T + , where  ≡ Cov[φ] indicates the covariance of φ. Thus, the covariance of the sensory inputs

1

In this work, based on the literature (Cybenko, 1989), we define the universality in thecomponent-wise sense as the condition wherein each element of x = F(s) is approximatedusing the two-layer network with an approximation error smaller than an arbitrary δ > 0,that is, sups [|Fj (s) − B j f (As + a)|] < δ for j = 1, . . . , Nx when N f is sufficiently large. Notethat Fj (s) is the jth element of F(s) and B j is the jth row vector of B. It is known that when
f is a sigmoidal function and parameters are selected appropriately, the approximation(Barron, 1993); hence, E[|Fj (s) − B j f (As +error can be upper-bounded by the order N−1fa)|2 ] ≤ O(N−1). This relationship holds true when N f  Ns  1, irrespective of Nx .
f

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

or, equivalently, x = B f (y + a) = B f . This expression using the two-layer
network is universal in the component-wise sense1 (Cybenko, 1989; Hornik
et al., 1989; Barron, 1993) and each element of x can represent an arbitrary
mapping x = F(s) as N f increases by adjusting parameters a, A, and B.
We further suppose that each element of A and a is independently generated from a gaussian distribution N [0, 1/Ns ], which retains its universality (Rahimi & Recht, 2008a, 2008b) as long as B is tuned to minimize themean squared error E[|B f − F(s)|2 ]. Here, E[•] describes the average overps (s). The scaling of A is to ensure that the argument of f is of order 1.
The Ns−1/2 order offset a is introduced to this model to express any generative process F(s); however, it is negligibly small relative to As for large Ns .
The whole system, including the generative process and neural network, is
depicted in Figure 1 (left). The corresponding equations are summarized in
Figure 1 (right).
For analytical tractability, we decompose nonlinear bases f into thesum of linear and nonlinear parts of hidden sources as follows: linear
components of the hidden sources in the bases are defined as Hs using a
coefficient matrix H that minimizes the mean squared error, H ≡arg minH E[| f − E[ f ] − Hs|2 ]. Such a coefficient matrix is computed as H =E[ f sT ]. The remaining part φ ≡ f − E[ f ] − Hs is referred to as nonlinear
components of the hidden sources, which are orthogonal (uncorrelated) tos (i.e., E[φsT ] = O). This definition of linear and nonlinear components is
unique. Thus, the sensory inputs (see equation 2.1) are decomposed intolinear and nonlinear transforms of the hidden sources:

1438

T. Isomura and T. Toyoizumi

Cov[x] ≡ BCov[ f ]BT can be decomposed into the signal and residual
covariances:
T T
Cov[x] = BHH
 B  +signal covariance

T
BB
 

.

(2.3)

residual covariance2

In this article, we refer to the neural network as the encoder because networks thatconvert the input data into a different (lower-dimensional) code or representation are
widely recognized as encoders in the literature on machine learning (Hinton & Salakhutdinov, 2006; Goodfellow, Bengio & Courville, 2016). However, one may think that thegenerative process encodes the hidden sources into the sensory inputs through nonlinear mixtures. From this viewpoint, one may term the neural network as the decoder that
decodes the inputs to recover the original sources.
3
Our numerical simulations suggest that the system behaves similarly for Nx ≥ N f ≥Ns > 1 in some cases, although our theorem holds mathematically only when Nx = N f 
Ns  1.

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Crucially, the signal covariance has only Ns nonzero eigenvalues when
N f > Ns owing to low-column-rank matrix H, whereas the eigenvalues of
the residual covariance are distributed in the N f -dimensional eigenspace.
This implies that if the norms of the linear and nonlinear components in the
inputs are in a similar order and Ns is sufficiently large, the eigenvalues ofthe signal covariance (i.e., linear components) are N f /Ns order times larger
than those of the residual covariance (i.e., nonlinear components). We will
prove in the following sections that when Ns  1, this property is derived
from the fact that elements of A are gaussian distributed and singular valuesof B are of order 1.
In what follows, we demonstrate that owing to the aforementioned property, the first Ns major principal components of the input covariance precisely match the signal covariance when the source dimensionality is large
and the input dimensionality is sufficiently larger than the source dimensionality. Consequently, projecting the inputs onto the subspace spannedby major eigenvectors can effectively extract the linear components in the
inputs. Moreover, the same projection can effectively filter out the nonlinear components in the inputs because the majority of the nonlinear components are perpendicular to the major eigenspace. Thus, applying PCA to
the inputs enables the recovery of a linear transformation of all the true hidden sources of the mixed sensory inputs with a small estimation error. This
property, termed asymptotic linearization, enables the reduction of the original nonlinear BSS problem to a simple linear BSS problem, consequently
satisfying requirements 1 to 5. In the remainder of this article, we mathematically validate this theorem for a wide range of nonlinear setups and thereby
demonstrate that linear encoder neural networks can perform the nonlinear BSS in a self-organizing manner.2 We assume that Nx = N f  Ns  1throughout the article.3

Achievability of Nonlinear Blind Source Separation

1439



H = diag[EN [ f ]]A + O Ns−3/2 = f A + O Ns−1 .

(2.4)

In the last equality,√ the scalar coefficient is expressed as f ≡f (ξ ) exp(−ξ 2 /2)dξ / 2π using the expectation over a unit gaussian variable ξ . The order Ns−1 error term includes the nongaussian contributions

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

2.2 PCA Can Extract Linear Projections of All True Hidden Sources.
In this section, we first demonstrate that the major principal componentsof the input covariance precisely match the signal covariance when N f 
Ns  1 by analytically calculating eigenvalues of HH T and  for the aforementioned system. For analytical tractability, we assume that f (•) is an odd
nonlinear function. This assumption does not weaken our proposed theorem because the presumed generative process in equation 2.1 remains universal. We further assume that ps (s) is a symmetric distribution.
Each element of—and any pair of two elements in—a vector y ≡ As isapproximately gaussian distributed for large Ns due to the central limit
theorem. The deviation of their marginal distribution p(yi , y j ) from thecorresponding zero-mean gaussian distribution pN (yi , y j ) ≡ N [0, ÃÃT ] isorder Ns−1 because the source distribution is symmetric, where Ã ∈ R2×Ns
indicates a submatrix of A comprising its ith and jth rows (see lemma
1 and its proof in section 4.1). This asymptotic property allows us tocompute H and  based on the expectation over a tractable gaussian
distribution pN (yi , y j )—as a proxy for the expectation over p(yi , y j )—asthe leading order for large Ns , despite the fact that s actually follows anongaussian distribution (ensure that p(yi , y j ) converges to pN (yi , y j )in the large
Ns limit). We denote the expectation over pN (yi , y j ) as
EN [•] ≡ •pN (yi , y j )dyi dy j to distinguish it from E[•]. The latter canbe rewritten as E[•] = EN [•(1 + G(yi , y j )/Ns )] using an order-one function
G(yi , y j ) that characterizes the deviation of p(yi , y j ) from pN (yi , y j ); thus,
the deviation caused by this approximation is negligibly smaller than theleading order in the following analyses.
Owing to this asymptotic property, the coefficient matrix H is computed as follows: the pseudo-inverse of A, A+ ≡ (AT A)−1 AT , satisfies
A+ A = I and A+T = AT+ . Thus, we have H = HAT A+T = E[ f yT ]A+T .
It can be approximated by EN [ f yT ]A+T as the leading order fromlemma
1. From the  integration by parts, we obtain EN [ f yT ] =− f pN (y)T dy AAT = diag[ f ]pN (y)dy AAT = diag[EN [ f ]]AAT . This
relationship is also known as the Bussgang theorem (Bussgang, 1952).
Thus, EN [ f yT ]A+T = diag[EN [ f ]]A is order Ns−1/2 for a generic odd function f (•) (as it is the product of order 1 diagonal matrix diag[EN [ f ]] andorder Ns−1/2 matrix A). According to lemma 1, the difference between H
and diag[EN [ f ]]A is order Ns−3/2 as it is Ns−1 order times smaller than theleading order (see remark 1 in section 4.1 for details). Hence, we obtain

1440

T. Isomura and T. Toyoizumi∞

CovN [φi , φ j ] =n=1

EN φi(n) EN φ (n)jn!

CovN [yi , y j ]n

(2.5)

for i = j. Here, CovN [yi , y j ]n is order Ns−n/2 . Because EN [φi(1) ] = O(Ns−1 ),EN [φi(2) ] = O(Ns−1/2 ) and |EN [φi(n) ]| ≤ O(1) for n ≥ 3 hold with a generic
odd function f (•) (see remark 2 in section 4.2 for details), CovN [φi , φ j ] isorder Ns−3/2 . In contrast, the ith diagonal element CovN [φi , φi ] is order 1.
These observations conclude that eigenvalues of CovN [φ]—and therefore
those of  according to lemma 1 (see remark 1 in section 4.1 for details)—
are not larger than order max[1, N f Ns−3/2 ]; thus, all nonzero eigenvaluesof HH T —which are order N f /Ns —are sufficiently greater than those of 
when N f  Ns  1.
One can further proceed to the calculation of  by explicitly computing the coefficients EN [φi(n) ] up to the fourth order. Because f (•) is an
odd function, using CovN [yi , y j ]n = ((AAT ) n )i j , equation 2.5 becomes
CovN [φi , φ j ] = EN [ f (3) (yi )]EN [ f (3) (y j )]{ai a j ((AAT ) 2 )i j /2 + ((AAT ) 3 )i j /6}+ O(Ns−5/2 ) (see remark 2 in section 4.2 for details). This analytical expression involves the Hadamard (element-wise) power of matrix AAT(denoted by ), for example, (AAT ) 3 ≡ (AAT ) (AAT ) (AAT ). WhenN f > Ns2  1, (AAT ) 3 has Ns major eigenvalues—all of which are 3N f /Ns2
(as the leading order) with the corresponding eigenvectors that match the

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

of y and the effect of the order Ns−1/2 deviation of diag[EN [ f ]] from f ,where the latter yields the order Ns−1 error owing to the product of orderNs−1/2 diagonal matrix (diag[EN [ f ]] − f I) and order Ns−1/2 matrix A. Dueto the low column rank of matrix A, AAT has Ns nonzero eigenvalues,
all of which are N f /Ns as the leading order, because elements of A are
independent and identically distributed variables sampled from a gaussian distribution N [0, 1/Ns ] (Marchenko & Pastur, 1967). The same scalingof eigenvalues also holds for HH T because the difference between thenonzero singular values of H and those of f A—caused by the order Ns−1/2
deviation of diag[EN [ f ]] from f I—is smaller than order (N f /Ns )1/2 .
Next, we characterize the covariance matrix  = Cov[φ] using theasymptotic property. The nonlinear components can be cast as a function ofy, φ(y) = f (y + a) − E[ f ] − E[ f yT ]A+T A+ y. For large Ns , the covariance of y,
Cov[y] = AAT , is close to the identity matrix—and the deviation (AAT − I)is order Ns−1/2 for each element—because elements of A are gaussian distributed. This weak correlation between yi and y j (i = j) leads to a weak correlation of their functions φi and φ j . Thus, by computing the Taylor expansion of CovN [φi , φ j ] ≡ EN [φi φ j ] − EN [φi ]EN [φ j ] with respect to the small
covariance CovN [yi , y j ] ≡ EN [yi y j ], we obtain (Toyoizumi & Abbott, 2011)(see lemma 2 in section 4.2):

Achievability of Nonlinear Blind Source Separation

1441

=

f2 − f

2



2

f (3)
I+
(AAT + aaT ) + .
2Ns

(2.6)

Here, from  = Cov[ f ] − HH T and equation 2.4, each diagonal element is2expressed as ii = f 2 − f + O(Ns−1/2 ) using f 2 ≡ f 2 (ξ ) exp(−ξ 2 /2)dξ /√2π, which yields the first term of equation 2.6 (where the error term is
involved in the third term). The second term is generated from equation
2.5 followed by lemma 3, where EN [ f (3) (yi )] is approximated by f (3) ≡√ (3)f (ξ ) exp(−ξ 2 /2)dξ / 2π up to order Ns−1/2 . The third term is the error matrix  that summarizes the deviations caused by the nongaussianity of p(yi , y j ) (see lemma 1; see also remark 1 in section 4.1), higher-orderterms of the Taylor series in equation 2.5, minor eigenmodes of (AAT ) 2 and(AAT ) 3 , and the effect of the order Ns−1/2 deviations of coefficients (e.g.,E[ fi2 ]) from their approximations (e.g., f 2 ). Due to its construction, eigenvalues of  are smaller than those of the first or second term of equation
2.6. This indicates that for large Ns , either the first or second term of equation 2.6 provides the largest eigenvalue of , which is order max[1, N f /Ns2 ].
Thus,  is negligible for the following analyses.
Accordingly, all nonzero eigenvalues of HH T —that are order N f /Ns —are much greater than the maximum eigenvalue of —that is, ordermax[1, N f /Ns2 ]—when N f  Ns  1. Thus, unless B specifically attenuates
one of Ns major eigenmodes of HH T or significantly amplifies a particulareigenmode of , we havemin eig[H T BT BH]  max eig[BBT ](2.7)

for Nx = N f  Ns  1. Here, eig[•] indicates a set of eigenvalues of •.
Because min eig[H T BT BH] is not smaller than order N f /Ns min eig[BT B],
while max eig[BBT ] is not larger than order max[1, N f /Ns2 ] max eig[BT B],
inequality 2.7 holds at least when N f /Ns and Ns are much greater thanmax eig[BT B]/ min eig[BT B]. This is the case, for example, if all the singular values of B are order 1. We call B sufficiently isotropic when

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

directions of AAT —and (N f − Ns ) minor eigenvalues that are negligibly
smaller than order N f /Ns2 (see lemma 3 in section 4.3). This property yields
an approximation (AAT ) 3 − diag[(AAT ) 3 ] = (3/Ns )(AAT − diag[AAT ])up to the negligible minor eigenmodes. Note that when Ns2 > N f  1,
off-diagonal elements of (AAT ) 3 are negligible compared to diagonal
elements of  (see below). Similarly, (AAT ) 2 has one major eigenvalue,
N f /Ns , in the direction of (1, ..., 1)T , while other minor eigenmodes of it
are negligible. These approximations are used to compute off-diagonal
elements of .
Hence,  is analytically expressed as

1442

T. Isomura and T. Toyoizumi

T
PM M PM BHH T BT .

(2.8)

The corrections of M and PM due to the presence of the residual covarianceare estimated as follows: by the first-order perturbation theorem (Griffiths,
2005), the correction of the ith major eigenvalue of Cov[x] is upper-boundedby max eig[BBT ], while the norm of the correction of the ith eigenvector
is upper-bounded by max eig[BBT ]/ min eig[H T BT BH]. These corrections
are negligibly smaller than the leading-order terms, when inequality 2.7
holds. Further details are provided in section 4.4.
Therefore, applying PCA to Cov[x] can extract the signal covarianceas the major principal components with high accuracy when Nx = N f 
Ns  1 and B is sufficiently isotropic. This property is numerically validated. The major principal components of Cov[x] suitably approximate
the signal covariance (see Figure 2A). In the simulations, elements of B
are sampled from a gaussian distribution with zero mean. The trace ratio tr[BBT ]/tr[H T BT BH] (i.e., the ratio of the eigenvalue sum of residual
covariance tr[BBT ] to that of signal covariance tr[H T BT BH]) retains thevalue of about 0.6 irrespective of Ns and Nx , indicating that the largescale systems we consider are fairly nonlinear (see Figure 2B). In contrast,
the eigenvalue ratio max eig[BBT ]/ min eig[H T BT BH] monotonically converges to zero when Nx increases; thus, inequality 2.7 holds for large Nx (seeFigure 2C). Consequently, PM approximates nonzero eigenmodes of the signal covariance accurately (see Figure 2D). These results indicate that PCA isa promising method for reliably identifying the linear components in sensory inputs. In what follows, we explicitly demonstrate that projecting the
inputs onto the major eigenspace can recover true hidden sources of sensory inputs accurately.
2.3 Asymptotic Linearization Theorem. We now consider a linear encoder comprising a single-layer linear neural network,u ≡ W (x − E[x]) ,(2.9)

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

max eig[BT B]/ min eig[BT B] can be upper-bounded by a (possibly large)
finite constant and focus on such B, because the effective input dimensionality is otherwise much smaller than Nx . In other words, one can redefine
Nx of any system by replacing the original sensory inputs with their compressed representation to render B isotropic.
When inequality 2.7 holds, the first Ns major eigenmodes of Cov[x] precisely match the signal covariance, while minor eigenmodes are negligibly
small. Hence, there is a clear spectrum gap between the largest Ns eigenvalues and the rest. This indicates that one can reliably identify the signal
covariance and source dimensionality based exclusively on PCA of Cov[x]in an unsupervised manner. Hence, when M ∈ RNs ×Ns is a diagonal matrix
that arranges the first Ns major eigenvalues of Cov[x] in descending order,
with the corresponding eigenvector matrix PM ∈ RNx ×Ns , we obtain

Achievability of Nonlinear Blind Source Separation

1443

where u ≡ (u1 , . . . , uNs )T are Ns -dimensional neural outputs and W ∈RNs ×Nx is a synaptic weight matrix. Suppose that by applying PCA to the
inputs, one obtains W, which represents a subspace spanned by the maTjor eigenvectors, W = M −1/2
M is an arbitrary Ns × Ns orM PM , where
thogonal matrix expressing an ambiguity. A normalization factor −1/2
M

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Figure 2: PCA can identify linear components comprising linear projections
from every true hidden source. Here, hidden sources s are independently generated from an identical uniform distribution with zero mean and unit variance; N f = Nx and f (•) = sign(•) are fixed; elements of A, a are independently
sampled from N [0, 1/Ns ]; and elements of B are independently sampled fromN [0, 1/N f ]. PCA is achieved via eigenvalue decomposition. (A) Comparison between the signal covariance and major principal components of sensory inputs,
when Ns = 10 and Nx = 103 . (B) Trace ratio tr[BBT ]/tr[H T BT BH]. (C) Eigenvalue ratio max eig[BBT ]/ min eig[H T BT BH]. (D) Error in extracting nonzero
eigenmodes of signal covariance as major components, which is scored byT
ULULT PM ]/Ns , where UL ∈ RNx ×Ns indicates the left-singular vectors of
1 − tr[PMTconverges to ULULT reliably and accurately. Solid
BH. As Nx increases, PM PMlines indicate theoretical values of estimation errors: see equation 4.24 for details. Circles and error bars indicate the means and areas between maximum
and minimum values obtained with 20 different realizations of s, A, B, a, wheresome error bars are hidden by the circles.

1444

T. Isomura and T. Toyoizumi

W(BH)+ .

(2.10)

Here, (BH)+ ≡ (H T BT BH)−1 H T BT is the pseudo inverse of BH, and is another arbitrary Ns × Ns orthogonal matrix. Error in approximating (BH)+
is negligible for the following calculations as long as inequality 2.7 holds(see section 4.4 for more details). Equation 2.10 indicates that the directionsof the linear components can be computed under the BSS setup—up to anarbitrary orthogonal ambiguity factor . Hence, we obtain the following
theorem:
Theorem 1 (Asymptotic Linearization). When inequality 2.7 holds, from equaTtions 2.2, 2.9, and 2.10, the linear encoder with optimal matrix W = M −1/2
M PM
can be analytically expressed asu=(s + ε)(2.11)

using a linearization error ε ≡ (BH)+ Bφ with the covariance matrix ofCov[ε] = (BH)+ BBT (BH)+T .

(2.12)The maximum eigenvalue of Cov[ε] is upper-bounded by max eig[BBT ]/min eig[H T BT BH], which is sufficiently smaller than one from inequality 2.7. In
particular, when ps (s) is a symmetric distribution, f (•) is an odd function, and B
is sufficiently isotropic, using equation 2.6, Cov[ε] can be explicitly computed as
Ns
Cov[ε] =
Nf



f2
f



− 1 (I +
2

)+

f (3)

2

2Ns f

2

I

(2.13)

as the leading order. Symmetric matrix ≡ UAT (BT B − I)2UA characterizes theanisotropy of BT B in the directions of the left-singular vectors of A, UA ∈ RN f ×Ns ,Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

is multiplied so as to ensure Cov[u] = I. Neural outputs with this W indeed express the optimal linear encoder of the inputs because it is thesolution of the maximum likelihood estimation that minimizes the loss
to reconstruct the inputs from lower-dimensional encoder u using a linear network under gaussian assumption (see Xu, 1993, and Wentzell, Andrews, Hamilton, Faber, & Kowalski, 1997, for related studies). In otherwords, when we assume that the loss follows a unit gaussian distribution,Targ minW E[|x − E[x] − W + u|2 ] = M −1/2
M PM holds under the constraint+of dim[u] = Ns and Cov[u] = I, where W indicates the pseudo inverseof W.
Crucially, equation 2.8 directly provides the key analytical expression torepresent a subspace spanned by the linear components:

Achievability of Nonlinear Blind Source Separation

1445

wherein max eig[ ] is upper-bounded by max eig[(BT B − I)2 ] = O(1). Together,
we conclude that

u=

s+O

Ns
Nf





1
.
+O √
Ns

(2.14)

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

The derivation detail of equation 2.13 is provided in section 4.5. Equation
2.13 indicates that the linearization error monotonically decreases as N f /Nsand Ns increase; thus, u converges to a linear mixture of all the true hidden
sources (u → s) in the limit of large N f /Ns and Ns . Only the anisotropyof BT B in the directions of UA increases the linearization error. In essence,
applying PCA to the inputs effectively filters out nonlinear components in
the inputs because the majority of nonlinear components are perpendicular
to the directions of the signal covariance.
Although the obtained encoder is not independent of each other becauseof the multiplication with , it is remarkable that the proposed approachenables the conversion of the original nonlinear BSS problem to a simple
linear BSS problem. This indicates that u can be separated into each independent encoder by further applying a linear ICA method (Comon, 1994;
Bell & Sejnowski, 1995, 1997; Amari et al., 1996; Hyvarinen & Oja, 1997)
to it, and these independent encoders match the true hidden sources up to
permutations and sign-flips. Zero-element-wise-error nonlinear BSS is attained in the limit of large N f /Ns and Ns . As a side note, PCA can indeed
recover a linear transformation of true hidden sources in the inputs evenwhen these sources have higher-order correlations if the average of these
correlations converges to zero with high source dimensionality. This property is potentially useful, for example, for identifying true hidden states oftime series data generated from nonlinear systems (Isomura & Toyoizumi,2021).
In summary, we analytically quantified the accuracy of the optimal linear
encoder—obtained through the cascade of PCA and ICA—in inverting thenonlinear generative process to identify all hidden sources. The encoder
increases its accuracy as N f /Ns and Ns increase and asymptotically attains
the true hidden sources.
The proposed theorem is empirically validated by numerical simulations. Each element of the optimal linear encoder obtained by the PCA–ICA
cascade represents a hidden source, wherein BSS errors (the difference between the true and estimated sources) decrease as Nx increases (see Figure
3A). The PCA–ICA cascade performs the nonlinear BSS with various typesof nonlinear basis functions (see Figures 3B and 3C). This is a remarkableproperty of the PCA–ICA cascade because these results indicate that it canperform the nonlinear BSS without knowing the true nonlinearity that characterizes the generative process. Although equation 2.13 is not applicable
to non-odd nonlinear basis function f (•), empirical observations indicate

1446

T. Isomura and T. Toyoizumi

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Figure 3: BSS error can be characterized by the source and input dimensionalities. In all panels, s is generated by an identical uniform distribution; N f = Nx
is fixed; and A, B, a are sampled from gaussian distributions as in Figure 2. Encoders ũ = (ũ1 , . . . , ũNs )T are obtained by applying the PCA–ICA cascade to sensory inputs. For visualization purpose, elements of ũ are permuted and signflipped to ensure that ũi encodes si . (A) Each element of ũ encodes a hidden
source, wherein an error in representing a source monotonically decreases when
Nx increases. The generative process is characterized with Ns = 100 and f (•) =sign(•). Shaded areas represent theoretical values of the standard deviation
computed using equation 2.13. As elements of B are gaussian distributed, = I
holds. Inset panels depict the absolute value of covariance matrix |Cov[ũ, s]|with gray-scale values ranging from 0 (white) to 1 (black). A diagonal covariance matrix indicates the successful identification of all the true hidden sources.
(B, C) Nonlinear BSS when the generative process is characterized by f (•) = (•)3or f (•) = ReLU(•). ReLU(•) outputs • for • > 0 or 0 otherwise. In panel C, theshaded area is computed using equation 2.12. (D) The quantitative relationshipbetween the source and input dimensionalities and element-wise BSS error isscored by the mean squared error E[|s − ũ|2 ]/Ns . Here, f (•) = sign(•) is supposed. Circles and error bars indicate the means and areas between maximum
and minimum values obtained with 20 different realizations of s, A, B, a, wheresome error bars are hidden by the circles. Solid and dashed lines represent
theoretical values computed using equations 2.12 and 2.13, respectively. Theselines fit the actual BSS errors when Nx  Ns  1, although some deviations
occur when Ns or Nx is small. Blue cross marks indicate actual BSS errors for
Ns = 10 when hidden sources are sampled from a symmetric truncated normal
distribution.

Achievability of Nonlinear Blind Source Separation

1447

2.4 Hebbian-Like Learning Rules Can Reliably Solve Nonlinear BSS
Problem. As a corollary of the proposed theorem, a linear neural network
can identify true hidden sources through Hebbian-like plasticity rules in the
nonlinear BSS setup under consideration. Oja’s subspace rule (Oja, 1989)—
which is a modified version of the Hebbian plasticity rule (see section 4.6)—
is a well-known PCA approach that extracts the major eigenspace withoutyielding a spurious solution or attaining a local minimum (Baldi & Hornik,
1989; Chen et al., 1998). Thus, with generic random initial synaptic weights,
this Hebbian-like learning rule can quickly and reliably identify an optimallinear encoder that can recover true hidden sources from their nonlinear
mixtures in a self-organizing or unsupervised manner.
Numerical experiments demonstrate that regardless of the random initialization of synaptic weight matrix WPCA ∈ RNs ×Nx , Oja’s subspace rule updates WPCA to converge to the major eigenvectors, that is, the directionsof the linear components. The accuracy of extracting the linear components increases as the number of training samples increases, and reaches
the same accuracy as an extraction via eigenvalue decomposition (see Figure 4A). Because the original nonlinear BSS problem has now been transformed to a simple linear BSS problem, the following linear ICA approach(Comon, 1994; Bell & Sejnowski, 1995, 1997; Amari et al., 1996; Hyvarinen& Oja, 1997) can reliably separate all the hidden sources from the features
extracted using Oja’s subspace rule. Amari’s ICA algorithm (Amari et al.,1996), which is another Hebbian-like rule (see section 4.6), updates synaptic
weight matrix WICA ∈ RNs ×Ns to render neural outputs independent of eachother. The obtained independent components accurately match the true
hidden sources of the nonlinear generative process up to their permutations

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

that the PCA–ICA cascade can identify the true hidden sources even with
non-odd function f (•) (see Figure 3C), as long as H = E[ f sT ] is nonzero.
The log-log plot illustrates that the magnitude of element-wise BSS
errors—scored by the mean squared error—decreases inversely proportional to Nx /Ns ; however, it saturates around Nx = Ns2 (see Figure 3D). Theseobservations validate equation 2.13 which asserts that the linearization error is determined by the sum of Ns /N f and 1/Ns order terms. Althoughequation 2.13 overestimates the BSS error when Ns = 10, this is because
each dimension of As significantly deviates from a gaussian variable due
to small Ns —as s is sampled from a uniform distribution in these simulations. We confirm that when hidden sources are generated from a distribution close to gaussian, actual BSS errors shift toward the theoretical valueof equation 2.13 even when Ns = 10 (see the cross marks in Figure 3D). Indeed, this deviation disappears for large Ns according to the central limit
theorem. Therefore, equation 2.13 is a good approximation of actual BSS
errors for Nx = N f  Ns  1, and the proposed theorem suitably predicts
the performance of the PCA–ICA cascade for a wide range of nonlinear BSS
setups.

1448

T. Isomura and T. Toyoizumi

and sign-flips (see Figure 4B). These results highlight that the cascade of
PCA and ICA—implemented via Hebbian-like learning rules—can selforganize the optimal linear encoder and therefore identify all the true hidden sources in this nonlinear BSS setup, with high accuracy and reliability
when Nx = N f  Ns  1.
3 Discussion
In this work, we theoretically quantified the accuracy of nonlinear BSS
performed using the cascade of linear PCA and ICA when sensory inputs are generated from a two-layer nonlinear generative process. First, wedemonstrated that as the dimensionality of hidden sources increases andthe dimensionalities of sensory inputs and nonlinear bases increase relative to the source dimensionality, the first Ns major principal components

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Figure 4: Hebbian-like learning rules can identify the optimal linear encoder.
As in Figures 2 and 3, s is generated by an identical uniform distribution; Ns =100, N f = Nx , and f (•) = sign(•) are fixed; and A, B, a are sampled from gaussian distributions. (A) Learning process of Oja’s subspace rule for PCA. Synaptic weight matrix WPCA ∈ RNs ×Nx —initialized as a random matrix—converges
T
up to multiplication of an orthogonal matrix from the left-hand side.
to PM
This indicates that Oja’s rule extracts the linear components according to the
proposed theorem. Dashed lines indicate the estimation errors when eigenvalue decomposition is applied (see Figure 2C). (B) Learning process of Amari’s
ICA algorithm. Synaptic weight matrix WICA ∈ RNs ×Ns —initialized as a randommatrix—learns to separate the compressed inputs WPCA (x − E[x]) into independent signals. Because PCA yields a linear transformation of hidden sources,
the ensuing independent encoder ũ ≡ WICAWPCA (x − E[x]) identifies all the true
hidden sources with a small BSS error. Elements of ũ are permuted and signflipped to ensure that ũi encodes si . Dashed lines are computed using equation
2.13. In panels A and B, learning rates are ηPCA = 10−3 and ηICA = 0.02, respectively. Solid lines represent the mean estimation errors, while shaded areas represent areas between maximum and minimum values obtained with 20 different
realizations of s, A, B, a.

Achievability of Nonlinear Blind Source Separation

1449

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

approximately express a subspace spanned by the linear projections fromall hidden sources of the sensory inputs. Under the same condition, wethen demonstrated that the optimal linear encoder obtained by projecting
the inputs onto the major eigenspace can accurately recover all the true
hidden sources from their nonlinear mixtures. This property is termed theasymptotic linearization theorem. The accuracy of the subspace extraction
increases as N f /Ns and Ns increase, because the gap between the minimum
eigenvalue of the linear (i.e., signal) components and maximum eigenvalueof the nonlinear (i.e., residual) components becomes significantly large.
Hebbian-like plasticity rules can also identify the optimal linear encoder
by extracting major principal components in a manner equal to PCA. Subsequent application of linear ICA on the extracted principal components
can reliably identify all the true hidden sources up to permutations andsign-flips. Unlike conventional nonlinear BSS methods that can yield spurious solutions (Hyvarinen & Pajunen, 1999; Jutten & Karhunen, 2004), thePCA–ICA cascade is guaranteed to identify the true hidden sources in theasymptotic condition because it successfully satisfies requirements 1 to 5
specified earlier.
The unique identification of true hidden sources s (up to permutations
and sign-flips) is widely recognized only under the linear BSS setup. In
contrast, it is well known that conventional nonlinear BSS approaches using nonlinear neural networks do not guarantee the identification of true
hidden sources under the general nonlinear BSS setup (Hyvarinen & Pajunen, 1999; Jutten & Karhunen, 2004). One may ask if nonlinear BSS methods may find a component-wise nonlinear transformation of the original
sources s = g(s) instead of s because s is still an independent source representation. This is observed when conventional nonlinear BSS approaches
based on the nonlinear ICA are employed because the nonlinear ICA finds
one of many representations that minimize the dependency among outputs.
However, such nonlinear BSS approaches do not guarantee the identification of true sources. This is because when g (s) is a component-wise nonlinear transformation that renders a nongaussian source gaussian distributed,
a transformation s = g(Rg (s))—characterized by some orthogonal matrix
R and component-wise nonlinear transformation g—can yield an arbitrary
independent representation s that differs from the original sources. The
former (s ) is a special case of the latter (s ). Even in the former case, findings transformed by a highly nonlinear, nonmonotonic function g(s) is problematic. Thus, the general nonlinear BSS is essentially an ill-posed problem.
Having said this, earlier works typically investigated the case in which thesources and inputs have the same dimensionality, while the other conditions have not been extensively analyzed. Thus, in the current work, we focused on the nonlinear BSS under the condition where Nx = N f  Ns  1.
A remarkable aspect of the proposed nonlinear BSS approach is that
it utilizes the linear PCA to extract the linear components of the original sources from the mixed sensory inputs. This is the process that

1450

T. Isomura and T. Toyoizumi

4This property can be understood as follows. Because x is generated through x =B f (As + a) with a sufficiently isotropic B that satisfies inequality 2.7, from the proposed−1/2 T
theorem, the encoder u = MPM (x − E[x]) = (s + ε) asymptotically becomes u → s
in the limit of large N f /Ns and Ns . In addition, we define s = g(s) ∈ RNs as a nonlinear
transformation of s. If there is another generative process x = B f (A s + a ) that can generate the observed x while satisfying inequality 2.7—where A and a are gaussian distributed matrix and vector, f is a nonlinear function (not the derivative of f , unlike inthe main text), and B is a sufficiently isotropic matrix—the encoder can also be expressed−1/2 T
as u = MPM (x − E[x]) = (s + ε ), which asymptotically becomes u → s in thelimit of large N f /Ns and Ns . Note that M and PM are the same as above because x does
not change. However, while u ≡ u holds by construction, s ≡ s for any orthogonal
matrices andbecause s is a nonlinear transformation of s. Thus, such a generative
process x = B f (A s + a ) does not exist in the limit of large N f /Ns and Ns .
5
It should be noted that unlike the condition we considered (i.e., Nx = N f  Ns  1),
when N f  Nx , one may observe a counterexample of the proposed theorem such thatthe major principal components are dissimilar to a linear transformation of the true hidden sources, as is well-known in the literature on nonlinear ICA (Hyvarinen & Pajunen,

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

enables the reduction of the original nonlinear BSS to a simple linear BSS
and consequently ensures the reliable identification of the true sources. In
other words, the proposed approach does not rely on the nonlinear ICA;
thus, no concerns exist regarding the creation of the already noted spurious
solutions. We mathematically demonstrated the absence of such spurious solutions when Nx = N f  Ns  1. Specifically, we adopted a minimal
assumption about the relationship between B and the system dimensionalities such that N f /Ns and Ns are much greater than max eig[BT B]/
min eig[BT B]. Our mathematical analyses demonstrate that this condition is
sufficient to asymptotically determine the possible form of hidden sources,
which can generate the observed sensory inputs x, in a unique manner. In
contrast, in the limit of large N f /Ns and Ns , no nonlinear transformation ofs (i.e., s or s ) can generate the observed sensory inputs while satisfying the
aforementioned condition.4 Thus, when inequality 2.7 holds, the true hidden sources s are uniquely determined up to permutations and sign-flips
ambiguities without component-wise nonlinear ambiguity. (Conversely, if
inequality 2.7 does not hold, it might not be possible to distinguish true hidden sources and their nonlinear transformations in an unsupervised manner.) Hence, under the condition we consider, the nonlinear BSS is formallyreduced to a simple linear BSS, wherein only permutations and sign-flips
ambiguities exist while no nonlinear ambiguity remains. This property is
crucial for identifying the true sources under the nonlinear BSS setup, without being attracted by spurious solutions such as s and s .
The nonlinear generative processes that we considered in this work aresufficiently generic. Owing to the universality of two-layer networks, eachelement of an arbitrary generative process x = F(s) can be approximatedusing equation 2.1 with a high degree of accuracy in the component-wise
sense when Nx = N f  Ns  1.5 This allows one to analyze the achievability of the nonlinear BSS using a class of generative processes comprisingAchievability of Nonlinear Blind Source Separation

1451

1999; Jutten & Karhunen, 2004). This is because when N f  Nx , a general B may not be
sufficiently isotropic. In other words, N f  Nx corresponds to the case wherein Nx is insufficient to ensure inequality, equation 2.7; thus, Nx needs to be greater. Hence, it is remarkable that the proposed theorem specifies the condition under which the achievabilityof the nonlinear BSS is mathematically guaranteed.

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

random basis functions, as a proxy for investigating x = F(s). As a general property of such generative processes, the eigenvalues of the signal
covariance (i.e., linear components) become significantly larger than thoseof the residual covariance (i.e., nonlinear components) when Nx = N f 
Ns  1. Sufficient system dimensionalities to exhibit this property are determined depending on matrix B that characterizes the generative process. Namely, this observation holds true at least when both N f /Ns and Ns
are much greater than max eig[BT B]/ min eig[BT B]. The existence of suchNs and N f is guaranteed if we consider sufficiently isotropic B wherein
max eig[BT B]/ min eig[BT B] is upper-bounded by a large, finite constant.
Meanwhile, with such a large constant, the two-layer networks can approximate each element of x = F(s) accurately. This in turn means that as N f /Nsand Ns increase, the proposed theorem becomes applicable to a sufficiently
broad class of nonlinear generative processes.
This asymptotic property can be understood as follows. When adding
a new random basis ( fk = Hk s + φk ) to the existing large-scale system,
the linear component of fk (Hk s) must be placed within the existing lowdimensional subspace spanned only by Ns ( N f ) linear projections ofsources. In contrast, its nonlinear component (φk ) is almost uncorrelatedwith other nonlinear components (φ), as shown in equation 2.6, because
these components are characterized by gaussian distributed A. Thus, theformer increases the eigenvalues of the signal covariance in proportion to
N f /Ns , while the latter adds a new dimension in the subspace of nonlinear components without significantly increasing the maximum eigenvalueof the residual covariance. Owing to this mechanism, the eigenvalues ofthe signal covariance become N f /Ns times larger than those of the residual
covariance when Ns  1. Hence, the linear PCA is sufficient to reduce thenonlinear BSS to a linear BSS.
Nonlinear variants of PCA such as autoencoders (Hinton & Salakhutdinov, 2006; Kingma & Welling, 2013) have been widely used for representation learning. Because natural sensory data are highly redundantand do not uniformly cover the entire input space (Chandler & Field,2007), finding a concise representation of the sensory data is essential tocharacterize its properties (Arora & Risteski 2017). In general, if a large,
nonlinear neural network is used, many equally good solutions will be produced (Kawaguchi, 2016; Lu & Kawaguchi, 2017; Nguyen & Hein, 2017).
In this case, there is no objective reason to select one solution over another if they have similar reconstruction accuracy. However, this property also leads to infinitely many spurious solutions if the aim is the1452

T. Isomura and T. Toyoizumi

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

identification of the true hidden sources (Hyvarinen & Pajunen, 1999; Jutten & Karhunen, 2004). Consequently, the outcomes of these approaches
using nonlinear neural networks are intrinsically ambiguous, and obtained
solutions highly depend on the heuristic design of the regularization parameters used in the networks and learning algorithms (Dahl, Yu, Deng,& Acero, 2012; Hinton, Srivastava, Krizhevsky, Sutskever, & Salakhutdinov, 2012; Wan, Zeiler, Zhang, LeCun, & Fergus, 2013). However, unlike
those nonlinear approaches, we proved that the cascade of linear PCA andICA ensures that the true hidden sources of sensory inputs are obtained if
there are sufficient dimensionalities in mappings from hidden sources tosensory inputs. Thus, this linear approach is suitable to solve the nonlinear BSS problem while retaining the guarantee to identify the true solution.
Its BSS error converges in proportion with the source-to-input dimensionality ratio for a large-scale system comprising high-dimensional sources. Furthermore, if needed, its outcomes can be used as a plausible initial condition
for nonlinear methods to further learn the forward model of the generative
process.
Because most natural data (including biological, chemical, and social
data) are generated from nonlinear generative processes, the scope of application of nonlinear BSS is quite broad. The proposed theorem provides atheoretical justification for the application of standard linear PCA and ICA
to nonlinearly generated natural data to infer their true hidden sources.
An essential lesson from the proposed theorem is that an artificial intelligence should employ a sufficient number of sensors to identify true hidden
sources of sensory inputs accurately, because the BSS error is proportional
to the source-to-input dimensionality ratio in a large-scale system.
In the case of computer vision, low-dimensional representations (typically 10–104 components) are extracted from up to millions of pixels of thehigh-dimensional raw images. Our theory implies that the linear PCA–ICA
cascade can be utilized to identify the true hidden sources of natural image data. This has been demonstrated using video images (Isomura & Toyoizumi, 2021). In particular, a combination of the proposed theorem withanother scheme that effectively removes unpredictable noise from sequential data (i.e., video sequence) is a powerful tool to identify the true hidden
sources of real-world data. By this combination, one can obtain data withminimal irrelevant or noise components—an ideal condition for the asymptotic linearization theorem—from the original noisy observations. Such a
combination enables the reliable and accurate estimations of hidden sourcesthat generate the input sequences, even in the presence of considerable observation noise. The source estimation performance was demonstrated byusing sequential visual inputs comprising hand-digits, rotating 3D objects,
and natural scenes. These results support the applicability of the proposed
theorem to the natural data and further highlight that the PCA–ICA cascade can extract true or relevant hidden sources when the natural data aresufficiently high-dimensional. (See Isomura & Toyoizumi, 2021, for details.)

Achievability of Nonlinear Blind Source Separation

1453

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

In addition, a few works have cited a preprint version of this article tojustify the application of their linear methods to real-world nonlinear BSS
tasks, in the contexts of wireless sensor networks (van der Lee, Exarchakos,& de Groot, 2019) and single-channel source separation (Mika, Budzik, &Józwik, 2020). These works demonstrated that linear methods work wellwith nonlinear BSS tasks, as our theorem predicts. The asymptotic linearization theorem is of great importance in justifying the application of theselinear methods to real-world nonlinear BSS tasks.
In terms of the interpretability of the outcomes of the PCA–ICA cascade, one may ask if applying the cascade to natural image patches simply yields less interpretable Gabor-filter-like outputs, which are usually notconsidered to be certain nonlinear, higher-order features underlying natural
images. However, when we applied the PCA–ICA cascade—featuring aseparate noise-reduction technique that we developed—to natural image
data with a high-dimensionality reduction rate, we obtained outputs or images that represent features relevant to hidden sources, such as the categories of objects (Isomura & Toyoizumi, 2021). Based on our observations,
three requirements to render the PCA–ICA cascade extract features relevant to hidden sources are considered. First, the dimensionality reduction
rate should be high relative to that of the conventional PCA application. Indeed, when the dimensionality reduction rate was low, we observed that theobtained encoders exhibit Gabor-filter-like patterns. Second, the PCA–ICA
cascade should be applied not to image patches but to entire images. Image
patches contain only a fraction of information about the original hidden
sources; thus, it is likely that the original hidden sources cannot be recovered from image patches. Third, PCA should be applied to denoised data
wherein observation noise is removed in advance using a separate scheme,
as in our case. When PCA is applied to noisy data, sufficient information forrecovering hidden sources cannot be extracted as the major principal components because PCA preferentially extracts large noise owing to its extra
variance. Hence, when the requirements are satisfied, the outcome of the
PCA–ICA cascade can exhibit a high interpretability, wherein the extracted
features are relevant to hidden sources.
As a side note, to utilize the asymptotic linearization theorem to estimate
the accuracy of source separation, one cannot use up-sampling techniques
to increase the input dimensionality beyond the resolution of the original data or images. This is because these techniques simply yield linear or
nonlinear transformations of the original data, which typically do not increase information about true hidden sources. In other words, the resolutions of the original data determine the effective input dimensionality, sothat these up-samplings merely render matrix B singular (thereby violating
the condition of isotropic B).
Furthermore, an interesting possibility is that living organisms might
have developed high-dimensional biological sensors to perform nonlinear BSS using a linear biological encoder because this strategy guarantees

1454

T. Isomura and T. Toyoizumi

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

robust nonlinear BSS. In particular, this might be the approach taken by thehuman nervous system using large numbers of sensory cells, such as approximately 100 million rod cells and 6 million cone cells in the retina, andapproximately 16,000 hair cells in the cochlea, to process sensory signals(Kandel, Schwartz, Jessell, Siegelbaum, & Hudspeth, 2013).
Neuronal networks in the brain are known to update their synapses byfollowing Hebbian plasticity rules (Hebb, 1949; Malenka & Bear, 2004), andresearchers believe that the Hebbian plasticity plays a key role in representation learning (Dayan & Abbott, 2001; Gerstner & Kistler, 2002) and BSS(Brown, Yamada, & Sejnowski, 2001) in the brain—in particular, the dynamics of neural activity and plasticity in a class of canonical neural networks can be universally characterized in terms of representation learningor BSS from a Bayesian perspective (Isomura & Friston, 2020). Thus, it follows that major linear BSS algorithms for PCA (Oja, 1982, 1989) and ICA(Bell & Sejnowski, 1995, 1997; Amari et al., 1996; Hyvarinen & Oja, 1997) areformulated as a variant of Hebbian plasticity rules. Moreover, in vitro neuronal networks learn to separately represent independent hidden sources in(and only in) the presence of Hebbian plasticity (Isomura, Kotani, & Jimbo,
2015; Isomura & Friston, 2018). Nonetheless, the manner in which the brain
can possibly solve a nonlinear BSS problem remains unclear, even though
it might be a prerequisite for many of its cognitive processes such as visual recognition (DiCarlo et al., 2012). While Oja’s subspace rule for PCA(Oja, 1989) and Amari’s ICA algorithm (Amari et al., 1996) were used in thisarticle, these rules can be replaced with more biologically plausible local
Hebbian learning rules (Foldiak, 1990; Linsker, 1997; Pehlevan, Mohan, &Chklovskii, 2017; Isomura & Toyoizumi, 2016, 2018, 2019; Leugering & Pipa,2018) that require only directly accessible signals to update synapses. A recent work indicated that even a single-layer neural network can perform
both PCA and ICA through a local learning rule (Isomura & Toyoizumi,2018), implying that even a single-layer network can perform a nonlinear
BSS.
In summary, we demonstrated that when the source dimensionality islarge and the input dimensionality is sufficiently larger than the source dimensionality, a cascaded setup of linear PCA and ICA can reliably identify
the optimal linear encoder to decompose nonlinearly generated sensory inputs into their true hidden sources with increasing accuracy. This is because
the higher-dimensional inputs can provide greater evidence about the hidden sources, which removes the possibility of finding spurious solutions
and reduces the BSS error. As a corollary of the proposed theorem, Hebbianlike plasticity rules that perform PCA and ICA can reliably update synaptic weights to express the optimal encoding matrix and can consequently
identify the true hidden sources in a self-organizing or unsupervised manner. This theoretical justification is potentially useful for designing reliable
and explainable artificial intelligence, and for understanding the neuronalmechanisms underlying perceptual inference.

Achievability of Nonlinear Blind Source Separation

1455

4 Methods
4.1 Lemma 1.

Lemma 1. When F(ỹ) is order 1, the difference between E[F(ỹ)] and EN [F(ỹ)] is
upper-bounded by order Ns−1 as a corollary of the central limit theorem. In particular, E[F(ỹ)] can be expressed as EN [F(ỹ)(1 + G(ỹ)/Ns )], where G(ỹ) = O(1) isa function that characterizes the deviation of p(ỹ) from pN (ỹ).
T

Proof. The characteristic function of p(ỹ) is given as ψ (τ ) ≡ E[eiτ ỹ ] as afunction of τ ∈ R2 , and its logarithm is denoted as (τ ) ≡ log ψ (τ ). The
first- to fourth-order derivatives of (τ ) with respect to τ are providedasT (1) = E[iỹeiτ ỹ ]/ψ,
 (2) = E[(iỹ)⊗2 eiτ ỹ ]/ψ − ( (1) )⊗2 ,T

 (3) = E[(iỹ)⊗3 eiτ ỹ ]/ψ − E[(iỹ)⊗2 eiτ ỹ ]/ψ ⊗  (1) −  (2) ⊗  (1)TT−  (1) ⊗  (2)= E[(iỹ)⊗3 eiτ ỹ ]/ψ − 2 (2) ⊗  (1) − ( (1) )⊗3 −  (1) ⊗  (2) ,T (4) = E[(iỹ)⊗4 eiτ ỹ ]/ψ − E[(iỹ)⊗3 eiτ ỹ ]/ψ ⊗  (1) − 2 (3) ⊗  (1)TT− 3( (2) )⊗2 −  (2) ⊗ ( (1) )⊗2 −  (1) ⊗  (2) ⊗  (1)− ( (1) )⊗2 ⊗  (2) −  (1) ⊗  (3) .

(4.1)

When τ = 0, they become (0) = 0,  (1) (0) = 0,  (2) (0) = −E[ỹ⊗2 ], (3) (0) = 0, and  (4) (0) = E[ỹ⊗4 ] − 3E[ỹ⊗2 ]⊗2 owing to the symmetricsource distribution ps (s).
They are further computed as E[ỹ⊗2 ]=Vec[ÃÃT ] and E[ỹ⊗4 ]=Vec[E[(ỹỹT )
 s
T ⊗2
4
⊗ (ỹỹT )]] = 3Vec[ÃÃT ]⊗2 + κVec[ N
k=1 (Ã•k Ã•k ) ], where κ ≡ E[sk ] − 3 indicates the kurtosis of the source distribution and Ã•k denotes the kth
⊗2
T
T
(4)
⊗4
column (2) (0)Nsof Ã. TMoreover,N·s τ T = −τ4 ÃÃ τ Tand2  (0) · τ−3/2=⊗2⊗4κVec[ k=1 (Ã•k Ã•k ) ] · τ= κ k=1 (τ Ã•k ) = 3κ (τ τ ) /Ns + O(Ns )
 s T
4
−1
becausehold. Here, a quartic function of τ , N
k=1 (τ Ã•k ) , is order Ns

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Definition 1. Suppose ps (s) is a symmetric distribution, ỹ ≡ (yi , y j )T ≡ Ãs is a
two-dimensional subvector of y = As, and Ã ≡ (Ai1 , . . . , AiNs ; A j1 , . . . , A jNs ) is
a 2 × Ns submatrix of A. The expectation of an arbitrary function F(ỹ) over ps (s)is denoted as E[F(ỹ)], whereas when ỹ is sampled from a gaussian distributionTpN (ỹ) ≡ N [0, ÃÃ ], the expectation of F(ỹ) over pN (ỹ) is denoted as EN [F(ỹ)] ≡F(ỹ)pN (ỹ)dỹ.

1456

T. Isomura and T. ToyoizumiÃ•k is sampled from N [0, 1/Ns ]. Thus, the Taylor expansion of (τ ) isexpressed as∞

(τ ) =n=0


 (n) (0) ⊗n
1κ(τ T τ )2 + O Ns−3/2 .
· τ = − τ T ÃÃT τ +
n!
2
8Ns

(4.2)

where λ denotes the Lebesgue measure. Here, G(ỹ) = O(1) indicates a function of ỹ that characterizes the difference between p(ỹ) and pN (ỹ). Hence,the expectation over p(ỹ) can be rewritten using that over pN (ỹ):
E[F(ỹ)] = F(ỹ)ps (s)ds = F(ỹ)p(ỹ)dỹ



G(ỹ)
= F(ỹ)pN (ỹ) 1 +
dỹ
Ns



G(ỹ)
.
= EN F(ỹ) 1 +
Ns


(4.4)



Remark 1. For equation 2.4, from the Bussgang theorem (Bussgang, 1952),
we obtain



G(ỹ)
E[ f (ỹ)ỹT ] = EN f (ỹ)ỹT 1 +
Ns




G (ỹ)TG(ỹ)+ f (ỹ)ÃÃT= EN diag[ f (ỹ)] 1 +
Ns
Ns


= diag[EN [ f (ỹ)]]Ã + O Ns−3/2 ÃT(4.5)
as Ã is order Ns−1/2 . Thus, we obtain H = E[ f yT ]A+T = diag[EN [ f ]]A +O(Ns−3/2 ).
For equation 2.5, because yi and y j are only weakly correlated witheach other, (1 + G(ỹ)/Ns ) can be approximated as (1 + G(ỹ)/Ns )  (1 +Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Hence, from ψ (τ ) = exp[−τ T ÃÃT τ /2 + κ (τ T τ )2 /8Ns + O(Ns−3/2 )] =T
Te−τ ÃÃ τ /2 (1 + κ (τ T τ )2 /8Ns + O(Ns−3/2 )), p(ỹ) is computed as follows (theinversion theorem):1Te−iτ ỹ ψ (τ )λ(dτ )p(ỹ) =
(2π )2 R2 −3/2κ
1−iτ T ỹ− 12 τ T ÃÃT τ
T 2
1
+
e
(τ
τ
)
+
O
N
λ(dτ )=s(2π )2 R2
8Ns


G(ỹ)
,
(4.3)
= pN (ỹ) 1 +NsAchievability of Nonlinear Blind Source Separation

1457

G(yi )/Ns )(1 + G(y j )/Ns ) as the leading order using newly defined functions
G(yi ) and G(y j ). Thus, we obtain
Cov[φi , φ j ] = E[(φi − E[φi ])(φ j − E[φ j ])]





G(y j )G(yi )
 EN (φi − E[φi ]) 1 +(φ j − E[φ j ]) 1 +
Ns
Ns

∞

=

EN φi∗(n) EN φ ∗(n)
j
n!

n=1

CovN [yi , y j ]n

(4.6)

for i = j as the leading order, where φi∗ ≡ (φi − E[φi ])(1 + G(yi )/Ns ) andφ ∗j ≡ (φ j − E[φ j ])(1 + G(y j )/Ns ) are modified nonlinear components. The
last equality holds according to lemma 2 (see below), wherein EN [φi∗(n) ] approximates E[φi(n) ]. As a corollary of lemma 1, the deviation of EN [φi∗(n) ] 
E[φi(n) ] from EN [φi(n) ] due to the nongaussianity of yi is order Ns−1 . Because CovN [yi , y j ]n = O(Ns−n/2 ), EN [φi(1) ] = O(Ns−1 ), EN [φi(2) ] = O(Ns−1/2 ),and |EN [φi(n) ]| ≤ O(1) for n ≥ 3 hold with a generic odd function f (•) (seeremark 2 for details), the difference between CovN [φi∗ , φ ∗j ] and CovN [φi , φ j ]is order Ns−5/2 . Thus, Cov[φi , φ j ] can be characterized by computing
CovN [φi , φ j ] as a proxy up to the negligible deviation; consequently, equation 2.6 is obtained.
4.2 Lemma 2.
Lemma 2. Suppose v and w are gaussian variables with zero mean, σv and σw
are their standard deviations, and their correlation c ≡ Cov[v, w]/σv σw is smaller
than 1. For an arbitrary function g(•),∞

Cov[g(v ), g(w)] =σvn σwn E[g(n) (v )]E[g(n) (w)]
n=1

cn
.
n!

(4.7)

Proof. Because√ the correlation between v and w is c, w can be cast as w =σw cv/σv + σw 1 − c2 ξ using a new a zero-mean and unit-variance gaussianvariable ξ that is independent of v. When we define the covariance betweeng(v ) and g(w) as (c) ≡ Cov[g(v ), g(w)], its derivative with respect to c isgiven as





σw v
σw cv
σw cξ
2
 (c) = Cov g(v ), g
+ σw 1 − c ξ−√σvσv
1 − c2






σw v
σw cv
σw cξ
+ σw 1 − c2 ξ−√
= E (g(v ) − E[g(v )])g.
σv
σv
1 − c2(4.8)

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

= CovN [φi∗ , φ ∗j ]

1458

T. Isomura and T. Toyoizumi

By applying the Bussgang theorem (Bussgang, 1952) with respect to v, we
obtain(4.9)

Similarly, applying the same theorem with respect to ξ yields


 

σw cv
2
+ σw 1 − c ξ ξ
E (g(v ) − E[g(v )])g
σv



 

σw cv
= E (g(v ) − E[g(v )])g
+ σw 1 − c2 ξ σw 1 − c2 .
σv
(4.10)

As equation 4.8 comprises equations 4.9 and 4.10,  (c) becomes




σw cv
 (c) = σv σw E g (v )g+ σw 1 − c2 ξ .
σv(4.11)

Hence, for an arbitrary natural number n, we obtain




σw cv
+ σw 1 − c2 ξ .
(n) (c) = σvn σwn E g(n) (v )g(n)σv(4.12)

When c = 0, (n) (c) = σvn σwn E[g(n) (v )]E[g(n) (σw ξ )] = σvn σwn E[g(n) (v )]E[g(n) (w)]holds, where we again regard σw ξ = w. Thus, from the Taylor expansion
with respect to c, we obtain∞

(c) =

(n) (0)
n=1

cn
=
n!

∞

σvn σwn E[g(n) (v )]E[g(n) (w)]

cn
.
n!

+T

+

n=1

(4.13)


Remark 2. For φ(y) = f (y + a) − E[ f ] − E[ f y ]A A y, the coefficients EN [φi(n) ] can be computed as follows (up to the fourth order): because applying the Bussgang theorem (Bussgang, 1952)
yields EN [ f yT ]A+T A+ = diag[EN [ f ]] = O(1), the difference between
E[ f yT ]A+T A+ and diag[EN [ f ]] is order Ns−1 (due to the nongaussianT

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021


 


σw cv
E (g(v ) − E[g(v )])g
+ σw 1 − c2 ξ v
σv




σw cv
= E g (v )g+ σw 1 − c2 ξ
σv




σw c 2σw cv+ (g(v ) − E[g(v )])g
+ σw 1 − c2 ξ
σv .
σv
σv

Achievability of Nonlinear Blind Source Separation

1459

contributions of y; see lemma 1). Thus, from φ (1) = diag[ f ] −E[ f yT ]A+T A+ , we have
EN [φ (1) ] = EN diag[ f ] − diag[EN [ f ]] + O(Ns−1 ) = O(Ns−1 ).

(4.14)

Moreover, for the second- to fourth-order derivatives, we obtain

EN [φi(3) ] = EN [ f (3) (yi + ai )] = EN [ f (3) (yi )] + O(Ns−1 ) = O(1),EN [φi(4) ] = EN [ f (4) (yi + ai )] = O(Ns−1/2 ).

(4.15)

Here, the Taylor expansion with respect to ai is applied, where
EN [ f (2) (yi )] = EN [ f (4) (yi )] = 0 owing to odd function f (•). For n ≥ 5,
|EN [φi(n) ]| ≤ O(1) holds. Thus, because CovN [yi , y j ]n = ((AAT ) n )i j is orderNs−n/2 , equation 2.5 becomes
1
(2)
(2)
T
T 2
CovN [φi , φ j ] = EN [φi(1) ]EN [φ (1)j ](AA )i j + EN [φi ]EN [φ j ]((AA ) )i j2O(Ns−5/2 )1
T 3
EN [φi(3) ]EN [φ (3)j ]((AA ) )i j
6
1
T 4
−5/2
+
)
EN [φi(4) ]EN [φ (4)j ]((AA ) )i j + O(Ns
24



+

O(Ns−3 )= EN [ f (3) (yi )]EN [ f (3) (y j )]


1
1
×
ai a j ((AAT ) 2 )i j + ((AAT ) 3 )i j + O(Ns−5/2 ).
2
6(4.16)

4.3 Lemma 3.
Lemma 3. When elements of A ∈ RN f ×Ns are independently sampled from a gaussian distribution N [0, 1/Ns ] and N f > Ns2  1, (AAT ) 3 ≡ (AAT ) (AAT )(AAT ) has Ns major eigenmodes—all of which have the eigenvalue of 3N f /Ns2 (asthe leading order) with the corresponding eigenvectors that match the directionsof AAT —and up to (N f − Ns ) randomly distributed nonzero minor eigenmodes
whose eigenvalues are order max[1, N f /Ns3 ] on average. Thus, the latter are negligibly smaller than the former when N f > Ns2  1.

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

EN [φi(2) ] = EN [ f (2) (yi + ai )] = EN [ f (3) (yi )]ai + O(Ns−3/2 ) = O(Ns−1/2 ),1460

T. Isomura and T. Toyoizumi

Proof. Suppose N f > Ns2  1. When the kth column vector of A is denoted
as A•k ≡ (A1k , . . . , AN f k )T ∈ RN f , (AAT ) 3 becomes
Ns
T

(AA )

3

Ns

Ns

=

(A•k

A•l

A•m )T(A•k2A•l )(A•k2A•l )TA•m )(A•k

A•l

A•l

A•m )(A•kk=1 l=1 m=1
Ns(A•k3 )(A•k3 )T + 3
k=1

+6

(A•k

A•l

k=l

A•m )T .

(4.17)

k<l<m

Ns

3
3 T
2
3
3 T
Here,
A•l )(A•k2 A•l )T
k=l (A•k
k=1 (A•k )(A•k ) = (A )(A ) and= {(A 2 )(A 2 )T } (AAT ) − (A 3 )(A 3 )T hold. Let QDQT be the eigenvalue
decomposition of (A 2 )(A 2 )T using a diagonal matrix that arranges eigenvalues in descending order D ∈ RNs ×Ns with the corresponding eigenvectors Q ∈ RN f ×Ns . We thus have an expression {(A 2 )(A 2 )T } (AAT ) =
Ns
1/2
1/2
T
k=1 diag(Dkk V•k )AA diag(Dkk V•k ). Thus, we obtain
Ns
1/2
T
3
3 T
diag(D1/2
kk V•k )AA diag(Dkk V•k ) − 2(A )(A )(AAT ) 3 = 3
k=1

+6

(A•k

A•l

A•m )(A•kA•lA•m )T .

(4.18)k<l<m
1/2
THence, the first entry in the first term, 3diag[D1/2
11 V•1 ]AA diag[D11 V•1 ] =
−1/2
−1/2
T
(3/Ns )diag[1 + O(Ns )]AA diag[1 + O(Ns )], corresponds to Ns largest
eigenmodes whose eigenvalues are 3N f /Ns2 (as the leading order) with thecorresponding eigenvectors that match the directions of AAT . This is because elements of A are independent and identically distributed variables
sampled from a gaussian distribution N [0, 1/Ns ], whereas the other entries(with k ≥ 2) in the first term approximately yield Ns (Ns − 1) different minor
eigenmodes whose eigenvalues are order N f /Ns3 . Then all the eigenmodesof the second term, −2(A 3 )(A 3 )T , can be expressed using eigenvectors involved in the first term, where all the eigenvalues of (A 3 )(A 3 )T are orderN f /Ns3 . Finally, the third term involves Ns (Ns − 1)(Ns − 2)/6 terms that arenearly uncorrelated with each other—which are also nearly uncorrelatedwith the first and second terms—leading to up to Ns (Ns − 1)(Ns − 2)/6 different minor eigenmodes whose eigenvalues are order N f /Ns3 on average.
Intuitively, A•k A•l A•m with k < l < m approximates each eigenmode
as elements of A are gaussian distributed. Ensure that A•k A•lA•mis nearly uncorrelated with A•k A•l A•m if at least one of (k , l , m )is different from (k, l, m). In summary, when Ns2 + Ns (Ns − 1)(Ns − 2)/6 <Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021=

Achievability of Nonlinear Blind Source Separation

1461

N f , (AAT ) 3 has (Ns (Ns − 1) + Ns (Ns − 1)(Ns − 2)/6) randomly distributed
nonzero minor eigenmodes whose eigenvalues are of order N f /Ns3 on average; otherwise, (AAT ) 3 has (N f − Ns ) randomly distributed nonzero minor
eigenmodes whose eigenvalues are of order 1 on average.

Remark 3. Because the nonzero minor eigenmodes are nearly uncorrelatedwith eigenmodes of AAT , their contributions are negligible when deriving
equation 2.13. It is straightforward to apply this proof to (AAT ) 2 .


BB = (UL , UN )
T

XLL

T
XNL

XNL

XNN



ULT
UNT


.

(4.19)

Here, XLL ≡ ULT BBT UL ∈ RNs ×Ns is a symmetric matrix and XNL ≡UNT BBT UL ∈ R(Nx −Ns )×Ns is a vertically long rectangular matrix. This decomposition separates BBT into components in the directions of UL andthe rest. From equations 2.3 and 4.19, we obtain the key equality thatlinks the eigenvalue decomposition with the decomposition into signal andresidual covariances:(PM , Pm )


M

O

O

m



T
PM

PmT


= (UL , UN )


S2L + XLL

T
XNL

XNL

XNN



ULT
UNT


.

(4.20)The left-hand side is the eigenvalue decomposition of Cov[x], where M ∈RNs ×Ns and m ∈ R(Nx −Ns )×(Nx −Ns ) are diagonal matrices of major and minor
eigenvalues, respectively, and PM ∈ RNx ×Ns and Pm ∈ RNx ×(Nx −Ns ) are theircorresponding eigenvectors. When inequality 2.7 holds, S2L + XLL is sufficiently close to a diagonal matrix and XNL is negligibly small relative
to S2L + XLL ; thus, the right-hand side can be viewed as the zeroth-order

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

4.4 First-Order Approximation of Major Principal Components. The
first Ns major eigenmodes of the input covariance approximates the signal covariance comprising the linear projections from every hidden source.
We first demonstrate this as the zeroth-order approximation by explicitly
decomposing the input covariance (see equation 2.3) into the product oforthogonal and diagonal matrices. We define the singular value decomposition of the linear components BH as BH = UL SLVLT using orthogonal matrices UL ∈ RNx ×Ns and VL ∈ RNs ×Ns and a diagonal matrix SL ∈ RNs ×Ns that
arranges singular values in descending order (note that Nx > Ns ). Moreover, we define an orthogonal matrix UN ∈ RNx ×(Nx −Ns ) such that it is perpendicular to UL and that multiplying it by the residual covariance BBT
from both sides diagonalizes BBT . Thus, (UL , UN )(UL , UN )T = I holds, and
XNN ≡ UNT BBT UN ∈ R(Nx −Ns )×(Nx −Ns ) is a diagonal matrix that arranges its
diagonal elements in descending order. Without loss of generality, BBT isdecomposed into three matrix factors:

1462

T. Isomura and T. Toyoizumi



=

I

ET

−E

I



S2L + XLL

T
XNL

XNL

XNN



I

−E T

E

I



T
S2L + XLL + E T XNL + XNL
E

T−(S2L + XLL )E T + XNL+ E T XNN

−E(S2L + XLL ) + XNL + XNN E

T−XNL E T − EXNL
+ XNN

+ O(|E|2 )(4.21)

up to the first order of E. Equation 4.21 is diagonalized only when −E(S2L +XLL ) + XNL + XNN E = O. Thus, we obtain−1
Vec[E] = I ⊗ (S2L + XLL ) − XNN ⊗ I
Vec[XNL ] ≈ I ⊗ S−2
Vec[XNL ]L⇐⇒ E ≈ XNL S−2
L(4.22)

as the first-order approximation. Substituting E = XNL S−2
L into equation
4.21 yields the first-order approximation of major and minor eigenvalues:2M = S2L I + S−2
L diag[XLL ] + O(|E| ) ,T2
m = XNN − 2diag[XNL S−2
L XNL ] + O(|E| ).

(4.23)

Moreover, the first-order approximation of their corresponding eigenvectors is provided as(PM , Pm ) = (UL , UN )

I

−E T



E
I

= UL + UN E, UN − UL E T
−2 T
2= UL + UN XNL S−2
L , UN − UL SL XNL + O(|E| ).

(4.24)

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

approximation of the eigenvalue decomposition of Cov[x]. Owing to theuniqueness of eigenvalue decomposition, we obtain (PM , Pm ) = (UL , UN ),M = S2L + diag[XLL ], and m = XNN as the leading order.
Next, we compute an error in the first order by considering a small perturbation that diagonalizes the right-hand side of equation 4.20 more accurately. Suppose (I, E T ; −E, I) ∈ RNx ×Nx as a block matrix that comprises Ns ×
Ns and (Nx − Ns ) × (Nx − Ns ) identity matrices and a rectangular matrix
E ∈ R(Nx −Ns )×Ns with small value elements. As (I, E T ; −E, I)(I, E T ; −E, I)T =I + O(|E|2 ), it is an orthogonal matrix as the first-order approximation. Because the middle matrix on the right-hand side of equation 4.20 has been
almost diagonalized, a small perturbation by (I, E T ; −E, I) can diagonalize
it with higher accuracy. By multiplying (I, E T ; −E, I) and its transpose by
the near-diagonal matrix from both sides, we obtain

Achievability of Nonlinear Blind Source Separation

1463

The accuracy of these approximations is empirically validated using numerical simulations (see Figure 2D), wherein an error in approximating
T
T
UL using PM is analytically computed as tr[E T E]/Ns = tr[S−2
L UL BB (I −−2
T
T
ULUL )BB UL SL ]/Ns as the leading order. This error is negligibly small
relative to the linearization error in equation 2.12 when inequality 2.7
holds.

Cov[ε] =

f2
f

2


+

+T

− 1 (BA) BB (BA)
T

+

f (3)

2

2Ns f

2

I

(4.25)

as the leading order. Here, aaT and  in equation 2.6 are negligibly smaller
than the leading-order eigenmodes when inequality 2.7 holds and B is sufficiently isotropic. If we define the left-singular vectors of A as UA ∈ RN f ×Ns ,because A has Ns singular values with the leading-order term of (N f /Ns )1/2
(Marchenko & Pastur, 1967), A = (N f /Ns )1/2UA holds approximately. Thus,(BA)+ = (Ns /N f )1/2 (UAT BT BUA )−1UAT BT holds as the leading order.
Further, as B is isotropic, UAT BT BUA is sufficiently close to the identity matrix; thus, (UAT BT BUA )−1 ≈ I − UAT (BT B − I)UA provides the first-order approximation. Similarly, UAT BT BBT BUA = I + UAT {2(BT B − I) + (BT B − I)2 }UA
is sufficiently close to the identity matrix. Hence,(BA)+ BBT (BA)+T ≈
Ns 
I + UAT (BT B − I)2UA
Nf

holds as the leading order. Thus, using
equation 2.13.

(4.26)≡ UAT (BT B − I)2UA , we obtain

4.6 Hebbian-Like Learning Rules. Oja’s subspace rule for PCA (Oja,
1989) is defined as
TẆPCA ∝ E uPCA (x − E[x] − WPCAuPCA )T ,(4.27)

where WPCA ∈ RNs ×Nx is a synaptic weight matrix, uPCA ≡ WPCA (x − E[x]) ∈RNs is a vector of encoders, and the dot over WPCA denotes a temporal
derivative. This rule can extract a subspace spanned by the first Ns princiT
pal components (i.e., WPCA → M PM
) in a manner equal to PCA via eigenvalue decomposition. Although Oja’s subspace rule does not have a costfunction, a gradient descent rule for PCA, based on the least mean squared

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

4.5 Derivation of Linearization Error Covariance. Equation 2.13 is de−1
rived as follows: from equation 2.4, (BH)+ = f (BA)+ + O(Ns−1 ) holds
when B is sufficiently isotropic. Thus, from equations 2.6 and 2.12, we
obtain1464

T. Isomura and T. Toyoizumi

ẆICA ∝ −
∂DKL T
WICAWICA = (I − E g(ũ)ũT )WICA ,∂WICA(4.28)

where g(ũ) ≡ −d log p0 (ũ)/dũ is a nonlinear activation function.
Acknowledgments
This work was supported by RIKEN Center for Brain Science (T.I. and T.T.),
Brain/MINDS from AMED under grant JP21dm020700 (T.T.), and JSPS
KAKENHI grant JP18H05432 (T.T.). The funders had no role in the study
design, data collection and analysis, decision to publish, or preparation of
the manuscript.
Data Availability
All relevant data are within the article. The Matlab scripts are available athttps://github.com/takuyaisomura/asymptotic_linearization.
References
Amari, S. I., Chen, T. P., & Cichocki, A. (1997). Stability analysis of learning algorithms for blind source separation. Neural Networks, 10(8), 1345–1351.
Amari, S. I., Cichocki, A., & Yang, H. H. (1996). A new learning algorithm for blind
signal separation. In D. S. Touretzky, M. C. Mozer, & M. E. Hasselmo (Eds.), Advances in neural information processing systems, 8 (pp. 757–763). Cambridge, MA:
MIT Press.
Arora, S., & Risteski, A. (2017). Provable benefits of representation learning. arXiv:1706.
04601.
Baldi, P., & Hornik, K. (1989). Neural networks and principal component analysis:
Learning from examples without local minima. Neural Networks, 2(1), 53–58.
Barron, A. R. (1993). Universal approximation bounds for superpositions of a sigmoidal function. IEEE Trans. Info. Theory, 39(3), 930–945.

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

error, is known (Xu, 1993), whose cost function is given as E[|x − E[x] −T
WPCAuPCA |2 ]. Indeed, this equals the maximum likelihood estimation ofx − E[x] using a lower-dimensional linear encoder uPCA under the assumption that the loss follows a unit gaussian distribution. The gradient descenton this cost function is equal to Oja’s subspace rule up to an additional termthat does not essentially change the behavior of the algorithm.
For BSS of uPCA , Amari’s ICA algorithm is considered (Amari et al.,1996). Using encoders ũ ≡ WICA uPCA ∈ RNs with synaptic weight matrix WICA ∈ RNs ×Ns , the ICA cost function is defined as the Kullback–Leibler divergence between the actual distribution of ũ, p(ũ), and its priorbelief p0 (ũ) ≡ i p0 (ũi ), DKL ≡ DKL [p(ũ)||p0 (ũ)] ≡ E[log p(ũ) − log p0 (ũ)].
The natural gradient of DKL yields Amari’s ICA algorithm,

Achievability of Nonlinear Blind Source Separation

1465

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Bell, A. J., & Sejnowski, T. J. (1995). An information-maximization approach to blindseparation and blind deconvolution. Neural Comput., 7(6), 1129–1159.
Bell, A. J., & Sejnowski, T. J. (1997). The “independent components” of natural scenes
are edge filters. Vision Res., 37(23), 3327–3338.
Brown, G. D., Yamada, S., & Sejnowski, T. J. (2001). Independent component analysis
at the neural cocktail party. Trends Neurosci., 24(1), 54–63.
Bussgang, J. J. (1952). Cross-correlation functions of amplitude-distorted gaussian signals(Technical Report 216). Cambridge, MA: MIT Research Laboratory of Electronics.
Calhoun, V. D., Liu, J., & Adali, T. (2009). A review of group ICA for fMRI data and
ICA for joint inference of imaging, genetic, and ERP data. NeuroImage, 45(1), S163–
S172.
Chandler, D. M., & Field, D. J. (2007). Estimates of the information content and dimensionality of natural scenes from proximity distributions. J. Opt. Soc. Am. A,
24(4), 922–941.
Chen, T., Hua, Y., & Yan, W. Y. (1998). Global convergence of Oja’s subspace algorithm for principal component extraction. IEEE Trans. Neural Netw., 9(1), 58–
67.
Cichocki, A., Zdunek, R., Phan, A. H., & Amari, S. I. (2009). Nonnegative matrix andtensor factorizations: Applications to exploratory multi-way data analysis and blindsource separation. Hoboken, NJ: Wiley.
Comon, P. (1994). Independent component analysis, a new concept? Sig. Process,
36(3), 287–314.
Comon, P., & Jutten, C. (2010). Handbook of blind source separation: Independent component analysis and applications. Orlando, FL: Academic Press.
Cybenko, G. (1989). Approximation by superpositions of a sigmoidal function. MathControl Signals Syst., 2(4), 303–314.
Dahl, G. E., Yu, D., Deng, L., & Acero, A. (2012). Context-dependent pretrained
deep neural networks for large-vocabulary speech recognition. IEEE Trans. Audio
Speech Lang. Proc., 20(1), 30–42.
Dayan, P., & Abbott, L. F. (2001). Theoretical neuroscience: Computational and mathematical modeling of neural systems. Cambridge, MA: MIT Press.
Dayan, P., Hinton, G. E., Neal, R. M., & Zemel, R. S. (1995). The Helmholtz machine.
Neural Comput., 7(5), 889–904.
DiCarlo, J. J., Zoccolan, D., & Rust, N.C. (2012). How does the brain solve visualobject recognition? Neuron, 73(3), 415–434.
Dinh, L., Krueger, D., & Bengio, Y. (2014). NICE: Non-linear independent components
estimation. arXiv:1410.8516.
Erdogan, A. T. (2007). Globally convergent deflationary instantaneous blind source
separation algorithm for digital communication signals. IEEE Trans. Signal Process., 55(5), 2182–2192.
Erdogan, A. T. (2009). On the convergence of ICA algorithms with symmetric orthogonalization. IEEE Trans. Signal Process., 57(6), 2209–2221.
Földiák, P. (1990). Forming sparse representations by local anti-Hebbian learning.
Biol. Cybern., 64(2), 165–170.
Friston, K. (2008). Hierarchical models in the brain. PLOS Comput. Biol., 4, e1000211.
Friston, K., Trujillo-Barreto, N., & Daunizeau, J. (2008). DEM: A variational treatmentof dynamic systems. NeuroImage, 41(3), 849–885.

1466

T. Isomura and T. Toyoizumi

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Gerstner, W., & Kistler, W. M. (2002). Spiking neuron models: Single neurons, populations, plasticity. Cambridge: Cambridge University Press.
Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep learning. (Cambridge, MA:
MIT Press).
Griffiths, D. J. (2005). Introduction to quantum mechanics (2nd ed.). Upper Saddle River,
NJ: Pearson Prentice Hall.
Hebb, D. O. (1949). The organization of behavior: A neuropsychological theory. New York:
Wiley.
Hinton, G. E., & Salakhutdinov, R. R. (2006). Reducing the dimensionality of datawith neural networks. Science, 313(5786), 504–507.
Hinton, G. E., Srivastava, N., Krizhevsky, A., Sutskever, I., & Salakhutdinov, R. R.
(2012). Improving neural networks by preventing co-adaptation of feature detectors.
arXiv:1207.0580.
Hornik, K., Stinchcombe, M., & White, H. (1989). Multilayer feedforward networks
are universal approximators. Neural Netw., 2(5), 359–366.
Hyvärinen, A., & Morioka, H. (2016). Unsupervised feature extraction by timecontrastive learning and nonlinear ICA. In D. Lee, M. Sugiyama, U. Luxburg,
I. Guyon, & R. Garnett (Eds.), Advances in neural information processing systems, 29(pp. 3765–3773). Red Hook, NY: Curran.
Hyvärinen, A. J., & Morioka, H. (2017). Nonlinear ICA of temporally dependent stationary sources. In Proceedings of the Conference on Machine Learn
Research.
Hyvärinen, A., & Oja, E. (1997). A fast fixed-point algorithm for independent component analysis. Neural Comput., 9(7), 1483–1492.
Hyvärinen, A., & Pajunen, P. (1999). Nonlinear independent component analysis:
Existence and uniqueness results. Neural Netw., 12(3), 429–439.
Isomura, T., & Friston, K. (2018). In vitro neural networks minimise variational freeenergy. Sci. Rep., 8, 16926.
Isomura, T., & Friston, K. (2020). Reverse engineering neural networks to characterize their cost functions. Neural Comput., 32(11), 2085–2121.
Isomura, T., Kotani, K., & Jimbo, Y. (2015). Cultured cortical neurons can perform
blind source separation according to the free-energy principle. PLOS Comput.
Biol., 11(12), e1004643.
Isomura, T., & Toyoizumi, T. (2016). A local learning rule for independent component
analysis. Sci. Rep., 6, 28073.
Isomura, T., & Toyoizumi, T. (2018). Error-gated Hebbian rule: A local learning rule
for principal and independent component analysis. Sci. Rep., 8, 1835.
Isomura, T., & Toyoizumi, T. (2019). Multi-context blind source separation by errorgated Hebbian rule. Sci. Rep., 9, 7127.
Isomura, T., & Toyoizumi, T. (2021). Dimensionality reduction to maximize prediction generalization capability. Nat. Mach. Intell., in press.
Jolliffe, I. T. (2002). Principal component analysis (2nd ed.). Berlin: Springer.
Jutten, C., & Karhunen, J. (2004). Advances in blind source separation (BSS) and independent component analysis (ICA) for nonlinear mixtures. Int. J. Neural. Syst.,
14(5), 267–292.
Kandel, E. R., Schwartz, J. H., Jessell, T., M., Siegelbaum, S. A., & Hudspeth, A. J.
(2013). Principles of neural science (5th ed.). New York: McGraw-Hill.

Achievability of Nonlinear Blind Source Separation

1467

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Karhunen, J. (2001). Nonlinear independent component analysis. In S. Roberts &R. Everson (Eds.), Independent component analysis: Principles and practice (pp. 113–
134). Cambridge: Cambridge University Press.
Kawaguchi, K. (2016). Deep learning without poor local minima. In D. Lee, M.
Sugiyama, U. Luxburg, I. Guyon, & R. Garnett (Eds.), Advances in neural information processing systems, 29. Red Hook, NY: Curran.
Khemakhem, I., Kingma, D., Monti, R., & Hyvarinen, A. (2020). Variational autoencoders and nonlinear ICA: A unifying framework. In Proceedings of the International Conference on Artificial Intelligence and Statistics 2207–2217. ML Research
Press.
Kingma, D. P., & Welling, M. (2013). Auto-encoding variational Bayes. arXiv:1312.6114.
Lappalainen, H., & Honkela, A. (2000). Bayesian non-linear independent component
analysis by multi-layer perceptrons. In M. Girolami, (Ed.), Advances in independent
component analysis (pp. 93–121). London: Springer.
Leugering, J., & Pipa, G. (2018). A unifying framework of synaptic and intrinsic plasticity in neural populations. Neural Comput., 30(4), 945–986.
Linsker, R. (1997). A local learning rule that enables information maximization for
arbitrary input distributions. Neural Comput., 9(8), 1661–1665.
Lu, H., & Kawaguchi, K. (2017). Depth creates no bad local minima. arXiv:1702.08580.
Malenka, R. C., & Bear, M. F. (2004). LTP and LTD: An embarrassment of riches.
Neuron, 44(1), 5–21.
Marchenko, V. A., & Pastur, L. A. (1967). Distribution of eigenvalues for some setsof random matrices. Mat. Sbornik, 114, 507–536.
Mika, D., Budzik, G., & Józwik, J. (2020). Single channel source separation with ICAbased time-frequency decomposition. Sensors, 20(7), 2019.
Nguyen, Q., & Hein, M. (2017). The loss surface of deep and wide neural networks.
arXiv:1704.08045.
Oja, E. (1982). Simplified neuron model as a principal component analyzer. J. Math.
Biol., 15(3), 267–273.
Oja, E. (1989). Neural networks, principal components, and subspaces. Int. J. Neural.
Syst., 1(10), 61–68.
Oja, E., & Yuan, Z. (2006). The FastICA algorithm revisited: Convergence analysis.
IEEE Trans. Neural Netw., 17(6), 1370–1381.
Papadias, C. B. (2000). Globally convergent blind source separation based on a multiuser kurtosis maximization criterion. IEEE Trans. Signal Process., 48(12), 3508–
3519.
Pearson, K. (1901). On lines and planes of closest fit to systems of points in space.
Philos. Mag., 2(11), 559–572.
Pehlevan, C., Mohan, S., & Chklovskii, D. B. (2017). Blind nonnegative source separation using biological neural networks. Neural Comput., 29(11), 2925–2954.
Rahimi, A., & Recht, B. (2008a). Uniform approximation of functions with randombases. In Proceedings of the 46th Annual Allerton Conference on Communication, Control, and Computing (pp. 555–561). Red Hook, NY: Curran.
Rahimi, A., & Recht, B. (2008b). Weighted sums of random kitchen sinks: Replacing
minimization with randomization in learning. In D. Koller, D. Schuurmans, Y.
Bengio, & L. Bottou (Eds.), Advances in neural information processing systems, 21(pp. 1313–1320). Cambridge, MA: MIT Press.

1468

T. Isomura and T. Toyoizumi

Received July 26, 2020; accepted December 23, 2020.

Downloaded from http://direct.mit.edu/neco/article-pdf/33/6/1433/1916364/neco_a_01378.pdf by guest on 02 June 2021

Sanger, T. D. (1989). Optimal unsupervised learning in a single-layer linear feedforward neural network. Neural Netw., 2(6), 459–473.
Toyoizumi, T., & Abbott, L. F. (2011). Beyond the edge of chaos: Amplification andtemporal integration by recurrent networks in the chaotic regime. Phys. Rev. E.,
84(5), 051908.
van der Lee, T., Exarchakos, G., & de Groot, S. H. (2019). In-network Hebbian plasticity for wireless sensor networks. In Proceedings of the International Conference on
Internet and Distributed Computing Systems (pp. 79–88). Cham: Springer.
Wan, L., Zeiler, M., Zhang, S., LeCun, Y., & Fergus, R. (2013). Regularization of neuralnetworks using DropConnect. In Proceedings of Machine Learning Research (28)3,
1058–1066.
Wentzell, P. D., Andrews, D. T., Hamilton, D. C., Faber, K., & Kowalski, B. R. (1997).
Maximum likelihood principal component analysis. J. Chemom., 11(4), 339–366.
Xu, L. (1993). Least mean square error reconstruction principle for self-organizing
neural–nets. Neural Netw., 6(5), 627–648.