Comparing representational geometries using whitened unbiased-distance-matrix similarity — part 5/6

Cleaned-text code points [62698, 78232). Bibliography omitted here; source.md includes it.

𝚺 ^ P = 1 M ⁡ ( N m − K m ) ​ ∑ m = 1 M 𝐑 m T ​ 𝐑 m , \hat{\bm{\Sigma}}_{P}=\frac{1}{M(N_{m}-K_{m})}\sum_{m=1}^{M}\mathbf{R}_{m}^{T}\mathbf{R}_{m}, (30)

where N m N_{m} is the number of observations, and K m K_{m} is the number of regressors of interest per partition. Oftentimes, we have the case that P > N P>N , which renders the estimate non-invertible. Even with N > P N>P , it is usually prudent to regularize the estimate, as it stabilizes the distance estimates. A practical way of doing this is to shrink the estimate of the covariance matrix to a diagonal version of the raw sample estimate:

𝚺 ~ P = h ​ diag ​ ( 𝚺 ^ P ) + ( 1 − h ) ​ 𝚺 ^ P . \tilde{\bm{\Sigma}}_{P}=h\,\text{diag}(\hat{\bm{\Sigma}}_{P})+(1-h)\,\hat{\bm{\Sigma}}_{P}. (31)

The scalar h determines the amount of shrinkage, with 0 corresponding to no shrinkage and 1 to only using the diagonal of the matrix (univariate prewhitening). Estimation methods for the optimal shrinkage coefficient have been proposed [ 30 ] , but in practice values in the range of h = [ 0.2 − 0.4 ] h=[0.2-0.4] perform well for fMRI data. The estimate is then used to obtain a spatially prewhitened versions of 𝐁 ^ m \hat{\mathbf{B}}_{m} :

𝐁 ∗ ^ m = 𝐁 ^ m 𝚺 ~ P − 1 / 2 \hat{\mathbf{B}^{*}}_{m}=\hat{\mathbf{B}}_{m}\tilde{\bm{\Sigma}}_{P}^{-1/2} (32)

Biased and unbiased estimates of the Mahalanobis distance can then be calculated via Eq. 4 , 5 , using 𝐁 ^ m ∗ \hat{\mathbf{B}}^{*}_{m} instead of 𝐁 ^ m \hat{\mathbf{B}}_{m} . To obtain a full expression of the variance-covariance matrix of this distance, we need to know the mean and covariance of the pre-whitened data. If whitening would work perfectly, the data would be independent across voxels. However, given that we operate with an estimate of 𝚺 P \bm{\Sigma}_{P} , this is not the case. Rather, the pre-whitened data will have matrix-normal distribution

𝐁 ∗ ^ m \displaystyle\hat{\mathbf{B^{*}}}_{m} ∼ ℳ 𝒩 ( 𝐁 𝚺 ~ P − 1 / 2 , 𝚺 K , 𝚺 R ) \displaystyle\sim\mathcal{MN}(\mathbf{B}\tilde{\bm{\Sigma}}_{P}^{-1/2},\bm{\Sigma}_{K},\bm{\Sigma}_{R}) (33) 𝚺 R \displaystyle\bm{\Sigma}_{R} = 𝚺 ~ P − 1 / 2 𝚺 P 𝚺 ~ P − 1 / 2 \displaystyle=\tilde{\bm{\Sigma}}_{P}^{-1/2}\bm{\Sigma}_{P}\tilde{\bm{\Sigma}}_{P}^{-1/2} (34)

The covariance matrix of these distance estimates is given by Eq. 27 , 29 , with 𝐁 \mathbf{B} replaced by 𝐁 𝚺 ~ P − 1 / 2 \mathbf{B}\tilde{\bm{\Sigma}}_{P}^{-1/2} , and 𝚺 P \bm{\Sigma}_{P} with 𝚺 R \bm{\Sigma}_{R} . The covariance structure under the assumption that 𝐁 = 𝟎 \mathbf{B}=\mathbf{0} , however, will be the same as for the Euclidean distance - therefore the whitened RDM correlation and the WUC can be use equivalently for the biased and unbiased estimates of the Mahalanobis distance.

### 4.5 Simulations

To evaluate different ways of comparing RDMs, we conducted a range of simulations, each from a known true model. Each model specified a second moment matrix. For every simulated data set, we first generated the matrix 𝐁 \mathbf{B} of true activity patterns from a multivariate normal distribution, such that the second-moment matrix ( 𝐁 T ​ 𝐁 / P \mathbf{B}^{T}\mathbf{B}/P ) was exactly 𝐆 ​ s \mathbf{G}s , with 𝐆 \mathbf{G} being the second moment matrix specified by the model and s s the scalar signal strength. The true patterns were independent across voxels. We then generated a set of M M measured activity patterns with

𝐁 ^ m = 𝐁 + 𝐄 m , \hat{\mathbf{B}}_{m}=\mathbf{B}+\mathbf{E}_{m},

where the noise 𝐄 m \mathbf{E}_{m} was drawn from a matrix normal distribution with mean zero, column covariance 𝚺 P \bm{\Sigma}_{P} and row covariance 𝚺 K \bm{\Sigma}_{K} . For most simulations, the noise was independent across channels 𝚺 P = 𝐈 P \bm{\Sigma}_{P}=\mathbf{I}_{P} , and observations 𝚺 K = 𝐈 K \bm{\Sigma}_{K}=\mathbf{I}_{K} . The influence of spatial noise correlations is studied in Appendix 6.3, and the influence of condition-wise correlations in the example in Figure 4a.

Each simulated data set was compared to the different candidate models using the different criteria. We then counted how often each method decided for the true (i.e., data generating) model over the competing model.

For the results shown in Figure 4, we used 2 simple models, each specifying different dissimilarities across 4 conditions. Data sets were generated from each of the two models with 2-12 partitions. The variance of the noise of the simulation was set to be proportional to the number of partitions, such that the variance of the averaged activity patterns was always constant. The noise was assumed to be independent across the P = 50 P=50 voxels. In simulation for Fig 4a, the measurement noise was correlated for neighboring conditions with r = 0.15 r=0.15 , for all other simulations it was independent.

For the simulations shown in Figure 6, we used three example experiments from published fMRI studies. The first two examples come from a paper investigating the representational structure of finger movements in primary motor and sensory cortex [ 9 ] . In Experiment 1, the activity patterns for K=5 fingers were measured. There were two candidate models, one that predicted the similarity structure based on the natural statistics of movement, the other based on the similarity of muscle activity patterns. The RDM correlation between the two models was relatively high ( r = 0.85 r=0.85 ).

Experiment 2 (Experiment 3 in [ 9 ] ), compared 31 different finger movements, which span the whole space of possible “piano-chord” combinations. Again, the predictions of the natural statistics and muscle activity model were compared.

The third example uses an experiment investigating the response of the human inferior temporal cortex to 96 images, including animate and inanimate objects [ 20 ] . The model predictions are derived from a convolutional deep neural network model – with each of the 7 layers providing a separate representational model. The bitmap images were presented to the deep neural network and the internal activity patterns used as representational models.

All data for Figure 6 were simulated with 8 runs, 160 voxels, and independent noise on the observations. The noise variance was set to σ 2 = 1 \sigma^{2}=1 . We first normalized the model predictions, such that the norm of the predicted squared Euclidean distances was 1. We then varied the strength of the signal systematically from 0 (pure noise data) to a level that achieved reasonably high accuracy. We generated 3,000 data sets for each experiment, parameter setting, and model. For Experiment 3, where there were 7 alternative models, we generated data sets from each of the models. We then decided whether the data was better fit by the data-generating or one of the alternative models. Accuracy was then averaged over all possible model pairs. Thus, for all 3 Experiments, chance performance was at 0.5. The code for the simulations presented in this paper is available at https://github.com/rsagroup/rsaModelComparison .

## 5 acknowledgements

This study was supported by the Discovery Grant from the National Science and Engineering Research Council (NSERC RGPIN-2016-04890) to J.D., the NSERC Discovery Grant (NSERC RGPIN-2019-06742) to M.M., the Canada First Research Excellence Fund (BrainsCAN) to Western University, and the German Research Foundation (DFG SCHU 3351/1-1) to H.S.

## 6 Appendix

### 6.1 Derivation of variance-covariance matrix of the distance estimate

#### 6.1.1 Expectations of products of random variables

The variance of distance estimates can be derived from the basic expectations of products of random variables. The expected value of the product any pair of random variable x , y x,y

E ⁡ ( x ⋅ y ) = E ⁡ ( x ) ⋅ E ⁡ ( y ) + cov ⁡ ( x , y ) \mathrm{E}\left(x\cdot y\right)=\mathrm{E}(x)\cdot\mathrm{E}(y)+{\rm cov}(x,y) (35)

If u , v , x , y u,v,x,y are jointly normally distributed, we have the following general expectations:

cov ⁡ ( x ​ y , u ​ v ) = E ⁡ ( x ) ​ E ​ ( u ) ​ cov ​ ( y , v ) + E ⁡ ( x ) ​ E ​ ( v ) ​ cov ​ ( y , u ) + E ⁡ ( y ) ​ E ​ ( u ) ​ cov ​ ( x , v ) + E ⁡ ( y ) ​ E ​ ( v ) ​ cov ​ ( x , u ) + cov ⁡ ( x , u ) ​ cov ​ ( y , v ) + cov ⁡ ( x , v ) ​ cov ​ ( y , u ) \begin{split}{\rm cov}\left(xy,uv\right)&=\mathrm{E}(x)\mathrm{E}(u){\rm cov}(y,v)+\mathrm{E}(x)\mathrm{E}(v){\rm cov}(y,u)\\ &+\mathrm{E}(y)\mathrm{E}(u){\rm cov}(x,v)+\mathrm{E}(y)\mathrm{E}(v){\rm cov}(x,u)\\ &+{\rm cov}(x,u){\rm cov}(y,v)+{\rm cov}(x,v){\rm cov}(y,u)\end{split} (36)

#### 6.1.2 Expectations of the product of normal matrices

From Eq. 36 , we can derive the basic expectations on the product of normal random matrices. Let us assume that vector 𝐱 \mathbf{x} has multi-variate normal distribution with mean 𝝁 \bm{\mu} and variance-covariance matrix 𝐕 \mathbf{V} . To derive Eq. 27 and 29 , we require the following results for the outer product 𝐱𝐱 T \mathbf{xx}^{T} . The mean is given by

E ⁡ ( 𝐱𝐱 T ) = 𝝁 ​ 𝝁 T + 𝐕 . \rm E(\mathbf{x}\mathbf{x}^{T})=\bm{\mu}\bm{\mu}^{T}+\mathbf{V}. (37)

The variance-covariance matrix of the diagonal of 𝐱𝐱 T \mathbf{xx}^{T} is

Var ⁡ ( diag ⁡ ( 𝐱𝐱 T ) ) = 4 ​ 𝝁 ​ 𝝁 T ∘ 𝐕 + 2 ​ ( 𝐕 ∘ 𝐕 ) \rm Var(\rm diag(\mathbf{xx}^{T}))=4\bm{\mu}\bm{\mu}^{T}\circ\mathbf{V}+2(\mathbf{V}\circ\mathbf{V}) (38)

These results can be easily extended to the distribution of the matrix product 𝐗𝐗 T \mathbf{X}\mathbf{X}^{T} , where 𝐗 \mathbf{X} is a random N × P N\times P matrix with independent normally-distributed columns, i.e. with matrix normal distribution 𝐗 ∼ ℳ ​ 𝒩 ​ ( 𝐌 , 𝐕 , 𝐈 ) \mathbf{X}\sim\mathcal{MN}(\mathbf{M},\mathbf{V},\mathbf{I}) .

E ⁡ ( 𝐗𝐗 T ) \displaystyle\rm E(\mathbf{X}\mathbf{X}^{T}) = \displaystyle= 𝐌𝐌 T + P ​ 𝐕 \displaystyle\mathbf{M}\mathbf{M}^{T}+P\mathbf{V} (39)

The full variance-covariance matrix of the diagonal 𝐝 \mathbf{d} of 𝐗𝐗 T \mathbf{X}\mathbf{X}^{T} is

Var ⁡ ( diag ⁡ ( 𝐗𝐗 T ) ) = 4 ​ 𝐌𝐌 T ∘ 𝐕 + 2 ​ P ​ ( 𝐕 ∘ 𝐕 ) . \rm Var(\rm diag(\mathbf{XX}^{T}))=4\mathbf{M}\mathbf{M}^{T}\circ\mathbf{V}+2P(\mathbf{V}\circ\mathbf{V}). (40)

Finally, we need to generalize these results to a situation, where the columns of 𝐗 \mathbf{X} are not independent, but have element-wise covariance of 𝚺 \bm{\Sigma} . Thus, we are interested in the joint distribution of the elements of the quadratic form 𝐗𝐗 T \mathbf{X}\mathbf{X}^{T} , where 𝐗 \mathbf{X} has matrix normal distribution 𝐗 ∼ ℳ ​ 𝒩 ​ ( 𝐌 , 𝐕 , 𝚺 ) \mathbf{X}\sim\mathcal{MN}(\mathbf{M},\mathbf{V},\bm{\Sigma}) .

E ⁡ ( 𝐗𝐗 T ) \displaystyle\rm E(\mathbf{X}\mathbf{X}^{T}) = \displaystyle= 𝐌𝐌 T + tr ⁡ ( 𝚺 ) ​ 𝐕 \displaystyle\mathbf{M}\mathbf{M}^{T}+\rm tr(\bm{\Sigma})\mathbf{V} (41) Var ⁡ ( diag ⁡ ( 𝐗𝐗 T ) ) \displaystyle\rm Var(\rm diag(\mathbf{XX}^{T})) = \displaystyle= 4 ​ 𝐌 ​ 𝚺 ​ 𝐌 T ∘ 𝐕 + 2 ​ t ​ r ​ ( 𝚺 ​ 𝚺 ) ​ ( 𝐕 ∘ 𝐕 ) \displaystyle 4\mathbf{M}\bm{\Sigma}\mathbf{M}^{T}\circ\mathbf{V}+2\rm tr(\bm{\Sigma\Sigma})(\mathbf{V}\circ\mathbf{V}) (42)

From this result, we can obtain Eq. 27 by considering that the mean of 𝐂 ​ 𝐁 ^ m \mathbf{C}\hat{\mathbf{B}}_{m} across partitions has variance 𝚵 / M \bm{\Xi}/M .

#### 6.1.3 Averaging across partitions

To derive the variance of the unbiased distances, we need to take into account the averaging of the estimated difference across the M M different crossvalidation folds. While data from different partitions can be assumed to be independent, the inner products across crossvalidation folds are not. This is because the partitions from one crossvalidation fold will be again included in other folds. The two pattern differences that enter the product in Eq. 5 come from a single partition (that is, A = m A={m} ), or from the set of all other partitions (that is, B = M ∖ m B=M\setminus{m} , which we will denote here in short by ∖ m \setminus{m} .

As a shorthand for the covariance between difference estimates i i and j j that are based on the set of partitions A A and B B , we introduce the symbol

Ξ i , j A , B = Cov ⁡ ( 𝜹 ^ i , A , 𝜹 ^ j , B ) . {\Xi}_{i,j}^{A,B}={\rm Cov}\big(\hat{\bm{\delta}}_{i,A},\,\hat{\bm{\delta}}_{j,B}\big). (43)

This is the covariance for each individual voxel. We now exploit the bilinearity of the covariance operator, that is,

Cov ⁡ ( ∑ m x m , ∑ n y n ) = ∑ m ∑ n Cov ⁡ ( x m , y n ) , {\rm Cov}\Big(\sum_{m}{x}_{m},\,\sum_{n}{y}_{n}\Big)=\sum_{m}\sum_{n}{\rm Cov}\big({x}_{m},\,{y}_{n}\big), (44)

to obtain the following general result:

Cov ⁡ ( d ~ i , d ~ j ) = 1 M 2 ​ P 2 ​ ∑ m ∑ n Cov ⁡ ( 𝜹 ^ i , m ​ 𝜹 ^ i , ∖ m T , 𝜹 ^ j , n ​ 𝜹 ^ j , ∖ n T ) = 1 M 2 ​ P 2 ​ ∑ m ∑ n P ​ 𝜹 i ​ 𝚺 P ​ 𝜹 j T ​ ( Ξ i , j ∖ m ∖ n + Ξ i , j ∖ m ​ n + Ξ i , j m ∖ n + Ξ i , j m ​ n ) + tr ​ ( 𝚺 P ​ 𝚺 P ) ​ ( Ξ i , j m ​ n ​ Ξ i , j ∖ m ∖ n + Ξ i , j m ∖ n ​ Ξ i , j ∖ m ​ n ) = 1 P 2 ​ { P ​ 𝜹 i ​ 𝚺 P ​ 𝜹 j T ​ S i , j + tr ⁡ ( 𝚺 P ​ 𝚺 P ) ​ N i , j } \begin{split}{\rm Cov}(\tilde{d}_{i},\,\tilde{d}_{j})=&\frac{1}{{M^{2}P^{2}}}\sum_{m}\sum_{n}{\rm Cov}\big(\hat{\bm{\delta}}_{i,m}\hat{\bm{\delta}}_{i,\setminus m}^{T},\,\hat{\bm{\delta}}_{j,n}\hat{\bm{\delta}}_{j,\setminus n}^{T}\big)\\ =&\frac{1}{M^{2}P^{2}}\sum_{m}\sum_{n}P\bm{\delta}_{i}\bm{\Sigma}_{P}\bm{\delta}_{j}^{T}\Big({\Xi}_{i,j}^{\setminus m\setminus n}+{\Xi}_{i,j}^{\setminus mn}+{\Xi}_{i,j}^{m\setminus n}+{\Xi}_{i,j}^{mn}\Big)\\ &+\text{tr}\big(\bm{\Sigma}_{P}\bm{\Sigma}_{P}\big)\Big({\Xi}_{i,j}^{mn}{\Xi}_{i,j}^{\setminus m\setminus n}+{\Xi}_{i,j}^{m\setminus n}{\Xi}_{i,j}^{\setminus mn}\Big)\ \\ =&\frac{1}{P^{2}}\Big\{P\bm{\delta}_{i}\bm{\Sigma}_{P}\bm{\delta}_{j}^{T}{S}_{i,j}+\rm tr\big(\bm{\Sigma}_{P}\bm{\Sigma}_{P}\big){N}_{i,j}\Big\}\end{split} (45)

where

S i , j = 1 M 2 ​ ∑ m ∑ n { Ξ i , j ∖ m ∖ n + Ξ i , j ∖ m ​ n + Ξ i , j m ∖ n + Ξ i , j m ​ n } {S}_{i,j}=\frac{1}{M^{2}}\sum_{m}\sum_{n}\Big\{{\Xi}_{i,j}^{\setminus m\setminus n}+{\Xi}_{i,j}^{\setminus mn}+{\Xi}_{i,j}^{m\setminus n}+{\Xi}_{i,j}^{mn}\Big\} (46)

and

N i , j = 1 M 2 ​ ∑ m ∑ n { Ξ i , j m ​ n ​ Ξ i , j ∖ m ∖ n + Ξ i , j m ∖ n ​ Ξ i , j ∖ m ​ n } . {N}_{i,j}=\frac{1}{M^{2}}\sum_{m}\sum_{n}\Big\{{\Xi}_{i,j}^{mn}{\Xi}_{i,j}^{\setminus m\setminus n}+{\Xi}_{i,j}^{m\setminus n}{\Xi}_{i,j}^{\setminus mn}\Big\}. (47)

This is the most general expression of the variance of the unbiased distance, which can even be used when the covariance structure of different partitions ( 𝚺 K \bm{\Sigma}_{K} ) differs from each other (see Appendix C).

For the case in which the difference estimates from all M M partitions can be assumed to have the same covariance, that is, 𝚵 i , j ≡ Cov ⁡ ( 𝜹 ^ i , m , 𝜹 ^ j , m ) \mathbf{\Xi}_{i,j}\equiv{\rm Cov}(\hat{\bm{\delta}}_{i,m},\,\hat{\bm{\delta}}_{j,m}) , we can simplify the expression dramatically. In this instance the best estimate of 𝜹 ∖ m \bm{\delta}_{\setminus m} is the average of all partitions except m m :

𝜹 ^ ∖ m = ∑ n ≠ m 𝜹 ^ n / ( M − 1 ) . \hat{\bm{\delta}}_{\setminus m}=\sum_{n\neq m}\hat{\bm{\delta}}_{n}/(M-1). (48)

Accordingly, for n ≠ m n\neq m , we have

𝚵 A , B = { 𝚵 if ​ A = m , B = m 𝟎 if ​ A = m , B = n 𝟎 if ​ A = m , B = ∖ m 𝚵 / ( M − 1 ) if ​ A = ∖ m , B = ∖ m 𝚵 / ( M − 1 ) if ​ A = m , B = ∖ n ( M − 2 ) ​ 𝚵 / ( M − 1 ) 2 if ​ A = ∖ m , B = ∖ n \bm{\Xi}^{A,B}=\left\{\begin{array}[]{lll}\bm{\Xi}&\mbox{if }A=m,\,B=m\\ \mathbf{0}&\mbox{if }A=m,\,B=n\\ \mathbf{0}&\mbox{if }A=m,\,B=\setminus m\\ \bm{\Xi}/(M-1)&\mbox{if }A=\setminus m,\,B=\setminus m\\ \bm{\Xi}/(M-1)&\mbox{if }A=m,\,B=\setminus n\\ (M-2)\bm{\Xi}/(M-1)^{2}&\mbox{if }A=\setminus m,\,B=\setminus n\\ \end{array}\right. (49) .

Substituting the elements of the appropriate representations of 𝚵 A , B \bm{\Xi}^{A,B} into Eq. 46 , 47 and summing up, we have