Tang, Huanwen

Tang, Yiyuan

Wang, Huili

Zhong, Mingjun

An expectation-maximization (EM) algorithm for independent component analysis in the presence of gaussian noise is presented. The estimation of the conditional moments of the source posterior can be accomplished by maximum a posteriori estimation. The approximate conditional moments enable the development of an EM algorithm for inferring the most probable sources and learning the parameter in noisy independent component analysis. Simulation results show that the proposed method can perform blind source separation of sub-Gaussian mixtures and super-Gaussian mixtures

Edinburgh Research Explorer

     Edinburgh Research Explorer                                      An EM Algorithm for Independent Component Analysis in thePresence of Gaussian NoiseCitation for published version:Zhong, M, Tang, H, Wang, H & Tang, Y 2004, 'An EM Algorithm for Independent Component Analysis in thePresence of Gaussian Noise' Neural Information Processing - Letters and Reviews, pp. 11-17.Link:Link to publication record in Edinburgh Research ExplorerDocument Version:Publisher's PDF, also known as Version of recordPublished In:Neural Information Processing - Letters and ReviewsGeneral rightsCopyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)and / or other copyright owners and it is a condition of accessing these publications that users recognise andabide by the legal requirements associated with these rights.Take down policyThe University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorercontent complies with UK legislation. If you believe that the public display of this file breaches copyright pleasecontact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately andinvestigate your claim.Download date: 05. Apr. 2019’An EM Algorithm for Independent Component Analysis in the Presenceof Gaussian NoiseMingjun Zhong and Huanwen TangInstitute of Computational Biology and Bioinformatics, Dalian University of TechnologyDalian 116023, People’s Republic of ChinaE-mail:sunxl zhong@yahoo.comHuili WangInstitute of Neuroinformatics, and Department of Foreign LanguagesDalian University of Technology, Dalian 116023, People’s Republic of ChinaE-mail:graduate@dlut.edu.cnYiyuan TangInstitute of Neuroinformatics, Dalian University of TechnologyDalian 116023, People’s Republic of ChinaLaboratory of Visual Information Processing and Key Lab for Mental HealthThe Chinese Academy of Sciences, Beijing 100101, People’s Republic of ChinaE-mail:yy2100@163.net(Submitted on November14,2003)Abstract—An expectation-maximization (EM) algorithm for independent component analysis in thepresence of gaussian noise is presented. The estimation of the conditional moments of the sourceposterior can be accomplished by maximum a posteriori estimation. The approximate conditionalmoments enable the development of an EM algorithm for inferring the most probable sources andlearning the parameter in noisy independent component analysis. Simulation results show that theproposed method can perform blind source separation of sub-Gaussian mixtures and super-Gaussianmixtures.Keywords—Independent component analysis, blind source separation, EM algorithm, maximum aposteriori1. IntroductionIndependent component analysis (ICA) [1, 2, 3, 4, 5] is a statistical technique whose goal is to recoverindependent components given only observations that are unknown linear mixtures of the unobserved independentsource signals. One may formulate ICA as the estimation of the following linear generative model for the data:x = As+ ε, (1)where x = (x1, · · · , xN )T is the vector of observed random variables, s = (s1, · · · , sM )T is the vector of theindependent components, A ∈ RN×M is the unknown constant mixing matrix, and the vector ε = (ε1, · · · , εN )Tis the noise modelled as Gaussian with zero mean and covariance matrix Σ. In the expectation-maximization (EM)algorithm [6] proposed in this paper, we assume the dimension of s equals the dimension of x, i.e., N = M , andthe noise covariance matrix Σ is known.It is shown in this paper that the conditional moments of the source posterior distribution can be estimatedby maximum a posteriori (MAP) estimation. Given these approximate conditional moments, an EM algorithm isLETTER Neural Information Processing-Letters and Reviews Vol. 2, No. 1, January 200411presented for learning the parameter A in equation (1). To illustrate the EM algorithm proposed in this paper, sim-ulation results are presented for blind source separation (BSS) [2, 7] of super-Gaussian mixture and sub-Gaussianmixture signals.2. The Approximate Conditional MomentsThe approximate conditional moments of the posterior can be derived by modelling the posterior as gaus-sian. Assuming that we have observed T data samples x = x(1), · · · , x(T )which are generated according to model(1), the approximate conditional moments for the gaussian posterior distribution, based on the Laplace estimation[8, 9], are given as (see appendix A):E{s|x(t)} = sˆ = argmaxsL(s) (2)E{ssT |x(t)} = H(sˆ)−1 + sˆsˆT , (3)where L(s) = log{p(s|x(t), A)}, where p denotes the probability distribution function (p.d.f.) which is usedthroughout this paper, is the log-posterior, and the Hessian of the approximate log-posterior computed at the MAPvalue sˆ is denoted as H(sˆ) = −∇s∇slog{p(s|x(t), A)} = ATΣ−1A − ∇s∇slog{p(sˆ)} where ∇s denotes thegradient with respect to s used throughout this paper.The estimation of the independent components can be accomplished by MAP estimation. According to thegenerative model (1), the log-posterior, i.e., L(s), is given as follows:L(s(1), · · · , s(T )) =T∑t=1{−12(x(t)−As(t))TΣ−1(x(t)−As(t)) + ϕ(s(t))}+ C, (4)where ϕ(s(t)) = log{p(s(t))} and C is a constant irrelevant to s(t). To estimate the independent components,taking the gradient of the log-posterior with respect to s(t) gives the following gradient learning rule:∇sL(s(t)) = ATΣ−1(x(t)−As(t)) + ϕ′(s(t)), (5)where ϕ′(s(t)) = ∇sϕ(s(t)) is called activation function in some literatures [7, 10]. Note that the log-posterior inequation (4) is essentially the same as the likelihood proposed by Hyva¨rinen [3]. It is indicated in the followingsection that the approximate conditional moments in equations (2) and (3) enable the development of an EMalgorithm for learning the parameter A in model (1).3. An EM Algorithm for Parameter EstimationTo derive a learning algorithm for inferring A in model (1), it is required to maximize the following datalikelihood given the model:p(x|A) =∫p(x|s,A)p(s)ds. (6)Rather than using some approximations to this intractable integral as in [8, 11], it is desirous to employ the EMframework for estimation and inference for this form of linear model, as this is a most natural method for maximiz-ing the data likelihood. Given the approximate conditional moments, the standard form of M-step for A emerges(see appendix B):Anew = (T∑t=1x(t)E{s|x(t)}T )(T∑t=1E{ssT |x(t)})−1. (7)Inserting equations (2) and (3) into (7) and setting Λ(t) = −∇s∇sϕ(sˆ(t)), we have:Anew = (T∑t=1x(t)sˆ(t)T )[T∑t=1{(ATΣ−1A+ Λ(t))−1 + sˆ(t)sˆ(t)T }]−1. (8)An EM Algorithm for Independent Component Analysis in the Presence M. Zhong, et. al.120 5000 10000−1−0.500.510 5000 10000−1−0.500.510 5000 10000−1−0.500.510 5000 10000−0.0500.050 5000 10000−0.1−0.0500.050.10 5000 10000−0.1−0.0500.050.10 5000 10000−0.1−0.0500.050.10 5000 10000−0.1−0.0500.050.10 5000 10000−0.0500.050.10.15(a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 1. (a-c)The original sources.(d-f)The noisy observations.(g-i)The inferred sources using the EM algorithmoutlined.Since sˆ(t) must be inferred in each M-step, a simple alternating variable method, which has already been usedin similar estimation tasks [3], should be derived for learning and inference. The method is based on first theoptimization of the objective function with respect to s(t) for fixed A, then the optimization with respect to Afor fixed s(t), and so on. The optimization with respect to s(t) for fixed A is accomplished by equation (5), andthe optimization with respect to A for fixed s(t) is accomplished by equation (8). Thus, the EM procedure of thefollowing form is obtained:(i)Taking some initial value for A0. Set s0 = (A0)+x where (A0)+ denotes the Moore-Penrose pseudo-inverse ofA0, and let k = 1. Normalize each column of A0 to have unit norm.(ii)Set sk(t) = sk−1(t) + η∇sL(sk−1(t)), using Ak−1 as the estimate of A.(iii)Compute Ak by equation (8), using Ak−1 as the estimate of A and substituting sk(t) for sˆ(t). Normalize eachcolumn of Ak to have unit norm.(iv)Set k = k + 1, and go back to step (ii) if not converged.In this algorithm η refers to some learning rate.4. SimulationsWe applied the EM algorithm of our paper for blind separation of independent components from noisymixtures. According to the dichotomy proposed by Hyva¨rinen [3], to illustrate the EM algorithm proposed in thispaper, it would be enough to use just two different densities in equation (4), one corresponds to a distribution thatis super-Gaussian (i.e., has a positive kurtosis) and another one corresponds to a distribution that is sub-Gaussian(i.e., has a negative kurtosis).For the inference and estimation in the case of super-Gaussian independent components, three speech sig-nals∗ are used. Each speech signal has a positive kurtosis. The activation function ϕ′(s) = −tanh(s) [2, 7, 10] isemployed in this simulation. The three speech signals were mixed by a randomly chosen matrix A, and noisy datawith zero mean and covariance matrix Σ = 0.01I were added to the mixtures. The EM algorithm proposed wasused for learning A and inferring the independent components s. The proposed algorithm procedure was randomly∗Where the three speech signals were obtained from http://www.cnl.salk.edu/˜ tewon/Over.Neural Information Processing-Letters and Reviews Vol. 2, No. 1, January 20041302468101214sourcesKL distances1 2 3 Figure 2. The KL distances between the super-Gaussian sources and the corresponding recovered sources. Theleft and right bars in each group represent the absolute values of the logarithms of the KL distances between thesources and the recovered sources obtained by the EM algorithm and the method reported in [3], respectively.initialized and converged in about 20 parameter updates. Figure 1 shows 10,000 samples of the original sources(a-c), the noisy observations (d-f) and the inferred sources (g-i). To assess the accuracy of the recovered sources,we compute the Kullback-Leibler (KL) distance [12] between the source density and the recovered source density.Since both densities are assumed to be super-Gaussian, the KL distance can be estimated via [13]:Km =∫p(sm)logp(sm)p(sˆm)dsm ≈T∑t=1logp(sm(t))p(sˆm(t)), (9)where p(s) ∝ cosh−1(s), and m = 1, · · · ,M . For the comparative purpose, we also computed the KL distancesbetween the source densities and the recovered source densities obtained by the method reported in [3] under thesame requirements. Figure 2 shows the absolute values of the logarithms of the KL distances. Note that since thelogarithms of the KL distances are negative values, the KL distances decrease respectively over those obtained bythe method reported in [3].In the case of sub-Gaussian independent components, the cubic nonlinear function, i.e., ϕ′(s) = −s3,has been a favorite choice as the activation function in equation (5) [7]. Binary and uniform distributions areused as the original distributions of the independent components. Each original signal shown in Figure 3 has anegative kurtosis. The “sign” nonlinearity which was originally used in [3] is employed for estimating A. Hence,the M-step in equation (8) has a new form of Anew = (∑Tt=1 x(t)sign(sˆ(t)T ))[∑Tt=1{(ATΣ−1A + Λ(t))−1 +sign(sˆ(t))sign(sˆ(t)T )}]−1, where the “sign” operator is applied separately on each component of its argumentin this simulation. The three signals were mixed by a randomly chosen matrix A, and noisy data with zero meanand covariance matrix Σ = 0.01I were added to the mixtures. Figure 3 shows the original sources (a-c), the noisyobservations (d-f), the inferred sources without using the “sign” nonlinearity (g-i) and the inferred sources usingthe “sign” nonlinearity (j-l). Figure 4 shows the absolute values of the logarithms of the KL distances. Note thatsince the logarithms of the KL distances are negative values, the KL distances decrease respectively over thoseobtained by the method reported in [3].These two simulations serve to show that the EM algorithm proposed in this paper is able to estimate theAn EM Algorithm for Independent Component Analysis in the Presence M. Zhong, et. al.140 50 100−2020 50 100−2020 50 100−2020 50 100−2020 50 100−2020 50 100−2020 50 100−0.200.20 50 100−0.200.20 50 100−0.200.20 50 100−0.200.20 50 100−0.200.20 50 100−0.200.2(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 3. (a-c)The original sources.(d-f)The noisy observations.(g-i)The inferred sources without using the “sign”nonlinearity. (j-l)The inferred sources using the “sign” nonlinearity.024681012sourcesKL distances1 2 3 Figure 4. The KL distances between the sub-Gaussian sources and the corresponding recovered sources. Theleft and right bars in each group represent the absolute values of the logarithms of the KL distances between thesources and the recovered sources obtained by the EM algorithm and the method reported in [3], respectively.Neural Information Processing-Letters and Reviews Vol. 2, No. 1, January 200415mixing matrix and infer the most probable independent components in noisy ICA.5. ConclusionsWe have proposed an EM algorithm for learning the parameter A and inferring the most probable indepen-dent components in noisy ICA. The approximate conditional moments of the posterior distribution estimated byMAP estimation are used to develop an EM algorithm for learning the parameter and inferring the most probableindependent components.Some work on ICA in the presence of Gaussian noise has been done in some literatures [3, 10]. One methodmore related to ours was introduced in [13], in which the EM algorithm was used for independent factor analysis.The proposed method is simple and computationally less demanding than the one proposed in [13]. In particular,the complexity of the method proposed in [13] is exponential as a function of the number of the independent com-ponents, and the algorithm in this paper has polynomial complexity. Future work will estimate the noise covariancematrix in this EM framework.Appendix A. Derivation of the Approximate Conditional MomentsBased on the Laplace estimation [8], the posterior distribution of the sources can be modelled as a gaussiandistribution:p(s|x(t)) ≈ (2pi)−M2 |H(sˆ)| 12 exp{−12(s− sˆ)TH(sˆ)(s− sˆ)}, (10)where sˆ is the MAP value and H(sˆ) is the Hessian of the approximate log-posterior computed at the MAP value.Thus, the approximate conditional moments in equations (2) and (3) can be derived by this gaussian posteriordistribution.Appendix B. Derivation of the M-stepTo derive the M-step for the parameter A, it is required to maximize the expected value of the completedata likelihood of the following form given the observed data x and the current model [6]:Q(Anew|Aold) =T∑t=1∫p(s|x(t), Aold)log{p(x(t), s|Anew)}ds, (11)where Aold is the parameter obtained in the previous iteration. Setting the gradient of Q(Anew|Aold) with respectto Anew to zero gives the new value of the parameter of the M-step in equation (7).Acknowledgments: The work was supported by NSFC(30170321,90103033), MOST(2001CCA00700) and MOE(KP0302).References[1] P.Comon, “Independent Component Analysis - A New Concept?”Signal Processing, Vol. 36, No. 3, pp. 287-314, April 1994.[2] A.J.Bell,T.J.Sejnowski, “An Information-maximization Approach to Blind Separation and Blind Deconvolu-tion,”Neural Computation, Vol. 7, No. 6, pp. 1129-1159, Nov. 1995.[3] A.Hyva¨rinen, “Independent Component Analysis in the Presence of Gaussian Noise by Maximizing JointLikelihood,”Neurocomputing, Vol. 22, No. 1-3, pp. 49-67, Nov. 1998.An EM Algorithm for Independent Component Analysis in the Presence M. Zhong, et. al.16[4] T-W.Lee,M.Girolami,T.J.Sejnowski, “Independent Component Analysis Using an Extended Infomax Algo-rithm for Mixed Subgaussian and Supergaussian Sources,”Neural Computation, Vol. 11, No. 2, pp. 417-441,Feb. 1999.[5] A.Hyva¨rinen, E.Oja, “Independent Component Analysis: Algorithms and Applications,”Neural Networks,Vol. 13, No. 4-5, pp. 411-430, June 2000.[6] A.P.Dempster,N.M.Laird,D.B.Rubin, “Maximum Likelihood from Incomplete Data via the EM Algo-rithm(with discussion),”Journal of the Royal Statistical Society B, Vol. 39, No. 1, pp. 1-38, Nov. 1977.[7] A.Cichocki,S.I.Amari, Adaptive Blind Signal and Image Processing: Learning Algorithms and Applications,Wiley, New York, 2002.[8] M.S.Lewicki,T.J.Sejnowski, “Learning Overcomplete Representations,”Neural Computation, Vol. 12, No. 2,pp. 337-365, Feb. 2000.[9] M.Girolami, “A Variational Method for Learning Sparse and Overcomplete Representations,”Neural Compu-tation, Vol. 13, No. 11, pp. 2517-2532, Nov. 2001.[10] A.Cichocki,S.C.Douglas,S.I.Amari, “Robust Technique for Independent Component Analysis with NoisyData,”Neurocomputing, Vol. 22, No. 1-3, pp. 113-129, Nov. 1998.[11] R.M.Neal, Bayesian Learning for Neural Networks, Lecture Notes in Statistics, Vol.118, Springer-Verlag,New York, 1996.[12] T.M.Cover,J.A.Thomas, Elements of Information Theory, Wiley, New York, 1991.[13] H.Attias, “Independent Factor Analysis,”Neural Computation, Vol. 11, No. 4, pp.803-951, May 1999.Mingjun Zhong studied mathematics and obtained his Master’s degree at the Dalian University of Technology,P.R.China. He is currently studying for a Ph.D. degree at the Institute of Computational Biology and Bioinfor-matics of the Dalian University of Technology, P.R.China. His research interests include independent componentanalysis, sparse coding, bioinformatics and neuroinformatics.Huanwen Tang graduated from the Dalian University of Technology in 1963, and currently a professor of math-ematics and management engineering at the Dalian University of Technology. His research interest includes com-putational models and algorithm of human cognition and neural information coding, bioinformatics and neuroin-formatics. (Home page: http://brain.dlut.edu.cn)Huili Wang graduated from the Shanghai International Studies University in 1989, and currently an assistantprofessor of school of foreign languages at the Dalian University of Technology. Her research interest includespsychal linguistics, neurolinguistics and neuroinformatics.Yiyuan Tang graduated from the Jinlin University in 1987, and currently a professor of neuroinformatics andneuroscience at the Dalian University of Technology. His research interest includes neuroimaging, cognition andemotion interaction, neuroinformatics. (Home page: http://brain.dlut.edu.cn)Neural Information Processing-Letters and Reviews Vol. 2, No. 1, January 200417

An EM Algorithm for Independent Component Analysis in the Presence of Gaussian Noise

https://www.pure.ed.ac.uk/ws/files/24656967/V02N01P2_11_17.pdf

An EM Algorithm for Independent Component Analysis in the Presence of Gaussian Noise

Abstract

Similar works

Full text

Available Versions

Edinburgh Research Explorer