Blind Source Separation of Fmri Signals Using Joint Diagonalization Algorithm

Blind Source Separation (BSS) is a model free source separation technique which decomposes observed mixture data into mixing matrix and source matrix both of which are unknown beforehand. One well known BSS algorithm is joint diagonalization which is from the algebraic class and in which mixing structures are recovered by jointly diagonalizing the source condition matrix. In this study we first review the existing joint diagonalizing algorithm and then propose a modified high order and exponential gradient of the algorithm. The proposed algorithm is tested on simulated images and synthetic fMRI signals. Quality and execution time of the extracted sources and time courses is compared with conventional JD algorithm.


INTRODUCTION
Different techniques are used for detection of brain functionality including Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT) and fMRI etc. fMRI is the latest and comparatively accurate non-invasive clinical and research technique in which the magnetic properties of the brain tissues associated with activity and nonactivity are recorded in the form of 3-D images of the brain (Ogawa et al., 1990).Since the requirement of oxygenated blood for active neurons is more than the inactive neurons, thus there is a considerable difference of blood flow and consequently magnetic properties of them are different.When a set of neurons becomes active, it needs more oxygenated blood as compared to nearby idle neurons.Oxygenated blood consists of more iron which change its magnetic properties.The change in magnetic properties of the oxygenated blood is recorded by the fMRI scanner in the form of Blood Oxygen Level Dependent (BOLD) signal (Ogawa et al., 1990).Statistical analysis of this multi image data is done for finding activity and non activity.The results of the analysis are normally presented on a statistical map.However due to small effect of BOLD variations in the signal intensity, SNR values of functional MR images are too low (Friman et al., 2001) and thus a de-noising step is normally performed before classification (Amir et al., 2012).
During fMRI experiment the subjects are normally asked to perform some activity like movement of hand etc and then quiet for some time.This cycle of activity and non-activity is repeated few times.BOLD signal is recorded by the fMRI scanner in the form of 3D images.During this trial period thousands of volume elements in each brain slice are acquired producing time series which are mixtures of activity, non activity, noises and structural information.It is not easy to localize active brain spots in this mixture (McKeown et al., 1998).Variations in BOLD signal has been analyzed by different approaches including time frequency analysis, statistical parametric mapping, correlation analysis, principal component analysis and Independent component analysis etc.
Since spectral characteristics of activity and non activity are different, thus regions of activity and non activity can be pointed out using time frequency analysis (Mitra et al., 1997).T-test and F-test are done to decide about each voxel in Statistical Parametric Mapping (SPM) which is a univariate approach (Friston, 1996).In case of correlation analysis decisions about voxels are made on the basis of experimental model.Each brain voxel is correlated with known a priori model and decision is made about the activity or non-activity of the voxel depending upon the correlation results (Bandettini et al., 1993).Methods discussed so far ie SPM, time frequency analysis and correlation analysis are known as hypothesis driven approaches because they are based on experimental model of the data.However data driven approaches of fMRI analysis do not require any experimental model like Principal Component Analysis (PCA) (Backfrieder, 1996) and Independent Component Analysis (ICA) (McKeown et al., 2003).In principal component analysis which is a data driven approach, eigen images are find out from the data having greatest variance.
PCA is also applied for data reduction in high dimensional data.There exists other multi-variate data processing approaches, but the most prominent approach is ICA, which is a powerful exploratory tool for fMRI data analysis (Qiu-Hua et al., 2010;Correa et al., 2007).ICA converts a multidimensional vector into statistically independent components.Therefore this algorithm is broadly used in other BSS problems as well.
However, ICA requires that the hidden source should be independent and which is the main limitation of ICA.Sometime this may not be the case and the sources may not be independent.In such cases other approaches are used for processing of fMRI data including NMF and PCA.However for PCA application the data is assumed to be uncorrelated spatially and temporally (Xiaoxiang et al., 2004).In case of NMF non-negitivity of data is the requirement (Lee and Seung, 1999).In this study we are going to apply joint diagonalization algorithm for the source extraction of functional MRI data with high order and exponential contrast function for speed and accuracy.

BSS problem:
Let ܺ is the observation matrix such that ܺ = ሾ‫ܣ‬ሿሾܵሿ.Normally the goal of BSS is to find the mixing matrix ‫ܣ‬ and source matrix ܵ .Now if ‫ܣ‬ is find out, ܵ can be estimated by ܵ = ሾܹሿሾܺሿ where ܹ = ‫ܣ‬ ற .
Here symbol † denotes pseudo-inverse of ‫,ܣ‬ which equals the inverse if numbers of rows and columns in ‫ܣ‬ are equal.Since in this case we are applying dimension reduction on fMRI data and thus rows becomes equal to columns therefore pseudo inverse is equal to the inverse.Thus the BSS problem is reduced to find ‫ܣ‬ only.
In case of fMRI data we normally assume that the observation matrix ܺ is a linear mixture of the source vector ܵ and the time courses matrix ‫ܣ‬ such that: In which ܺ is the observed data of dimension n×m, ‫ܣ‬ is a mixing matrix of dimension k n× and ܵ consist of independent sources of dimension k×m.
The aim of any BSS algorithm is to find the unmixing matrix ܹ such that: where, ܹ = ‫ܣ‬ ିଵ with the constraint that ܹ ் ܹ = ‫.ܫ‬In most of the BSS algorithms, it is required that the data should be pre processed (centered and white).The mean is subtracted from the data to make it centered as in Eq. ( 3 where, X and X' represents the observed data and centered data, respectively. Eigen decomposition is used for whitening the data.A whitening matrix v is defined and multiplied with the observed data so that the correlation and covariance matrix of sources becomes identity matrix i.e., E [YY T ] = I.To find v, eigen decomposition is used (Xiaoxiang et al., 2004) and is given as under: where, Since fMRI is a high dimensional data.Therefore, another preprocessing step is also required i.e., the dimension reduction step.Here we have used Singular Value Decomposition (SVD) as a dimension reduction technique (Xiaoxiang et al., 2004).Now the data is ready and the independent sources can be extracted by using the BSS algorithm.
JD algorithm: Joint diagonalization based BSS algorithms normally employ diagonalization techniques on some source conditioning matrix for the identification of mixing matrix.Now if a source of symmetric matrices is given like ܺ = ‫ݔ{‬ ଵ , ‫ݔ‬ ଶ , ‫ݔ‬ ଷ, ‫ݔ‬ ସ … … }, then the requirement is to find the invertible orthogonal matrix ‫ܣ‬ such that: Such mixing matrix ‫ܣ‬ can be found by minimizing ݂ሺ‫ܣ‬ሻ = ‫ܣ‪݂݂ሺ‬‬ ் ܴ ௫ ‫ܣ‬ሻ with respect to ‫.ܣ‬Here ܴ ௫ is the autocovariance matrix of ܺ and off denotes the off diagonal elements (Fabian et al., 2005).
A global minima of such a cost function is called joint diagonalizer of ܴ ௫ .Application of this cost function for BSS problem can be seen in Cardoso and Souloumiac (1995) and Yeredor (2002).
Proposed higher order modified JD: Since X = ሾAሿሾSሿ Taking covariance of both sides: where, R ୶ and R ୱ are covariance matrices of observed data matrix and hidden source matrix.Searching for sources S where S = ሾWሿሾXሿ leads to R ୱ = WR ୶ W .
We know that all sources are uncorrelated therefore EሾSS ሿ = R ୱ = I.
To find W iteratively here we suggest to use fixed point iteration algorithm with 3 rd order contrast function.High order contrast functions converge the algorithm rapidly and the global maxima or minima is achieved with high accuracy and speed:

Simulated fMRI data:
To evaluate the performance of the suggested algorithm, synthetic fMRI data is used.This data is freely available on web http:// mlsp.umbc.edu/simulated_fmri_data.html.This synthetic data was initially synthesized by Correa et al. (2007) for testing the performance of their algorithm.This data consists of eight sources and their corresponding time series as shown in Fig. 2. The first source and corresponding time series is quitec prominent in Fig. 2 which represents activated voxels.Source 2 and 6 represents the physiological functions, while remaining sources and time courses show artifacts and noises.Observed data is the mixture of these eight sources and time courses.Thus forming 100 images.Four sample mixed images are shown in Fig. 3.The mixed data is then preprocessed for further analysis.

RESULTS AND DISCUSSION
Four Simulated signals and its mixtures are shown in Fig. 1.Our goal is to recover the source signals from its mixtures using conventional JD algorithm and then using the proposed algorithms and see the accuracy of the results by comparing the extracted sources with the simulated sources using correlation.First row of Fig. 4 shows the recovered sources by conventional JD algorithm.It can be seen that source 2 and 4 are Fig.4: Recovered sources by Conventional, Proposed 3 rd order and Exponential JD algorithm recovered 99 and 96% accuracy, respectively, however source 1 and 3 are recovered with 24 and 0.3% accuracy.The overall accuracy of conventional JD is 56%.Second row of Fig. 4 shows the extracted sources using 3 rd order proposed JD algorithm.It can be seen that the results have one to one correspondence with the simulated sources of Fig. 1.Accuracy of source 1, 2, 3, 4 99% for all sources and thus average accuracy of all sources is also 99%.Third row of Fig. 4 shows sources extracted by exponential contrast function JD algorithm with the accuracy of 99% as well.Correlation results of the extracted sources with the simulated sources is shown in Table 1.It is now clear that proposed algorithm accuracy is better than conventional JD algorithm, therefore it can now be applied to simulated fMRI data.For this purpose a standard simulated fMRI    2. On average the algorithm recovers sources and time courses with the accuracy of 68 and 48%, respectively.Figure 6 and 7 shows the extracted sources and Time courses recovered by proposed 3 rd order and Exponential contrast function algorithm.Accuracy of the 3 rd order and exponential algorithm for sources/time courses is 77/76 and 80/76%, respectively.

CONCLUSION
In this study joint diagnalization algorithm is reviewed and applied to BSS problem.Third order and exponential contrast function is used for faster convergence and accuracy.Conventional and proposed algorithms are applied to simulated signals and standard simulated fMRI data and the recovered results were compared for accuracy and time of execution.By looking into Table 1 and 2 it can be concluded that the proposed algorithms performed well not only for general but also for fMRI data source separation.

Fig. 3 :
Fig. 3: Simulated fMRI mixtures Simulated signals: First row of Fig. 1 shows four simulated source signals.These signals are mixed by linearly combining them with random signals thus making eight mixed signals as shown in row 2 and 3 of Fig. 1.

Fig. 5 :
Fig. 5: Recovered fMRI sources and time courses by conventional JD

Fig. 7 :
Fig. 7: Recovered sources and time courses by the proposed exponential JD algorithm data is used having eight sources and corresponding time courses.These sources are mixed with Time course matrix thus forming hundred images, four of them are shown in Fig. 3. To recover fMRI sources and time courses from these mixture images, first conventional JD algorithm is applied, the results of which are shown in Fig. 5. Correlation results of the extracted sources and time courses with actual sources and time courses of Fig. 2 are shown in Table

Table 1 :
Execution time and correlation results of extracted sources with actual sources (simulated signals)

Table 2 :
Execution time and correlation results of extracted sources/time courses with actual sources/time courses (fMRI simulated data)