In this paper we derive an independent-component analysis (ICA) method for analyzing two or more data sets simultaneously. Our model permits there to be components individual to the various data sets, and others that are common to all the sets. We explore the assumed time autocorrelation of independent signal components and base our algorithm on prediction analysis. We illustrate the algorithm using a simple image separation example. Our aim is to apply this method to functional brain mapping using functional magnetic resonance imaging (fMRI). 1

Hansen, Lars Kai

Lukic, Ana S.

Strother, Stephen C.

Wernick, Miles N.

Online Research Database In Technology

AN ICA ALGORITHM FOR ANALYZING MULTIPLE DATA SETSAna S. Lukic1,  Miles N. Wernick1, Lars Kai Hansen2, Stephen C. Strother31Illinois Institute of Technology, Chicago, 2Technical University of Denmark, Lyngby, Denmark3VA Medical Center and University of Minnesota, MinneapolisABSTRACTIn this paper we derive an independent-componentanalysis (ICA) method for analyzing two or more datasets simultaneously.  Our model permits there to becomponents individual to the various data sets, andothers that are common to all the sets.  We explore theassumed time autocorrelation of independent signalcomponents and base our algorithm on predictionanalysis. We illustrate the algorithm using a simple imageseparation example.   Our aim is to apply this method tofunctional brain mapping using functional magneticresonance imaging (fMRI).1. INTRODUCTIONThere are several approaches to the problem of blindseparation of linear signal mixtures. Some of them arebased on non-Gaussianity [1, 3]  while others exploittemporal correlations [2, 5]. The technique of Molgedeyand Schuster [5] is especially attractive since it offers anon-iterative solution. The limitations of the originalalgorithm, such as restriction to square non-singularmixing matrices and inherent erroneous complex-valuedresults, were eliminated in [4].In this paper, we extend the Molgedey-Schuster ICAmethod to analyze more than one data set simultaneously.There are many potential applications for simultaneousindependent-component analysis (ICA) of multiple datasets.  For example, one might simultaneously analyze theaudio and video portions of a video sequence.Our interest is in developing methods to analyzefunctional magnetic resonance images (fMRI) of thebrain, to map spatial and temporal patterns of brainactivation.  In this application, the data can be dividedinto multiple data sets by considering different subjectsseparately, or by considering image data separately fromother measurements made during the fMRI study (e.g.,reaction time of the subject).2. EXTENDED BLIND SIGNAL SEPARATIONPROBLEMIn the classical blind signal separation problem it isassumed that there are N independent source signals inmatrix S observed through a mixing matrix A, i.e.,X AS= , (1)to obtain a matrix of observations X. The goal in thisclassical ICA problem is to recover the sources in S andthe mixing matrix A without any prior knowledge exceptthe assumption of independence of the sources in S.To analyze two data sets X and Y simultaneously, wecould lump them together into one matrix, and use themodel in (1); however, this would force the description ofthe two matrix to be based on the same set of sources.  Toallow the two data sets to have some commoncomponents, and other distinct components, we proposethe following model:X A S B TY A S B TL T xL NxN TxL K K TP T yP MyM TyP K K T× × × × ××× × × ×= += +0 5 0 5 0 5 0 5 0 50 5 0 5 0 5 0 5 0 5(2)where X and Y are two sets of observations; Sx  and Syare independent sources specific to X and Y, respectively;T is a matrix of K independent sources common to both Xand Y; xA , yA , xB , and yB  are mixing matrices.  Thenumbers of observations in X and Y are L and P,respectively; the numbers of sources in Sx  and Sy  are Mand N, respectively.   In this paper, we discuss thealgorithm in terms of only two data sets, but it is easilyextended to more than two.3. MIXING MATRICES AS EIGENVECTORMATRICESIn this section we show that the mixing matrices in thisICA problem can be obtained as eigenvectors of linearminimum mean-square error (LMMSE) predictionmatrices.Assuming the independent signals have non-vanishing time autocorrelation functions we supposethere is a LMMSE prediction matrix W such that:XYWXYττ! "$# =! "$#+ × +L P L P0 5 0 51 6(3)where τX and τY  are delayed versions of X and Yrespectively. By combining (2) and (3) we can write:A S B TA S B TW WW WA S B TA S B Tx x xy y yxxL LxyL PyxP LyyP Px x xy y yτ ττ τ++! "$#=!  "$##++! "$#× ×× ×0 5 0 50 5 0 5(4)( ) ( ) ( )( ) ( )( )( )1T T TT T T T T TT T 1 T1ˆ212t tt tt--= += +W UDV UDV UDV UDV UDV UDVUD V V V V D Uwhere we have divided the matrix W into four blocksWxx , Wxy , Wyx  and Wyy ; Sxτ , Syτ  and Tτ  are delayedversions of Sx , Sy  and T , respectively.Right multiplying (4) by SxT yields:A S S B T SA S S B T SW WW WA S S B TSA S S B TSx x x x xy y x y xxx xyyx yyx x x x xy y x y xτ ττ τT TT TT TT T++!  "$##=! "$#++! "$#(5)Because the sources in Sx , Sy  and T  are independent,S Sy xτT= 0, S Sy xT = 0 , T Sτ xT= 0 and TSxT = 0. Thus,  (5)becomes:A C W WW WA Cx x xx xyyx yyx xτ0 5 0 5000! "$# =! "$#! "$#, (6)where C S Sx x xτ τ0 5 = T  and C S Sx x x00 5 = T . The first row of(6) yieldsW A A C Cxx x x x x=−τ0 5 0 50 1 (7)which is an eigenvector equation for Wxx  sinceC Cx xτ0 5 0 50 1−  is diagonal.  The second row of (6) yields0 = W Ayx x (8)since Cx 00 5  is a full rank positive diagonal matrix.Similarly, by right-multiplying (4) by SyT we obtain:0 0= W A Cxy y y0 5 (9)W A A C Cyy y y y y=−τ0 5 0 50 1 (10)Finally, after right multiplying (4) by TT  we obtain:BB CW WW WBB CT Txyxx xyyx yyxy! "$# =! "$#! "$#τ0 5 0 50(11)or, rearranging,WBBBB C CT Txyxy! "$# =! "$#−τ0 5 0 50 1 (12)Equations (7), (10), and (12) are eigenvectorequations that provide a recipe for identifying the variousmixing matrices as eigenvectors of prediction matrices.4. THE ALGORITHMWe divide the algorithm in two parts. First, we estimatethe prediction matrix W and then we use it to estimatemixing matrices and independent sources. It is wellknown [6] that the LMMSE prediction matrix is:$ $ $W R RXY XYτ τ0 5 0 5 0 5= −0 1, (13)where $RXY 00 5  is an estimate of the autocorrelationfunction at lag 0,$RXYX YXY 00 5 = ! "$#T T , (14)and $RXY τ0 5  is an estimate of the autocorrelation functionat lag τ ,$RXYX YXYX YXY ττττ τ0 5 = ! "$# +! "$#12T T T T (15)where Xτ and Yτ  are cyclic permutations by τ steps, i.e.Xτ τ= +x n T0 51 6= B, x k0 5  denotes the kth column of X and⋅0 5T  denotes the argument modulo T. Equation (15)guarantees that $RXY τ0 5  is symmetric.  Consider thesingular value decompositions (SVDs):XYUDV! "$# =T (16)XYU D Vτττ τ τ! "$# =T (17)where eigenvector matrices contain only componentswhose corresponding eigenvectors are nonzero. SinceXX X XT T= τ τ  it follows that U Uτ =  and D Dτ = . Then:(18)We will estimate $Was the average of $W τ0 5 over a set ofvalues for τ:$ $W W UD V V V V D U= = += =−∑ ∑1 121 11t tjjtjtj jτ τ τ2 7 4 9T T T (19)where t is the number of different values of τ  we want toaverage over and τ j j t, = 1K  are those values. We arenow ready to outline the algorithm:PART I: Calculating prediction matrix $WStep 1: Perform the SVD in (16), keep only thecomponents whose eigenvalues are nonzero.Step 2: Calculate K V V V V= +=∑ τ τj jjtT T4 91Step 3: Calculate $W UDKD U= −121tTPART II: Estimating mixing matricesStep 1: Estimate $A x  by selecting only the eigenvectors ofWxx  satisfying (8) and having corresponding nonzeroeigenvalues. Denote the number of such eigenvectors(columns of $A x ) by N.Step 2: Estimate $A y  by selecting only the eigenvectors ofWyy  satisfying (9) and having corresponding nonzeroeigenvalues. Denote the number of such eigenvectors(columns of $A y ) by M.Step 3: Estimate $ $ $B B B= xT yTT as the matrix of Keigenvectors of W corresponding to the K largesteigenvalues of W where:K rank N MT T= − −X Y3 8 (20)PART III: Estimating independent sourcesStep 1: From X A BST=! "$#x xx calculate:$$$$$ $ $$STABA B ABXxxxxx xxx! "$#=! "$#! "$#−TTTT1(21)Step 2: Calculate:$$$$$ $$$STABA BABYyyyyy yyy!  "$## =!  "$##!  "$##−TTTT1 (22)Step 3: Calculate $$ $TT T=+x y2, since $Tx  and $Ty  should bethe same.5. RESULTSAt this writing, we have nearly completed ademonstration of the method using fMRI data; however,the results are too preliminary to report and thesubmission deadline is upon us.  We plan to have theanalysis complete by the time of the conference.In the meantime, to illustrate the algorithm simply weextend the toy example given in [4]. The top row of Fig. 1shows the four image components used to produce thedata set in the bottom row:Xxxxsst N=!   "$###= −−!   "$###! "$# +!   "$###+1231210 10 12 1111xx(23)Yyys t N= ! "$# = −! "$# +! "$# +1211111y1(24)0  0.10.20.30.40.50.60.70.80.91  Sx 120 40 60 8020406080-1  -0.8-0.6-0.4-0.20   0.2 0.4 0.6 0.8 1   Sx 220 40 60 8020406080-1  -0.8-0.6-0.4-0.20   0.2 0.4 0.6 0.8 1   T20 40 60 8020406080-1  -0.8-0.6-0.4-0.20   0.2 0.4 0.6 0.8 1   Sy 120 40 60 8020406080-2-10 1 2 3 X 120 40 60 8020406080-2-10 1 2 3 X 220 40 60 8020406080-2  -1.5-1  -0.50   0.5 1   1.5 2   2.5 3   X 320 40 60 8020406080-2-10 1 2 3 Y 120 40 60 8020406080-2-10 1 2 3 Y 220 40 60 8020406080Figure 1. Top row shows four independent  image components used to produce the images in the bottom row. Bottomrow shows two sets of mixtures X and Y of the four components with white Gaussian noise added. The first two components appearonly in the mixture X, the third one is common to both while the last one is specific to Y only-2.5-2  -1.5-1  -0.50   0.5 1   1.5 2   2.5 Sx 120 40 60 8020406080-3  -2.5-2  -1.5-1  -0.50   0.5 1   1.5 2   Sx 220 40 60 8020406080-3-2-10 1 2 3 4 T 120 40 60 8020406080-3-2-10 1 2 Sy 120 40 60 8020406080Figure 2. Results of the proposed algorithm applied to the X and Y mixtures shown in the bottom row of Fig.1. In this reconstructed set of independent components the first two changed places. Also, the first, second and forthcomponents are inverted compared to the original set.where x i i, = 1 3K  and yi i, = 1 2K  are two observed sets,s s t sx x y11 2 1, , ,  are four independent image components andN is the matrix of zero-mean Gaussian white noise withvariance σ n2 0 3= .  The results of our algorithm areshown in Fig. 2. The estimated mixing matrices are:$Ax-0.7260 -0.31340.6862 -0.41480.0443 0.8542=!   "$###, $ ..Ay =−! "$#0 71420 6999, $.4539.4475.4398Bx =!   "$###000, $ .4482.4466By =! "$#00(25)Comparing the original independent sources fromFig. 2 and those produced by the algorithm it can be seenthat $sx1  is the inverted and scaled version of sx2 ,  $sx2  is theinverted and scaled version of sx1 , and that $sy1  is theinverted version of sy1 . We can write:Xsst= −− −!   "$###−−!   "$###+!   "$###! "$#0 6268 10 8296 1 27241 6084 0 088612120 90780 89500 879612211. .4520. .. .$$...$xx(26)Y s t=−! "$# −! "$# +! "$#! "$#1 01000 9898120 89640 893212 1..$..$y1(27)where the estimated mixing matrices are closer to theoriginals. This illustrates that the algorithm recovers theindependent sources up to a scale and permutation, asexpected.6. DISCUSSIONThe algorithm presented while simple in theory might bedifficult to implement especially in the presence of noise.The condition crucial to determining the number ofindependent components in Sx is (8). Due to the presenceof noise and an error made by the prediction none of thecolumns in the product will be exactly zero. To overcomethe problem we set the following threshold and took onlythe columns that satisfy:W EW Eyx xyx xp l l ilLpPp l l ilLpPiLi L, ,, ,.=====∑∑∑∑∑< =111110 04 1K (28)where Ex  is the eigen vector matrix of Wxx. This meansthat we took the absolute value of W Eyx x , calculated thesum of elements in each column, normalized by dividingby the overall sum of elements, and thresholded theresulting score. We repeated the same procedure toestimate Ay .Another important fact is that (20) does not hold inthe presence of white noise so principal componentanalysis (PCA) is required as a preprocessing step, amethod proven successful in [1]. In this example we firstsmoothed the data by performing the PCA on pooled dataset X YT T  and eliminating the principal componentswhose eigenvectors do not satisfy:hni L Pi > = +2 1 122. K (29)where hi  is the ith eigen value and n2 is number  of pixelsin each image.  To demonstrate the robustness of thealgorithm to the noise we repeated the experiment withhigher noise level, σ n2 1= , and obtained the estimates ofmixing matrices similar to those in (25).7. CONCLUSIONSIn this paper we extended the blind source separationproblem and presented an ICA algorithm based on non-vanishing source autocorrelation functions for solving it.We assumed that there are more than one set ofobservations and that independent sources can becommon to all observations or specific to someobservations only. The algorithm developed is anextension of the original Molgedey-Schuster algorithmand is capable of distinguishing between the independentcomponents  specific to particular observations and thosecommon to all of them. It is very robust to the amount ofwhite noise added although PCA is required as apreprocessing step. It can easily be extended to the casewhen there are more than two sets of observations, eachcontaining a mixture of components common to allobservations and components specific to the particularsets only.8. REFERENCES[1] A. Bell and T. J. Sejnowski, “An Information-Maximization Approach to Blind Separation and BlindDeconvolution,” Neural Computation, vol. 7, pp. 1129-1159, 1995.[2] Cardoso, J.-F., A. Belouchrani, K. Abed-Meriam, and E.Moulines, “Blind Source Separation Technique UsingSecond-Order Statistics,” IEEE Transactions on SignalProcessing, vol. 45, pp. 434-444, 1997.[3] Comon, P., “Independent Component Analysis: A NewConcept?,” Signal Processing, vol. 36, pp. 287-314, 1994.[4] Larsen, J., L. K. Hansen, and T. Kolenda, “On IndependentComponent Analysis for Multimedia Signals,” inMultimedia Image and Video Processing, G. L., K. S.Y.,and J. Larsen, Eds.: CRC Press, 2000, pp. 175-200.[5] Molgedey, L. and H. Schuster, “Separation of IndependentSources Using Time-Delayed Correlations,” PhysicalReview Letters, vol. 72, pp. 3634-3637, 1994.[6] Proakis G. John and M. G. Dimitris, Digital SignalProcessing: Principles, Algorithms and Applications, Thirded. Upper Saddle River, New Jersey: Prentice-Hall, 1996.

An ica algorithm for analyzing multiple data sets

Abstract

Similar works

Full text

Available Versions

Online Research Database In Technology

CiteSeerX