6 Discussion

We have presented an iterative rank-1 tensor-PICA decomposition for the analysis of single group FMRI data. The method was derived from the three-way PARAFAC model by adding additional maximum non-Gaussianity constraints to the estimated spatial maps. The result of this constraint for estimates in the spatial domain is that the tensor-PICA approach no longer treats all modes of variation as equal. This is an important aspect of the tensor-PICA model, since in FMRI there are substantially different numbers of observations available in the different domains, i.e. a typical FMRI group study involves $10-30$ subjects, with $50-300$ volumes and $25000$ - $45000$ intra-cranial voxels (after co-registration into a common space). The tensor-PICA approach, unlike PARAFAC, places stronger statistical constraints on the spatial domain where plentiful data is available.

The approach differs from existing group-ICA methodology [Calhoun et al., 2001,Lukic et al., 2002,Svensén et al., 2002,Leibovici and Beckmann, 2001 (sec. 4)] in that it does not simply concatenate the data in space or time in order to perform a single two-dimensional ICA decomposition followed by some meta-analysis to estimate the variation between subjects. Instead, the tensor-PICA approach directly estimates separate modes in the three domains by iterating between estimating a 2-D PICA model on the data concatenated in time and a rank-1 decomposition of the resulting estimate of the mixing matrix $\mbox{\protect\boldmath$M$}^c$. Due to the iterative nature, the three-dimensional nature of the data is represented within the estimation stage. This eliminates the need for heuristic post-processing of a set of time-courses (or a set of spatial maps) in order to express the variation across subjects. The technique is fully automated, including the estimation of the model order $R$.

In this paper, we have concentrated on the case of a single group. The methodology can, however, be extended to higher dimensions: under the model of equation 1 where we assume the existence of a single group there is only one single non-zero Eigenvalue for each matrix ${\mbox{\protect\boldmath$M$}}^c_r$ and a rank-1 approximation is appropriate. In practical applications with finite observations and in the presence of noise, the matrices ${\mbox{\protect\boldmath$M$}}^c_r$ will be of full rank. If a sufficient number of observations in the session/subject domain is available, we can apply model order selection techniques and estimate the number of time courses which combine to represent the temporal characteristics of the sources. In the case of 2 groups with similar spatial signal but sufficiently different temporal characteristics, a rank-2 approximation to each matrix ${\mbox{\protect\boldmath$M$}}^c_r$ will then result in a 4-way decomposition of the data. The typically small number of observations in the session/subject domain, however, makes estimation of the model order from the data very difficult. It is possible, however, to impose the number of different time courses that are used to represent the temporal characteristics of all sessions/subjects/groups, e.g. one can use a rank-2 approximation to ${\mbox{\protect\boldmath$M$}}^c_r$ in order to estimate different temporal responses from 2 different subgroups in the population.

For the generative model of equation 1 we have demonstrated, on a set of artificial group data sets, that tensor-PICA can successfully estimate multiple processes in the spatial, temporal and subject/session domain. Compared to a PARAFAC decomposition, the tensor-PICA estimation shows significant improvements in the form of: (i) an increased accuracy for primary 'activation' maps, (ii) reduced cross-talk between the different estimated spatial maps and (iii) an increased robustness against deviation from the model assumptions and against estimating the model order $R$ incorrectly.

All of these improvements are a direct consequence of the optimisation for maximally non-Gaussian spatial source distributions. Typical FMRI 'activation' is sparse in the spatial domain and the estimated linear regression coefficients (spatial maps) will contain only a few 'significantly large' values embedded in random Gaussian distributed 'background noise'. An optimisation for non-Gaussianity of estimated spatial maps optimises for the largest possible separation of the first set of regression coefficients ('activation') from all other remaining coefficients ('background'). Unlike the sum-of-squares error function associated with the PARAFAC model, the error function associated with an optimisation for maximum non-Gaussianity does not improve from 'splitting' components and/or having multiple components which 'explain' the same signals. In particular, an optimisation for jointly maximal non-Gaussian spatial maps implies a minimization of statistical dependence and cross-talk in the spatial domain (see [Hyvärinen and Oja, 1997] for a clear account of the relation between statistical independence and non-Gaussianity).

As an additional benefit over PARAFAC, the tensor-PICA decompositions each required significantly less computation (between 1/10 and 1/100 times the number of floating point operations) compared to PARAFAC in order to converge to a solution¹¹. Again this is a consequence of the fact that the cost function in a Tensor-PICA decomposition is more sensitive to the particular signal characteristics in the spatial domain. As a result, the number of iterations until convergence is much reduced (the number of floating-point operations per iteration is actually greater when using Tensor-PICA as opposed to PARAFAC).

Using real FMRI data we have demonstrated that tensor-PICA can extract plausible spatial maps and time courses. The main activation pattern of the multi-session decomposition identified cortical regions which correspond to what has been reported in [McGonigle et al., 2000]. Furthermore, the tensor-PICA decomposition gives a rich description of additional processes in the data. For example, the tensor-PICA decomposition separated negative (de-)activation into different plausible spatial maps with associated time-courses and variation across sessions. The final decomposition does suggest that there are at least two distinct processes which contribute to the negative $Z$-scores: plausible ipsilateral de-activation in the primary motor cortex consistent across sessions and de-activation in the superior occipital lobule which only appears in a few sessions, possibly due to visual fixation. In addition, the tensor-PICA decomposition identified nuisance effects like artefactual RF signal components or stimulus correlated motion at the group level.

We believe that the tensor-PICA approach can provide simple and useful representations of multi-subject/multi-session FMRI data that can aid the interpretation and optimisation of group FMRI studies.