3 Tensor PICA

The set of equations 3-5 can alternatively be expressed as simple matrix products, e.g. the set of equations 5 can be expressed as:

$$\mbox{\protect\boldmath$X$}_{IK\times J} = (\mbox{\protect\boldmath$C$}\,\vert\!\!\otimes\!\!\vert\,\mbox{\p... ...ox{\scriptsize\textit{\sffamily {t}}}}+\widetilde{\mbox{\protect\boldmath$E$}}.$$ (8)

Here, $\mbox{\protect\boldmath$X$}_{IK\times J}$ denotes the $IK\times J$ data matrix formed by concatenating all $K$ different data sets in the temporal domain and $\mbox{\protect\boldmath$A$}, \mbox{\protect\boldmath$B$}$ and $\mbox{\protect\boldmath$C$}$ are matrices containing the $R$ individual temporal, spatial and session/subject factors in their columns. The first factor $(\mbox{\protect\boldmath$C$}\,\vert\!\!\otimes\!\!\vert\,\mbox{\protect\boldmath$A$})$ denotes the Khatri-Rao product of $\mbox{\protect\boldmath$A$}$ and $\mbox{\protect\boldmath$C$}$, i.e. a $IK\times R$ matrix formed by $K$ copies of $\mbox{\protect\boldmath$A$}$ stacked and column-wise scaled by diagonal matrices formed from rows of $\mbox{\protect\boldmath$C$}$ such that $(\mbox{\protect\boldmath$C$}\,\vert\!\!\otimes\!\!\vert\,\mbox{\protect\boldmat... ...\scriptsize\textit{\sffamily {t}}}})^{\mbox{\scriptsize\textit{\sffamily {t}}}}$  [Bro, 1998].

From equation 8, the matrix of spatial factors, $\mbox{\protect\boldmath$B$}$, has least-squares estimates of

$$\widehat{\mbox{\protect\boldmath$B$}}^{\mbox{\scriptsize\tex... ...oldmath$A$})^\dagger \mbox{\protect\boldmath$X$}_{IK\times J}.$$ (9)

Similarly, $\widehat{\mbox{\protect\boldmath$A$}}$ and $\widehat{\mbox{\protect\boldmath$C$}}$ can be estimated from equations 3 and 4. Thus, ALS estimates for each of the three matrices can be calculated via a linear projection of the three-way data, reshaped to 3 different 2-D matrices.

The PARAFAC model and the ALS algorithm for estimation treat all three domains equally and do not utilise any domain-specific information. Section 4 demonstrates how this can lead to PARAFAC results which are difficult to interpret, mainly due to significant cross-talk between estimated spatial maps.

In order to address this, we formulate a tensor-PICA model which incorporates the assumption of maximally non-Gaussian distributions of estimated spatial maps, $\mbox{\protect\boldmath$B$}$: equation 8 is identical to a standard (2-D) factor analysis or noisy ICA model [Beckmann and Smith, 2004], where the matrix $(\mbox{\protect\boldmath$C$}\,\vert\!\!\otimes\!\!\vert\,\mbox{\protect\boldmath$A$})$ denotes the 'mixing' matrix and $\mbox{\protect\boldmath$B$}^{\mbox{\scriptsize\textit{\sffamily {t}}}}$ contains the set of spatial maps as row vectors. Unlike the single subject (2-D) PICA model, however, the mixing matrix now has a special block structure which can be used to identify the factor matrices $\mbox{\protect\boldmath$A$}$ and $\mbox{\protect\boldmath$C$}$. Given the first matrix factor in equation 8, it is easy to recover the two underlying matrices $\mbox{\protect\boldmath$A$}$ and $\mbox{\protect\boldmath$C$}$: each of the $R$ columns in $(\mbox{\protect\boldmath$C$}\,\vert\!\!\otimes\!\!\vert\,\mbox{\protect\boldmath$A$})$ is formed by $K$ scaled repetitions of a single column from $\mbox{\protect\boldmath$A$}$, i.e. when reshaped into a $I\times K$ matrix is of rank 1. Thus, we can transform each column $r$ into a $I\times K$ matrix and calculate its (single) non-zero left Eigenvector of length $I$, together with a set of $K$ factor loadings (projections of the matrix onto the left Eigenvector), using a Singular Value Decomposition (SVD) and use these to re-constitute a column of the underlying factor matrices $\mbox{\protect\boldmath$A$}$ and $\mbox{\protect\boldmath$C$}$. This needs to be repeated for each of the $R$ columns separately and the matrices $\mbox{\protect\boldmath$A$}$ and $\mbox{\protect\boldmath$C$}$ are proportional to the $R$ different Eigenvectors and factor loadings respectively, i.e. the values obtained by projecting the matrix of $I\times K$ matrix of time courses onto the Eigenvector of the SVD.

This gives the following algorithm for a rank-1 tensor PICA decomposition of three-way data $\mbox{\protect\boldmath$X$}$:

1. perform an iteration step for the decomposition of the full data
   

   $$\mbox{\protect\boldmath$X$}_{IK\times J} = \mbox{\protect\bo... ...it{\sffamily {t}}}}+ \widetilde{\mbox{\protect\boldmath$E$}}_1$$ (10)

   
   
   using the 2-D PICA approach for the decomposition into a compound mixing matrix $\mbox{\protect\boldmath$M$}^c$ of dimension $IK\times R$ and a $R\times J$ matrix containing the associated spatial maps $\mbox{\protect\boldmath$B$}^{\mbox{\scriptsize\textit{\sffamily {t}}}}$.
2. decompose the estimated mixing matrix $\mbox{\protect\boldmath$M$}^c$ such that
   

   $$\mbox{\protect\boldmath$M$}^c = (\mbox{\protect\boldmath$C$}... ...rotect\boldmath$A$})+\widetilde{\mbox{\protect\boldmath$E$}}_2$$ (11)

   
   
   via a column-wise rank-1 Eigenvalue decomposition: for each column $r$ of $\mbox{\protect\boldmath$M$}^c$, form the $I\times K$ matrix $\widetilde{\mbox{\protect\boldmath$M$}}_r=(\mbox{\protect\boldmath$M$}^c_{r,1},\dots,\mbox{\protect\boldmath$M$}^c_{r,K})$. Under the model, this matrix contains $K$ scaled repetitions of a single temporal factor $\mbox{\protect\boldmath$A$}_r$ which can be found by calculating its least-squares rank-1 approximation using SVD. Along with the temporal factor $\mbox{\protect\boldmath$A$}_r$ (left Eigenvectors), the SVD provides the individual scalings (right Eigenvectors) which define the corresponding vector in $\mbox{\protect\boldmath$C$}_r$. This needs to be repeated for each column of $\mbox{\protect\boldmath$M$}^c$ and provides estimates for $\mbox{\protect\boldmath$A$}$ and $\mbox{\protect\boldmath$C$}$.

3. iterate decomposition of $\mbox{\protect\boldmath$X$}_{IK\times J}$ and $\mbox{\protect\boldmath$M$}^c$ until convergence⁵.

Note that, like PARAFAC, the rank-1 tensor PICA decomposition estimates factor matrices for the generative model of equation 1. The estimated matrices, however, provide a different structural representation of the three-way data $\mbox{\protect\boldmath$X$}$. Note also, that the singular value decomposition of each matrix $\widetilde{\mbox{\protect\boldmath$M$}}_r$ not only provides the left and right Eigenvectors which form the relevant columns in $\mbox{\protect\boldmath$A$}_r$ and $\mbox{\protect\boldmath$C$}_r$ but also gives a set of Eigenvalues. The ratio of the largest Eigenvalue and the sum of all Eigenvalues can be used to assess the quality of the rank-1 approximation: if the matrix $\widetilde{\mbox{\protect\boldmath$M$}}_r$ is not well approximated by the outer product of the left and right Eigenvectors the corresponding ratio will be low, i.e. only represent a small amount of variability in $\widetilde{\mbox{\protect\boldmath$M$}}_r$.