Barber, D

Chiappa, S

English

We discuss a method to extract independent dynamical systems underlying a single or multiple channels of observation. In particular, we search for one-dimensional subsignals to aid the interpretability of the decomposition. The method uses an approximate Bayesian analysis to determine automatically the number and appropriate complexity of the underlying dynamics, with a preference for the simplest solution. We apply this method to unfiltered EEG signals to discover low-complexity sources with preferential spectral properties, demonstrating improved interpretability of the extracted sources over related methods. © 2007 IEEE

UCL Discovery

SIGNAL PROCESSING LETTERS 1Bayesian Factorial Linear Gaussian State-SpaceModels for Biosignal DecompositionSilvia Chiappa and David BarberAbstract—We discuss a method to extract independent dynam-ical systems underlying a single or multiple channels of observa-tion. In particular, we search for one dimensional subsignals toaid the interpretability of the decomposition. The method usesan approximate Bayesian analysis to determine automatically thenumber and appropriate complexity of the underlying dynamics,with a preference for the simplest solution. We apply this methodto unﬁltered EEG signals to discover low complexity sourceswith preferential spectral properties, demonstrating improvedinterpretability of the extracted sources over related methods.I. INTRODUCTIONDecomposing a multivariate time-series vnt , t = 1,...,T,n = 1,...,V into a set of C simpler subsignals (sources) isa central goal in signal processing and is of particular interestin the analysis of biomedical signals. The goal of this paperis to introduce a model which can automatically determinethe number of sources underlying the observations and inwhich we can bias the sources to be in certain frequencyranges. Furthermore, we are interested in taking into accountthe temporal structure of the time-series which can help inobtaining a good decomposition, especially when C > V .More speciﬁcally, our criterion for the decomposition is thatindependent dynamical systems generate the sources which,under linear noisy mixing, give rise to the observations. Forany two scalar sources sit and sjt and all times t, we seek amodel of statistically independent dynamics p(si1:T,sj1:T) =p(si1:T)p(sj1:T). Furthermore, the aim is to ﬁnd a matrix Wthat relates the sources st = vert(s1t,...,sCt ) to observationsvt = vert(v1t,...,vVt ) through noisy mixing1. This is a formof Independent Components Analysis (ICA) [1] although itdiffers from the more standard assumption of independenceat each time step, that is p(si1:T,sj1:T) =QTt=1 p(sit)p(sjt).We consider a Linear Gaussian State-Space Model (LGSSM),which is a powerful, yet interpretable and tractable, model.We constrain the LGSSM in order that independent dynamicalprocesses can be identiﬁed and furthermore that scalar sourcescan be extracted from the signal. To determine the correctnumber of underlying processes and bias the solution towardsa certain dynamics, we use a Variational Bayesian analysiswhich deﬁnes a prior distribution over the model parameters.There are several existing decomposition methods whichencode constraints such as desired frequencies of the inde-pendent sources (see for example [2], [3]). However, thesemethods do not automatically determine the correct numberof underlying sources nor do they consider the dynamics ofthe signal in the model structure. A closely related technique1vert(a,b,c) is the matrix formed by vertically stacking a, b and c.to ours is (Non) Linear Dynamical Factor Analysis (NDFA)[4], [5]. Whilst being an attractive and powerful method,standard NDFA places no constraint that the observations areformed from mixing independent scalar sources, which makesinterpretation of the resulting sources difﬁcult. Furthermore,NDFA does not directly force the sources to contain particularfrequencies but rather attempts to bias the discovered sourcesby careful initialization [5]. In addition, NDFA uses nonlinearstate dynamics (and mixing), which hampers inference andmakes the incorporation of known constraints more complex.Inference in the Variational Bayesian LGSSM has previ-ously been achieved using Belief Propagation, and differsfrom inference in the Kalman ﬁltering/smoothing literature,for which highly efﬁcient and stabilized procedures exist. Acentral contribution of this paper is to show how inference canbe performed using the standard Kalman ﬁltering/smoothingrecursions by augmenting the original model.II. FACTORIAL LINEAR GAUSSIAN STATE-SPACE MODELSIn LGSSMs [6], the hidden state vectors h1:T and the visibleobservations v1:T are linearly related by:ht = Aht−1 + ηht , h1 ∼ N (µ,Σ), ηht ∼ N (0H,ΣH)vt = Bht + ηvt , ηvt ∼ N(0V ,ΣV ),where N denotes a Gaussian distribution. The notation 0Dstands for a D × 1 zero vector. Probabilistically:p(v1:T,h1:T) = p(v1|h1)p(h1)T Yt=2p(vt|ht)p(ht|ht−1),with p(vt|ht)=N(Bht,ΣV ) and p(ht|ht−1)=N(Aht−1,ΣH).To make independent dynamical subsystems we use blockdiagonal transition and state noise matrices A, ΣH and Σ,where each block c has dimension Hc. A one dimensionalsource sct for each independent dynamical subsystem isformed from sct = 1Tchct, where 1c is a unit vector and hct isthe state of the dynamical system c. Combining the sources,we can write st = Pht, where P = diag(1T1,...,1TC),ht = vert(h1t,...,hCt ). The resulting emission matrix isconstrained to be of the form B = WP, where W is theV × C mixing matrix and P is a C × H projection, withH =Pc Hc. Such a constrained form for B is required toprovide interpretable scalar sources.BAYESIAN FACTORIAL LINEAR GAUSSIAN STATE-SPACEMODELSIn our Bayesian treatment of learning we deﬁne the priorsp(A|α) and p(W|β), where α and β are hyperparameters. WeSIGNAL PROCESSING LETTERS 2do not deﬁne any prior for ΣH, ΣV , µ and Σ, which willformally be considered as hyperparameters2. The total set ofhyperparameters is Θ = {α,β,ΣH,ΣV ,µ,Σ}. Therefore:p(v1:T|Θ) =ZA,Wp(v1:T|A,W,Θ)p(A|α)p(W|β)dAdW . (1)Here we take the ML-II (‘evidence’) framework, which in-volves maximizing p(v1:T|Θ) with respect to Θ [4], [7].Ideally, the number of sources effectively contributing to theobserved signal should be small. This suggests the prior:p(W|β) =C Yj=1￿βj2π￿V/2e−βj2PVi=1 W2ij .We can bias A to be close to a desired transition ˆ A (possiblyzero) by using:p(Ac|αc) =￿αc2π￿H2c/2e−αc2PHci,j=1(Acij− ˆ Acij)2for each component c, so that p(A|α) =Qc p(Ac|αc)3.Variational BayesOptimizing Eq. (1) with respect to Θ is difﬁcult due to theintractability of the integrals. Instead we consider the lowerbound [4], [7]:L = logp(v1:T|Θ) ≥Hq(A,W,h1:T) (2)+  logp(v1:T,h1:T,A,W) q(A,W,h1:T) ,where we dropped the explicit dependence on Θ on the rhs4.The notation Hd(x) signiﬁes the entropy of the distributiond(x), and    d(x) denotes the expectation operator. For certainvariational distributions q, we hope to achieve a tractablebound, which we may then optimize with respect to q andΘ. The key approximation in Variational Bayes (VB) isq(A,W|h1:T) ≡ q(A,W). Since A and W separate in therhs of Eq. (2), optimally q(A,W) = q(A)q(W), hence:L ≥ − D(q(A),p(A)) − D(q(W),p(W)) + Hq(h1:T)+  logp(v1:T,h1:T|A,W) q(h1:T)q(A)q(W) ≡ F .D(q(x),p(x)) is the KL divergence  logq(x)/p(x) q(x). TheVB procedure iteratively performs co-ordinate wise ascent ofF with respect to q(W),q(A),q(h1:T) and Θ.Determining q(W)By examining F, optimally, q(W) is a Gaussian. The co-variance [ΣW]ij,kl ≡￿￿Wij −  Wij ￿￿Wkl −  Wkl ￿￿(aver-ages wrt q(W)) is given by the inverse of:￿Σ−1W￿ij,kl =￿Σ−1V￿ikXtD˜ hjt˜ hltEq(ht)+ βjδikδjl ,2A Bayesian treatment of ΣH, ΣV , µ and Σ is straightforward usingconjugate priors (see [7], [8]) but is not taken here for space restrictionsand since we have little preference for constraining these parameters.3For dimensional reasons, we can also assume a Gaussian prior on thecolumns of W with exponent − 12βjWTj Σ−1V Wj. This simpliﬁes the statisticsof q(W) and Eq. (4). The same holds for A. This is also convenient whenwe assume a prior for ΣH and ΣV , since it ensures conjugacy [7], [8].4Strictly we should write throughout q(·|v1:T). We omit the dependenceon v1:T for notational convenience.where ˜ ht = Pht and δij is the Kronecker delta function. Themean is given by: Wij  =Xk,l,n,t[ΣW]ij,kl￿Σ−1V￿k,nD˜ hltEq(ht)vnt .Determining q(A)Optimally we have a factorized distribution q(A) = Qc q(Ac), where q(Ac) is Gaussian with inverse covariancegiven by (dropping the dependence on c):￿Σ−1A￿ij,kl =￿Σ−1H￿ikT Xt=2Dhjt−1hlt−1Eq(ht−1)+ αδikδjl .The mean is: Aij =Xk,l[ΣA]ij,kl α ˆ Akl+Xn￿Σ−1H￿knT Xt=2￿hlt−1hnt￿q(ht−1:t)!.Inference on q(h1:T)Optimally q(h1:T) is Gaussian since its log is quadratic inh1:T, being namely5:−12T Xt=1D(vt − WPht)TΣ−1V (vt − WPht)Eq(W)(3)−12T Xt=2D(ht − Aht−1)T Σ−1H (ht − Aht−1)Eq(A).We can carry out the averages over A and W since q(A)and q(W) are Gaussian and the above is quadratic in theparameters A and W. In order to compute the requiredstatistics  ht q(ht) and￿ht−1hTt￿q(ht−1:t), our aim is to rep-resent Eq. (3) as the log ˜ q(h1:T|˜ v1:T) of a LGSSM with somesuitable parameters. To do that we use a mean + ﬂuctuationdecomposition:D(vt − Bht)TΣ−1V (vt − Bht)Eq(W)= (vt −  B ht)TΣ−1V (vt −  B ht) + hTtPTSWPht ,where  B  ≡  W P and the ﬂuctuation is by determined by:[SW]jl =V Xi,k=1[ΣW]ij,kl￿Σ−1V￿ik , j,l ∈ 1,...,C. (4)Similarly:D(ht − Aht−1)TΣ−1H (ht − Aht−1)Eq(A)= (ht −  A ht−1)TΣ−1H (ht −  A ht−1) + hTt−1SAht−1 ,[SA]jl =H Xi,k=1[ΣA]ij,kl￿Σ−1H￿ik , j,l ∈ 1,...,H.To represent Eq. (3) as a LGSSM ˜ q(h1:T|˜ v1:T), we augmentvt and B as:˜ vt = vert(vt,0H,0C), ˜ B = vert( B ,UA,UWP),5For simplicity, we ignore the contribution from h1 and a constant term.SIGNAL PROCESSING LETTERS 30 50 100 150 200 250 300(a)0 50 100 150 200 250 300(b)0 50 100 150 200 250 300(c)0 50 100 150 200 250 300(d)Fig. 1. (a) Original sources st. (b) Observations resulting from mixing theoriginal sources, vt = Wst + ηvt , ηvt ∼ N(0V ,I). (c) Recovered sourcesusing the Bayesian FLGSSM. (d) Sources found with MAP FLGSSM. Theretained sources have been rescaled to aid visualization.where UA is the Cholesky decomposition of SA, so thatUTAUA = SA. Similarly, UW is the Cholesky decomposition ofSW. The equivalent LGSSM is then completed by specifying˜ A ≡  A , ˜ ΣH ≡ ΣH, ˜ ΣV ≡ diag(ΣV ,I,I), ˜ µ ≡ µ, ˜ Σ ≡ Σ6.In this way any standard inference routines in the literaturemay be applied to compute q(ht) = ˜ q(ht|˜ v1:T), includingthose speciﬁcally addressed at improving numerical stability[9]. In the experiments, we used the standard predictor-corrector ﬁltering and Rauch-Tung-Striebel smoothing [9]. Aminor modiﬁcation to the standard predictor-corrector ﬁlteringroutine may be applied for computational efﬁciency (see [8]for details). This method is considerably simpler and moregeneral than the procedure given in [7], which is based onBelief Propagation and do not correspond to any of thestandard forms in the Kalman ﬁltering/smoothing literature.Finding the Optimal ΘDifferentiating F with respect to Θ we ﬁnd that, optimally:βj =VPi￿W2ij￿q(W), αc =H2cPi,jD[Ac− ˆ Ac]2ijEq(Ac),ΣV =1TT Xt=1D(vt−WPht)(vt−WPht)TEq(W)q(ht),ΣcH =1T−1T Xt=2D￿hct−Achct−1￿￿hct−Achct−1￿TEq(Ac)q(hct−1:t),Σ =D(h1−µ)(h1−µ)TEq(h1), µ =  h1 q(h1) .A. DemonstrationIn a proof of concept experiment, we used a FLGSSM togenerate 3 sources sct with random 5 × 5 transition matricesAc, µ = 0H and Σ ≡ ΣH ≡ I, see Fig. 1a. The sources6Strictly, we need a time-dependent emission ˜ Bt = ˜ B, for t = 1,...,T −1. For time T, ˜ BT has the Cholesky factor UA replaced by a zero matrix.were mixed into three observations vt = Wst + ηvt , for Wchosen with elements from a zero mean unit variance Gaussiandistribution, and ΣV = I (Fig.1b). We then trained a BayesianFLGSSM with 5 sources and 7 × 7 transition matrices Ac.To bias the model to ﬁnd the simplest sources, we used zeromatrices ˆ Ac for all sources. In Fig. 1c we plot the estimatedsources from our method after convergence. Two of the 5sources have been removed, and the remaining three are areasonable estimation of the original sources. Another possibleapproach for introducing prior knowledge is to use a Maxi-mum a Posteriori (MAP) procedure by adding a prior term tothe original log-likelihood logp(v1:T|A,W,Θ)+logp(A|α)+logp(W|β). However, it is not clear how to reliably ﬁnd thehyperparameters α and β in this case. One solution is toestimate them by optimizing the new objective function jointlywith respect to the parameters and hyperparameters (this is theso-called joint map estimation – see for example [10]). Thecomplexity of this approach is similar to the unaugmentedBayesian LGSSM, although in this case solving a Sylvesterequation is required for updating the parameters. A typicalresult of using this joint MAP approach on the artiﬁcial datais presented in Fig. 1d. The joint MAP does not estimate thehyperparameters well, so that an incorrect number of sourcesis found, and the sources are not as well estimated as in theBayesian procedure.B. Application to EEG AnalysisIn Fig. 2a we plot three seconds of EEG data recordedfrom 4 channels (located in the right hemisphere) while aperson is performing imagined movement of the right hand.As is typical in EEG, each channel shows drift terms below1 Hz which correspond to artifacts of the instrumentation,together with the presence of 50 Hz mains contamination.These effects mask the rhythmical activity related to themental task, mainly centered at 10 and 20 Hz, which wewant to extract. Standard ICA methods such as FastICAdo not ﬁnd satisfactory sources based on raw ‘noisy’ data,and preprocessing with band-pass ﬁlters is usually required.Additionally, in EEG research, ﬂexibility in the number ofrecovered sources is important, since there may be manyindependent oscillators of interest underlying the observationsand we would like some way to automatically determine theireffective number. To preferentially ﬁnd sources at particularfrequencies, we speciﬁed a block diagonal matrix ˆ Ac witheach block being a rotation at the desired frequency ω: ￿cos(2πω/N) −sin(2πω/N)sin(2πω/N) cos(2πω/N)￿, where N is the numberof samples per second. In order to extract the dominant driftsbelow 1 Hz, the mains contaminations and the oscillationsrelated to the mental task, we deﬁned the following 16groups of frequencies ω: [0.5], [0.5], [0.5], [0.5]; [10,11],[10,11], [10,11], [10,11]; [20,21], [20,21], [20,21], [20,21];[50], [50], [50], [50]. Hence, the total hidden dimension of theFLGSSM is H = 48. The temporal evolution of the sourcesobtained after training the Bayesian FLGSSM is shown inFig. 2b (grouped by frequency range). This method removed4 unnecessary sources from the mixing matrix W, that is one[10,11] Hz and three [20,21] Hz sources. We can see that theSIGNAL PROCESSING LETTERS 40   1 2 3(b)0   1 2 3(c)0   1 2 3(a)0   1 2 3(d)Fig. 2. (a) Three seconds of unﬁltered EEG data recorded from 4 electrodes.(b) The 16 sources st estimated by the Bayesian FLGSSM. (b) Sourcesestimated by the MAP FLGSSM. (c) The 16 factors estimated by NDFA.The retained sources have been rescaled to aid visualization.ﬁrst 4 sources (counting from the top down) contain dominantlow frequency drift, source 5, 6 and 8 contain [10,11] Hz,while source 10 contains [20,21] Hz centered activity. Outof the 4 sources initialized to 50 Hz, only 2 retained 50 Hzactivity, while the Ac of the other two have changed to modelother frequencies present in the EEG.The MAP FLGSSM approach is presented in Fig. 2c. Wecan see that none of the [10,11] Hz sources has been removed,even if contribution of source 8 to the observations is relativelysmall. One Ac biased at [10,11] Hz includes other frequenciesin addition to 10 Hz (source 7). As in the Bayesian case, onlyone [20,21] Hz component is retained. There are two dominant50 Hz components, however none of 50 Hz sources has beenremoved. In conclusion, the Bayesian FLGSSM seems betterable to remove unnecessary components and gives cleanersources at the desired frequencies.To asses the advantage of using prior frequencies for ex-tracting task-related information and the potential limitationsof using a linear model, we compared our method with NDFA[4]. We extracted 16 factors using a NDFA model in whichboth MLPs had one hidden layer of 30 neurons. In Fig. 2dwe show the temporal evolution of the resulting factors. Theﬁrst 10 factors from the top give the strongest contribution tothe observations. In agreement with the Bayesian FLGSSM,there are 2 main 50 Hz sources (ﬁrst two factors), although asmall 50 Hz activity is present also in other factors, namely7, 11, 12 and 14. The slow drift has not been isolated andis present in almost all factors. The information related tohand movement, namely [10,20] Hz activity, is spread overfactors 3, 4, 9, 10 and 13, which however contain also otherfrequencies. The prior speciﬁcation of independent dynamicalprocesses at particular frequencies has therefore helped theBayesian FLGSSM to better isolate the activity of interestinto a smaller number of sources and, among these sources,to separate the contribution of oscillators at 10 Hz and 20 Hz.III. CONCLUSIONWe presented a method to identify independent dynamicalsources in noisy temporal data, based on a Bayesian procedurewhich automatically biases the solution to ﬁnding a smallnumber of sources with preferential dynamics. This proce-dure is closely related to others previously proposed in theliterature, but has the property that the sources are themselvesprojections from higher dimensional independent linear dy-namical systems. Here we concentrated on the projection to asingle dimension since this aids interpretability of the signals,being of particular importance for applications in biomedicalsignal analysis. An advantage of our linear dynamics approachis tractability of inference, and we demonstrated how thestatistics of the hidden variables in the Bayesian LGSSM canbe estimated by using any Kalman ﬁltering/smoothing routine.The method is able then to automatically extract signals, forexample, biased towards particular frequencies.REFERENCES[1] A. Hyv¨ arinen, J. Karhunen, and E. Oja, Independent Component Anal-ysis. John Wiley and Sons, 2001.[2] K. Knuth, “Bayesian source separation and localization,” in SPIE’98:Bayesian Inference for Inverse Problems, 1998, pp. 147–158.[3] J. S¨ arel¨ a and H. Valpola, “Denoising source separation,” Journal ofMachine Learning Research, vol. 6, pp. 233–272, 2005.[4] H. Valpola and J. Karhunen, “An unsupervised ensemble learningmethod for nonlinear dynamic state-space models,” Neural Computation,vol. 14, pp. 2647–2692, 2002.[5] J. S¨ arel¨ a, H. Valpola, R. Vig´ ario, and E. Oja, “Dynamical FactorAnalysis of Rhythmic Magnetoencephalographic Activity,” in 3rd In-ternational Conference on Independent Component Analysis and BlindSignal Separation, 2001, pp. 457–462.[6] Y. Bar-Shalom and X. Li, Estimation and Tracking: Principles, Tech-niques and Software. Artech House, 1998.[7] M. Beal, “Variational Algorithms for Approximate Bayesian Inference,”Ph.D. Thesis, Gatsby Computational Neuroscience Unit, UniversityCollege London, UK, 2003.[8] S. Chiappa, “Analysis and classiﬁcation of EEG signals using prob-abilistic models for brain computer interfaces,” Ph.D. Thesis, EPFL,Lausanne, Switzerland, 2006.[9] M. Grewal and A. Andrews, Kalman Filtering: Theory and PracticeUsing MATLAB. John Wiley and Sons, Inc., 2001.[10] S. Saquib, C. Bouman, and K. Sauer, “ML parameter estimation forMarkov random ﬁelds with applicationsto Bayesian tomography,” IEEETransactions on Image Processing, vol. 7, pp. 1029–1044, 1998.

An unsupervised ensemble learning method for nonlinear dynamic state-space models,”

Analysis and classiﬁcation of EEG signals using probabilistic models for brain computer interfaces,” Ph.D. Thesis,

Bayesian source separation and localization,”

Denoising source separation,”

Dynamical Factor Analysis of Rhythmic Magnetoencephalographic Activity,”

Independent Component Analysis.

Kalman Filtering: Theory and Practice Using MATLAB.

ML parameter estimation for Markov random ﬁelds with applicationsto Bayesian tomography,”

Variational Algorithms for Approximate Bayesian Inference,”

Bayesian factorial linear Gaussian state-space models for biosignal decomposition

Bayesian factorial linear Gaussian state-space models for biosignal decomposition

Abstract

Similar works

Full text

Available Versions

UCL Discovery