# Estimating Neural Representation Alignment from Sparsely Sampled Inputs and Features

## Abstract

In both artificial and biological systems, the centered kernel alignment (CKA) has become a widely used tool for quantifying neural representation similarity. While current CKA estimators typically correct for the effects of finite stimuli sampling, the effects of sampling a subset of neurons are overlooked, introducing notable bias in standard experimental scenarios. Here, we provide a theoretical analysis showing how this bias is affected by the representation geometry. We then introduce a novel estimator that corrects for both input and feature sampling. We use our method for evaluating both brain-to-brain and model-to-brain alignments and show that it delivers reliable comparisons even with very sparsely sampled neurons. We perform within-animal and across-animal comparisons on electrophysiological data from visual cortical areas V1, V4, and IT data, and use these as benchmarks to evaluate model-to-brain alignment. We also apply our method to reveal how object representations become progressively disentangled across layers in both biological and artificial systems. These findings underscore the importance of correcting feature-sampling biases in CKA and demonstrate that our bias-corrected estimator provides a more faithful measure of representation alignment. The improved estimates increase our understanding of how neural activity is structured across both biological and artificial systems. Keywords: CKA; Representation Alignment; Disentanglement; Estimation

## Introduction

Over the past decade, the concept of representation similarity has emerged as a powerful framework for comparing complex neural and computational systems ( Kriegeskorte et al., (2008) ; Kriegeskorte & Kievit, (2013) ) . One of the most popular tools is centered kernel alignment (CKA), a measure originally adapted from kernel-based independence metrics but now widely employed across fields such as machine learning, neuroscience, and cognitive science. In machine learning, CKA has become the de facto standard for quantifying how similarly different layers—or even entirely different architectures—encode the same input data ( Kornblith et al., (2019) ) . In neuroscience, CKA has emerged as a core analytical tool to assess whether neural populations, either within or across brain regions and species, produce similar activity patterns in response to identical stimuli ( Yamins & DiCarlo, (2016) ; Schrimpf et al., (2018) ) . Its popularity arises from two advantages: (1) CKA is invariant to orthogonal transformations, making it robust to small perturbations in feature space, and (2) it normalizes for overall variation in activity levels, facilitating meaningful comparisons based on alignment.

Despite these strengths, there is a growing consensus that conventional CKA estimators overlook a critical limitation in many experimental settings, that the sampling of only a subset of “features” as with experimental recordings, can introduce a systematic bias. This is particularly significant in neuroscience, where only a fraction of the total neural population is recorded. Existing CKA estimators often assume that if enough data points (e.g., stimuli or input images) are provided, the measure becomes reliable. Yet this assumption ignores the additional requirement for sufficiently large samples along the “feature” dimension. Consequently, researchers risk drawing misleading conclusions about the alignment of brain regions and neural network layers ( Murphy et al., (2024) ; Cloos et al., (2024) ; Han et al., (2023) ; Sucholutsky et al., (2023) ) .

In this work, we address this pressing concern. First, we show analytically how the geometry of high-dimensional representations contributes to spurious underestimation of similarity when only limited neuronal or model “units” are observed. Second, building on these insights, we introduce a novel estimator designed to remain consistent even with limited feature samples. By systematically correcting for finite-sample effects in both inputs and features, our method offers a more faithful gauge of representation alignment.

We demonstrate the real-world impact of our estimator through analyses of convolutional neural networks and multi-electrode electrophysiological recordings in visual cortices V1, V4, and IT ( Papale et al., (2025) ) . Our new estimator enables more accurate model-to-brain comparisons, revealing alignment trends that are otherwise corrupted by sampling biases. Beyond model-to-brain alignment, we show that the improved CKA estimator illuminates how object-category representations become increasingly disentangled along the primate ventral stream, similar to the observations in deep neural networks. Taken together, these results show that accounting for feature sampling is indispensable for robust representation analysis. Our work thus not only strengthens the theoretical foundations of CKA but also expands its practical utility for probing neural and computational representations.

## Problem Statement

There are many scenarios in which the accurate quantification of the similarity between representations is of central interest ( Figure 1 ). For example, for a given stimulus set, one may want to compare the representations of the brain and a neural network model. To set a benchmark for brain-to-model similarity, it is also useful to measure brain-to-brain similarities across two individual animals or across trials in one animal. In a separate scenario, one may be interested in comparing the neural representations of two stimuli sets, e.g. two image categories, in a single brain region of one individual animal. In all these setups, CKA can be used to quantify the similarities.

Suppose one measures two separate neural representation matrices, 𝐗 {\mathbf{X}} and 𝐘 {\mathbf{Y}} , where each row corresponds to a distinct stimulus and each column represents a recorded neuron. Here, assume the data has been preprocessed so that each neuron has zero mean across all the recorded stimuli. Typically, in its most intuitive form, the CKA is then computed as

tr ​ ( 𝐗𝐗 ⊤ ​ 𝐘𝐘 ⊤ ) tr ​ ( ( 𝐗𝐗 ⊤ ) 2 ) ​ tr ​ ( ( 𝐘𝐘 ⊤ ) 2 ) \frac{\text{tr}({\mathbf{X}}{\mathbf{X}}^{\top}{\mathbf{Y}}{\mathbf{Y}}^{\top})}{\sqrt{\text{tr}\left(\left({\mathbf{X}}{\mathbf{X}}^{\top}\right)^{2}\right)\text{tr}\left(\left({\mathbf{Y}}{\mathbf{Y}}^{\top}\right)^{2}\right)}} (1) In a hypothetical scenario where one has access to all neurons and stimuli, CKA takes the value 1 1 when the representations are perfectly aligned and 0 0 when they span orthogonal subspaces. The normalization factor in the denominator keeps CKA invariant to the scaling of the representations.

In reality, the representation matrices are actually sampled submatrices of an unobserved larger matrix. Due to practical limitations, we often do not get to observe all neurons in a brain region, and it is not possible to present all possible stimuli from a stimulus set, e.g. natural images. Due to these sampling effects, computing the naive CKA from the measurement matrices is heavily biased. Currently, the estimator derived from Song et al. ((2012)) is widely used in the literature to correct the bias contributed by stimulus sampling. In this paper, we show how to correct the bias in the CKA estimator when both the stimuli and neurons are sampled from the underlying large matrix.

## Our Contribution

We make the following contributions in this paper: • We find that the popular CKA estimator based on Song et al. ((2012)) can take arbitrarily small values under finite sampling of both the stimuli and neurons, even when the underlying representations are perfectly aligned. Through theoretical analysis, we show that the bias increases with the intrinsic dimensionality of the representation.

• To address this issue, we develop a more general CKA estimator that corrects the bias contributed by both the stimulus and neuron sampling. We demonstrate its reliability in both synthetic and neurophysiological data.

• We demonstrate our estimator enables a novel application of CKA that quantifies the representation disentanglement in the brain from real data.

## Definitions

### Centered Kernel Alignment

A neural population ϕ ⁡ ( x ) {\bm{\phi}}(x) defines a representation of stimuli x ∈ 𝒳 x\in\mathcal{X} with a large number of neurons. The associated centered kernel function k ⁡ ( x , x ′ ) k(x,x^{\prime}) measures the similarity of two stimuli in the representation space of ϕ {\bm{\phi}} and is defined as: k ⁡ ( x , x ′ ) ≔ ⟨ ϕ ⁡ ( x ) − 𝔼 ⁡ [ ϕ ] , ϕ ⁡ ( x ′ ) − 𝔼 ⁡ [ ϕ ] ⟩ , \displaystyle k(x,x^{\prime})\coloneqq\Braket{\,{\bm{\phi}}(x)-\mathbb{E}[{\bm{\phi}}]\;,\;{\bm{\phi}}(x^{\prime})-\mathbb{E}[{\bm{\phi}}]\,}, (2) where ⟨ ⋅ , ⋅ ⟩ \braket{\cdot,\cdot} denotes an appropriate inner product. Here, 𝔼 ⁡ [ ϕ ] \mathbb{E}[{\bm{\phi}}] denotes the expectation of ϕ {\bm{\phi}} over the stimuli space 𝒳 \mathcal{X} and is used to center the population activity.

For two distinct neural populations ϕ ( a ) ​ ( x ) {\bm{\phi}}^{(a)}(x) and ϕ ( b ) ​ ( x ) {\bm{\phi}}^{(b)}(x) , we measure their similarity using the Hilbert-Schmidt Independence Criterion (HSIC), a popular metric in machine learning ( Gretton et al., (2005) ) . HSIC compares two populations based on the correlations of their kernel functions k ( a ) k^{(a)} and k ( b ) k^{(b)} , and yields, what we call an ℋ \mathcal{H} -value: ℋ ⁡ ( k ( a ) , k ( b ) ) ≔ 𝔼 x , x ′ ​ [ k ( a ) ​ ( x , x ′ ) ​ k ( b ) ​ ( x , x ′ ) ] . \displaystyle\mathcal{H}(k^{(a)},k^{(b)})\coloneqq\mathbb{E}_{x,x^{\prime}}[k^{(a)}(x,x^{\prime})k^{(b)}(x,x^{\prime})]. (3) A low ℋ \mathcal{H} -value indicates less similarity and is zero if and only if the populations ϕ ( a ) {\bm{\phi}}^{(a)} and ϕ ( b ) {\bm{\phi}}^{(b)} vary independently.

The centered kernel alignment (CKA) is essentially the normalized version of HSIC ( Cortes et al., (2012) ; Kornblith et al., (2019) ) and is defined as CKA ​ ( k ( a ) , k ( b ) ) ≔ ℋ ⁡ ( k ( a ) , k ( b ) ) ℋ ⁡ ( k ( a ) , k ( a ) ) ​ ℋ ​ ( k ( b ) , k ( b ) ) . \displaystyle\text{CKA}(k^{(a)},k^{(b)})\coloneqq\frac{\mathcal{H}(k^{(a)},k^{(b)})}{\sqrt{\mathcal{H}(k^{(a)},k^{(a)})\mathcal{H}(k^{(b)},k^{(b)})}}. (4) CKA is normalized to the interval [ 0 , 1 ] [0,1] and invariant to the overall magnitude of population activations.

### Measurement matrix

Greek indices ( α , β , ⋯ \alpha,\beta,\cdots ) denote neurons; Latin indices ( i , j , ⋯ i,j,\cdots ) denote stimuli, and Latin indices ( a , b , ⋯ a,b,\cdots ) denote distinct populations. We use the letter Q Q for the number of neurons and P P for the number of stimuli. Quantities with a hat ( ^ \;\widehat{\;}\; ) indicate empirical estimates.

For each population ϕ ( a ) {\bm{\phi}}^{(a)} , we observe Q a Q_{a} neurons ϕ α ​ ( x ) {\bm{\phi}}_{\alpha}(x) for α = 1 , ⋯ , Q a \alpha=1,\cdots,Q_{a} . Each neuron is measured on the same set of stimuli { x i } \{x_{i}\} for i = 1 , ⋯ , P i=1,\cdots,P , where P P is the number of stimuli. The corresponding measurement matrix 𝚽 ( a ) ∈ ℝ P × Q a {\bm{\Phi}}^{(a)}\in\mathbb{R}^{P\times Q_{a}} with elements 𝚽 i ​ α ( a ) = ϕ α ( a ) ​ ( x i ) {\bm{\Phi}}^{(a)}_{i\alpha}={\bm{\phi}}^{(a)}_{\alpha}(x_{i}) denotes the response of each neuron to stimulus x i x_{i} .

We define the empirical uncentered kernels (also called Gram matrices) by 𝐊 ( a ) = 1 Q a ​ 𝚽 ( a ) ​ 𝚽 ( a ) ⊤ , {\mathbf{K}}^{(a)}=\frac{1}{Q_{a}}{\bm{\Phi}}^{(a)}{{\bm{\Phi}}^{(a)}}^{\top}, and the empirical centered kernels defined in Equation 2 by 𝐊 ¯ ( a ) = 𝐇𝐊 ( a ) ​ 𝐇 , 𝐇 = 𝐈 − 1 P ​ 𝟏𝟏 ⊤ , \displaystyle\bar{\mathbf{K}}^{(a)}={\mathbf{H}}{\mathbf{K}}^{(a)}{\mathbf{H}},\quad{\mathbf{H}}={\mathbf{I}}-\frac{1}{P}\mathbf{1}\mathbf{1}^{\top}, (5) where 𝐇 {\mathbf{H}} is the centering matrix.

## Existing CKA estimators

With these definitions, we now construct estimators of ℋ \mathcal{H} -value from finite data. Here, we convey the overall ideas and provide detailed analyses in Appendix A .

### Naive Estimator

Following the definition of ℋ \mathcal{H} -value in Equation 3 , the naive estimator is ℋ ^ 0 ​ ( k ( a ) , k ( b ) ) \displaystyle\widehat{\mathcal{H}}_{0}(k^{(a)},k^{(b)}) = 1 P 2 ​ ∑ i ≠ j 𝐊 ¯ i ​ j ( a ) ​ 𝐊 ¯ j ​ i ( b ) + 1 P 2 ​ ∑ i 𝐊 ¯ i ​ i ( a ) ​ 𝐊 ¯ i ​ i ( b ) , \displaystyle=\frac{1}{P^{2}}\sum_{i\neq j}\bar{\mathbf{K}}^{(a)}_{ij}\bar{\mathbf{K}}^{(b)}_{ji}+\frac{1}{P^{2}}\sum_{i}\bar{\mathbf{K}}^{(a)}_{ii}\bar{\mathbf{K}}^{(b)}_{ii}, (6) which approaches the true ℋ \mathcal{H} -value as P → ∞ P\to\infty . In this limit, the second term falls as 𝒪 ⁡ ( 1 / P ) \mathcal{O}(1/P) and hence does not contribute to the true ℋ \mathcal{H} -value, demonstrating that the naive estimator ℋ ^ 0 \widehat{\mathcal{H}}_{0} for finite P P is biased 1 1 1 The centering operation is another source of bias for ℋ ^ 0 \widehat{\mathcal{H}}_{0} , which is discussed in Appendix B .

We define CKA ^ 0 ​ ( k ( a ) , k ( b ) ) \widehat{\text{CKA}}_{0}(k^{(a)},k^{(b)}) to be the naive estimator of CKA obtained by replacing ℋ ^ 0 \widehat{\mathcal{H}}_{0} in Equation 4 . Observe that this naive estimator is biased due to the biases of ℋ ^ 0 \widehat{\mathcal{H}}_{0} which appears both in numerator and denominator.

### Stimulus-Corrected Estimator

The cause of the bias in the Equation 6 was the term where two independent stimuli indices ( i , j ) (i,j) coincided. In general, unbiased estimates are obtained by averaging over all indices where each stimulus is used at most once ( Hoeffding, (1948) ) .

An unbiased estimator for ℋ \mathcal{H} -value correcting for these terms was developed by Song et al. ((2012)) and is given by ℋ ^ S ​ ( k ( a ) , k ( b ) ) ≔ 1 P ⁡ ( P − 3 ) \displaystyle\widehat{\mathcal{H}}_{S}(k^{(a)},k^{(b)})\coloneqq\frac{1}{P\left(P-3\right)} (7) × [ tr ​ ( 𝐊 ~ ( a ) ​ 𝐊 ~ ( b ) ) − 2 ​ 𝟏 ⊤ ​ 𝐊 ~ ( a ) ​ 𝐊 ~ ( b ) ​ 𝟏 P − 2 + 𝟏 ⊤ ​ 𝐊 ~ ( a ) ​ 𝟏𝟏 ⊤ ​ 𝐊 ~ ( b ) ​ 𝟏 ( P − 1 ) ​ ( P − 2 ) ] , \displaystyle\times\left[\text{tr}\left(\tilde{\mathbf{K}}^{(a)}\tilde{\mathbf{K}}^{(b)}\right)-2\frac{\mathbf{1}^{\top}\tilde{\mathbf{K}}^{(a)}\tilde{\mathbf{K}}^{(b)}\mathbf{1}}{P-2}+\frac{\mathbf{1}^{\top}\tilde{\mathbf{K}}^{(a)}\mathbf{1}\mathbf{1}^{\top}\tilde{\mathbf{K}}^{(b)}\mathbf{1}}{\left(P-1\right)\left(P-2\right)}\right], where 𝐊 ~ i ​ j = 𝐊 i ​ j ​ ( 1 − δ i ​ j ) \tilde{\mathbf{K}}_{ij}={\mathbf{K}}_{ij}(1-\delta_{ij}) is the uncentered kernel with diagonals removed.

The sample corrected estimator for CKA, denoted by CKA ^ S ​ ( k ( a ) , k ( b ) ) \widehat{\text{CKA}}_{S}(k^{(a)},k^{(b)}) , is obtained by replacing ℋ ^ S \widehat{\mathcal{H}}_{S} in Equation 4 . This is the current version of CKA that is used both in deep learning ( Nguyen et al., (2021) ; Raghu et al., (2021) ; Davari et al., (2022) ) and neuroscience ( Murphy et al., (2024) ) .

Note that this estimator removes the bias due to coinciding stimuli indices. However, the inputs to ℋ ^ S ​ ( k ( a ) , k ( b ) ) \widehat{\mathcal{H}}_{S}(k^{(a)},k^{(b)}) are kernels which involve a sum over neurons. Therefore, if its inputs have shared neurons as in the case of the denominator of CKA, a similar bias discussed in Equation 6 occurs but for neuron indices. Therefore, assuming that the populations ϕ ( a ) {\bm{\phi}}^{(a)} and ϕ ( b ) {\bm{\phi}}^{(b)} have independent neurons, the numerator of CKA ^ S ​ ( k ( a ) , k ( b ) ) \widehat{\text{CKA}}_{S}(k^{(a)},k^{(b)}) remains unbiased, but its denominator is biased due to finite neuron sampling effects (See Section B.2 , Equation S7 ).

### Bias of the Existing Estimator

We find that the representation geometry affects the bias in the widely-used stimulus-corrected estimator CKA ^ S \widehat{\text{CKA}}_{S} . Assuming the variance of the ℋ ^ S \widehat{\mathcal{H}}_{S} estimates are negligible, and the activation variance is normalized neuron-wise, the bias of CKA ^ S \widehat{\text{CKA}}_{S} can be approximated as 𝔼 ⁡ [ CKA ^ S ​ ( k ( a ) , k ( b ) ) ] − CKA ​ ( k ( a ) , k ( b ) ) ≈ [ 1 ( 1 + γ a − 1 Q a ) ​ ( 1 + γ b − 1 Q b ) − 1 ] ​ CKA ​ ( k ( a ) , k ( b ) ) \mathbb{E}\left[\widehat{\text{CKA}}_{S}(k^{(a)},k^{(b)})\right]-\text{CKA}(k^{(a)},k^{(b)})\approx\\ \left[\frac{1}{\sqrt{\left(1+\frac{\gamma_{a}-1}{Q_{a}}\right)\left(1+\frac{\gamma_{b}-1}{Q_{b}}\right)}}-1\right]\text{CKA}(k^{(a)},k^{(b)}) (8) where γ a \gamma_{a} and γ b \gamma_{b} are the intrinsic dimensionalities of the two underlying representations, quantified by the participation ratio of the eigenvalues { λ i } \{\lambda_{i}\} of 𝐊 ¯ \bar{\mathbf{K}} in the infinite data limit: γ = ( ∑ i λ i ) 2 ∑ i λ i 2 . \gamma=\frac{\left(\sum_{i}\lambda_{i}\right)^{2}}{\sum_{i}\lambda_{i}^{2}}. From Equation 8 , we make the following observations:

• CKA ^ S \widehat{\text{CKA}}_{S} always underestimates the true CKA on average.

• The scale by which CKA ^ S \widehat{\text{CKA}}_{S} underestimates is independent of the alignment but only dependent on the intrinsic dimensionalities and neuron sample sizes.

• Having larger intrinsic dimensionalities contributes to the greater underestimation of CKA and requires more neuron samples to mitigate the bias.

Therefore, if intrinsic dimensionality is much greater than the number of feature samples, γ ≫ Q \gamma\gg Q , CKA ^ S \widehat{\text{CKA}}_{S} can take an arbitrarily small value even when two representations are perfectly aligned. On the other hand, our estimator is not affected by this issue. For more details, including the bias of CKA ^ 0 \widehat{\text{CKA}}_{0} , see Section B.2 .

## Stimulus-Neuron-Corrected Estimator

Unlike the previous two estimators, an unbiased estimator correcting for both finite stimulus and neuron sampling must be a function of populations ϕ ( a ) {\bm{\phi}}^{(a)} rather than their kernels.

In this paper, we develop a novel estimator ℋ ^ C \widehat{\mathcal{H}}_{C} that corrects for finite neuron sampling effects when two populations have correlated, e.g., identical neurons. We leave its derivation to Appendix A and report the final estimator here ℋ ^ C ​ ( ϕ ( a ) , ϕ ( b ) ) = 1 P 3 ​ ( P − 3 ) ​ Q a ​ ( Q b − δ a ​ b ) × ∑ i ​ j ​ l ​ m [ c i ​ j ​ i ​ j − 2 ​ P P − 2 ( c i ​ j ​ j ​ l − c i ​ i ​ i ​ j P − c i ​ j ​ j ​ j P + c i ​ i ​ i ​ i 2 ​ P ) + P 2 ( P − 1 ) ​ ( P − 2 ) ( c i ​ j ​ l ​ m − c i ​ i ​ j ​ l P − c j ​ l ​ i ​ i P + c i ​ i ​ j ​ j P 2 ) ] \widehat{\mathcal{H}}_{C}({\bm{\phi}}^{(a)},{\bm{\phi}}^{(b)})=\frac{1}{P^{3}\left(P-3\right)Q_{a}\left(Q_{b}-\delta_{ab}\right)}\\ \times\sum_{ijlm}\left[c_{ijij}-\frac{2P}{P-2}\left(c_{ijjl}-\frac{c_{iiij}}{P}-\frac{c_{ijjj}}{P}+\frac{c_{iiii}}{2P}\right)\right.\\ \left.+\frac{P^{2}}{\left(P-1\right)\left(P-2\right)}\left(c_{ijlm}-\frac{c_{iijl}}{P}-\frac{c_{jlii}}{P}+\frac{c_{iijj}}{P^{2}}\right)\right] (9)

where we introduced the tensor c i ​ j ​ l ​ m c_{ijlm} defined as c i ​ j ​ l ​ m ≔ ∑ α ​ β 𝚽 i ​ α ( a ) ​ 𝚽 j ​ α ( a ) \displaystyle c_{ijlm}\coloneqq\sum_{\alpha\beta}{\bm{\Phi}}_{i\alpha}^{(a)}{\bm{\Phi}}_{j\alpha}^{(a)} 𝚽 l ​ β ( b ) ​ 𝚽 m ​ β ( b ) − δ a ​ b ​ ∑ α 𝚽 i ​ α ( a ) ​ 𝚽 j ​ α ( a ) ​ 𝚽 l ​ α ( b ) ​ 𝚽 m ​ α ( b ) . \displaystyle{\bm{\Phi}}_{l\beta}^{(b)}{\bm{\Phi}}_{m\beta}^{(b)}-\delta_{ab}\sum_{\alpha}{\bm{\Phi}}_{i\alpha}^{(a)}{\bm{\Phi}}_{j\alpha}^{(a)}{\bm{\Phi}}_{l\alpha}^{(b)}{\bm{\Phi}}_{m\alpha}^{(b)}. (10)

The second term in c i ​ j ​ l ​ m c_{ijlm} encodes the finite neuron correction only when two populations have identical neurons ( a = b a=b ). Finally, we define our unbiased CKA estimator as CKA ^ C ​ ( ϕ ( a ) , ϕ ( b ) ) = ℋ ^ C ​ ( ϕ ( a ) , ϕ ( b ) ) ℋ ^ C ​ ( ϕ ( a ) , ϕ ( a ) ) ​ ℋ ^ C ​ ( ϕ ( b ) , ϕ ( b ) ) . \displaystyle\widehat{\text{CKA}}_{C}\left({\bm{\phi}}^{(a)},{\bm{\phi}}^{(b)}\right)=\frac{\widehat{\mathcal{H}}_{C}({\bm{\phi}}^{(a)},{\bm{\phi}}^{(b)})}{\sqrt{\widehat{\mathcal{H}}_{C}({\bm{\phi}}^{(a)},{\bm{\phi}}^{(a)})\widehat{\mathcal{H}}_{C}({\bm{\phi}}^{(b)},{\bm{\phi}}^{(b)})}}. (11) Several remarks are in order: • Our estimator ℋ ^ C \widehat{\mathcal{H}}_{C} reduces to the stimulus-corrected estimator ℋ ^ S \widehat{\mathcal{H}}_{S} of Song et al. ((2012)) when two populations ϕ ( a ) {\bm{\phi}}^{(a)} and ϕ ( b ) {\bm{\phi}}^{(b)} have independent neurons, generalizing previous results.

• When two distinct populations compared, our CKA estimator CKA ^ C \widehat{\text{CKA}}_{C} corrects the bias in the denominator of CKA ^ S \widehat{\text{CKA}}_{S} and their numerators remain identical. However, in certain cases, e.g. trial-to-trial similarity between two measurements of a single population, the numerator of CKA ^ S \widehat{\text{CKA}}_{S} also becomes biased.

• While the estimator looks complicated, it can be efficiently implemented in terms of einsum operations. An implementation is provided in Appendix D .

• Note that while the ℋ \mathcal{H} -estimator we derive is unbiased, the proposed CKA-estimator is still biased due to non-linear operations (multiplication, square-root) involving ℋ ^ C \widehat{\mathcal{H}}_{C} . Nevertheless, our CKA estimator is the least biased among the available estimators. Full analysis on this matter can be found in Section B.3 . Also, as shown in the rest of the paper, we empirically find that the effect of this bias is much smaller compared to the effect of bias due to finite neuron sampling. One can also mitigate this bias by empirically averaging ℋ ^ \widehat{\mathcal{H}} over multiple trials and using that to estimate CKA. If N N estimates of ℋ ^ \widehat{\mathcal{H}} are available from independent trials, then this bias falls like 𝒪 ⁡ ( 1 / N ) \mathcal{O}\left(1/N\right) ( Section B.3 ). We will denote the CKA estimator derived from the empirical average of N N number of ℋ ^ \widehat{\mathcal{H}} by CKA N ^ \widehat{\text{CKA}^{N}} .

## Linear example

Next, we numerically test all three estimators on a simple synthetic dataset with a known CKA value. We consider d − d- dimensional stimuli which are drawn from the distribution 𝐱 ∼ 𝒩 ⁡ ( 0 , 𝐈 d ) {\mathbf{x}}\sim\mathcal{N}(0,{\mathbf{I}}_{d}) . We define a linear population of the form ϕ ( a ) ​ ( 𝐱 ) = 𝐱 ⊤ ​ 𝐰 ( a ) {\bm{\phi}}^{(a)}({\mathbf{x}})={\mathbf{x}}^{\top}{\mathbf{w}}^{(a)} , where each weight is drawn from the distribution 𝐰 ( a ) ∼ 𝒩 ⁡ ( 0 , 𝚺 a ) {\mathbf{w}}^{(a)}\sim\mathcal{N}(0,{\bm{\Sigma}}_{a}) and corresponds to a single neuron. Similarly, we define the population ϕ ( b ) ​ ( 𝐱 ) = 𝐱 ⊤ ​ 𝐰 ( b ) {\bm{\phi}}^{(b)}({\mathbf{x}})={\mathbf{x}}^{\top}{\mathbf{w}}^{(b)} where 𝐰 ( b ) ∼ 𝒩 ⁡ ( 0 , 𝚺 b ) {\mathbf{w}}^{(b)}\sim\mathcal{N}(0,{\bm{\Sigma}}_{b}) . Since the entire stimuli and population distributions are known, the true CKA in Equation 4 can be evaluated exactly: CKA = tr ⁡ ( 𝚺 a ​ 𝚺 b ) tr ⁡ 𝚺 a 2 ​ tr ⁡ 𝚺 b 2 . \displaystyle\text{CKA}=\frac{\tr(\bSigma_a \bSigma_b)}{\sqrt{\tr{\bm{\Sigma}}_{a}^{2}\tr{\bm{\Sigma}}_{b}^{2}}}. (12) We consider 𝚺 a = 𝚺 b = 𝐈 d {\bm{\Sigma}}_{a}={\bm{\Sigma}}_{b}={\mathbf{I}}_{d} in which case both populations perfectly align with a CKA = 1 \text{CKA}=1 . As indicated by Equation 8 , however, despite the perfect alignment of populations, the stimulus-corrected estimator CKA ^ S \widehat{\text{CKA}}_{S} can be arbitrarily small depending on the number of neurons sampled Q Q , and the dimensionality γ \gamma , which is exactly d d in this example. We numerically confirm these findings in Figure 2 and find that both estimators become highly sensitive when a limited number of neurons are observed. On the other hand, our estimator is able to recover the true CKA even at small neuron samples.

In this example, we observe that CKA N ^ \widehat{\text{CKA}^{N}} ’s and CKA ^ \widehat{\text{CKA}} ’s almost completely overlap in all three estimators ( Figure 2 lighter vs. darker lines), except when Q Q is very small, CKA N ^ C \widehat{\text{CKA}^{N}}_{C} and CKA ^ C \widehat{\text{CKA}}_{C} diverge a bit. This means that the bias correction for ℋ \mathcal{H} contributes more to the bias correction of CKA than correcting the bias from the non-linear operations on the ℋ \mathcal{H} estimates. The next section shows that this pattern is observed in neural data as well.

## Practical applications in neuroscience

In most of the brain recordings, the observed neurons are samples of much larger population. Here, we use electrophysiological data of the three key cortical regions for visual processing, V1, V4, and IT, in the order of processing cascade. It has been shown that the neurons in V1 are simple filters, but the later regions V4 and IT are sensitive to semantic information ( Hubel & Wiesel, (1962) ; Majaj et al., (2015) ; Hung et al., (2005) ; Cadena et al., (2024) ) . Papale et al. ((2025)) presented 20,000 natural images of 1,854 object categories from the THINGS dataset ( Hebart et al., (2019) ) to two monkeys and recorded neural responses with multiple electrode arrays over V1, V4, and IT, totaling 1046 electrodes in one monkey and 960 in the other. The measurement value is the average spiking voltage levels of neurons adjacent to a given electrode averaged over a small time window immediately following an image presentation. We view the electrodes and images as rows and columns, respectively.

### Brain-to-Brain alignment

Here use CKA to benchmark the similarity of a given brain region across individual animals (see Figure 1 a). Such brain-to-brain comparisons provide an essential reference for evaluating model-to-brain alignments. In our analysis, we first compare two disjoint sets of electrodes from the same brain region within a single animal. The ground truth in this within-region comparison is a perfect alignment (CKA = 1) on average. However, as shown in Figure 3 a, both the naive estimator CKA ^ 0 \widehat{\text{CKA}}_{0} and the stimulus-corrected estimator CKA ^ S \widehat{\text{CKA}}_{S} significantly underestimate the true similarity, with their estimates strongly dependent on the number of neurons sampled. In contrast, our proposed estimator CKA ^ C \widehat{\text{CKA}}_{C} closely approximates the true value even with very limited neuron sampling, thereby highlighting the risk that conventional CKA estimates may fail to detect even perfect alignment when undersampling is present.

We further extend the analysis by measuring the similarity of each brain region across animals. As depicted in Figure 3 b, CKA ^ C \widehat{\text{CKA}}_{C} yields reliable estimates over a range of neuron sampling sizes, whereas the conventional estimators remain highly sensitive to the sample size. Our results show that V1 representations are relatively conserved across individuals (with CKA around 0.6), IT representations exhibit lower similarity (around 0.3), and V4 lies in between (approximately 0.45). The similarity decreases over the layers of visual processing.

### Brain-to-Model alignment

The alignment between brain representations and artificial neural network layers is a topic of considerable interest. It has been repeatedly observed in the literature that the early visual regions, e.g., V1, have representations that are strikingly similar to the early layers of CNNs, and the later regions, e.g., IT, are aligned with the deeper layers ( Yamins et al., (2014) ; Schrimpf et al., (2018) ; Nonaka et al., (2021) ) . As depicted in Figure 1 b, here we use the CKA estimators to compare the representations between brain regions (V1, V4, and IT) and layers of convolutional neural networks (CNNs), namely ResNet18 and CorNet-s ( He et al., (2016) ; Kubilius et al., (2018) ) . We then see how these CKA estimates compare against the brain-to-brain CKA ^ C \widehat{\text{CKA}}_{C} from the previous section.

We make the same observation of the V1-early layer alignment and IT-late layer alignment from all of the estimators. In the absolute scale, the V1-early layer alignment is generally higher, but the IT-late layer alignment is generally lower ( Figure 4 left three plots). However, CKA ^ C \widehat{\text{CKA}}_{C} estimates that the alignment between animals also decreases by a similar factor. This indicates that in the relative scale, the model-to-brain alignments for all three regions are similarly close to the animal-to-animal alignments, an insight that cannot be reliably found with the other estimators.

We then test the sensitivity of the estimators to the number of electrodes. In Figure 4 (right three plots), we pick the best alignment layer for each brain region and see how this alignment value changes with the number of electrodes. We observe that all estimators are similar for comparison with V1, meaning that there is only a small bias in CKA ^ 0 \widehat{\text{CKA}}_{0} and CKA ^ S \widehat{\text{CKA}}_{S} . However, for the comparison with V4 and IT, the gap between the estimators is generally large, meaning we need our estimator CKA ^ C \widehat{\text{CKA}}_{C} for a reliable measurement. The CKA ^ 0 \widehat{\text{CKA}}_{0} and CKA ^ S \widehat{\text{CKA}}_{S} values are not converged even when all available electrodes are used.

### Object disentanglement

Beyond inter-region comparisons, our estimator enables a novel application of CKA: quantifying the disentanglement of semantic information across object categories in the brain. Similar to the representation disentanglement observed in deeper layers of CNNs, the IT cortical region is believed to encode high-level semantic features that distinguish between objects ( Yamins et al., (2014) ) . In this analysis, we estimate the CKA between pairs of object categories within a single brain region, where the columns of the measurement matrix correspond to stimulus images and the rows correspond to the shared neurons (as illustrated in Figure 1 b). A lower CKA in this context suggests greater semantic disentanglement.

We observe that the estimated CKA between object categories generally decreases over V1, V4, and IT, indicating a gradual semantic disentanglement over these regions ( Figure 5 ). Interestingly, natural object pairs get separated faster than artificial object pairs. This pattern is more strongly observed in CKA ^ C \widehat{\text{CKA}}_{C} than in CKA ^ S \widehat{\text{CKA}}_{S} . Interestingly, there are many cases where CKA ^ S \widehat{\text{CKA}}_{S} indicates high disentanglement (small CKA) but ours CKA ^ C \widehat{\text{CKA}}_{C} indicates low disentanglement (large CKA). This is observed in V1 for all pairwise comparisons: Our estimator consistently estimates CKA value near 1, suggesting V1 does not encode semantic information, whereas the heavily biased CKA ^ S \widehat{\text{CKA}}_{S} spuriously suggests otherwise. Also, in IT, CKA ^ C \widehat{\text{CKA}}_{C} is near 1 for artificial objects pairs, but the bias in CKA ^ S \widehat{\text{CKA}}_{S} spuriously suggests the objects are disentangled.

## Trial-to-Trial Similarity

As our final application, we evaluate the trial-to-trial similarity of neural recordings from the same brain region on the same stimuli set as depicted in Figure 1 c. We use the neural recordings from Papale et al. ((2025)) of visual cortical areas V1, V2, and IT over 30 30 trials on 100 100 images. This case differs from the previous analyses since a single population ϕ ( a ) {\bm{\phi}}^{(a)} is compared across two trials ( ϕ 1 ( a ) {\bm{\phi}}^{(a)}_{1} and ϕ 2 ( a ) {\bm{\phi}}^{(a)}_{2} ). In this case, CKA ^ S \widehat{\text{CKA}}_{S} becomes biased in its numerator as well as its denominator, which our estimator accounts for correctly.

In Figure 6 , we calculate the CKA between all pair-wise single trial measurements for each brain region and report the mean of each estimator. As before, our estimator yields consistent estimates across different neuron samplings Q Q , while the others remain highly biased. Furthermore, we observe that trial-to-trial similarities for regions V1 and V4 are significantly higher than the region IT. This indicates that the trial-to-trial variability in IT is larger than in earlier visual regions.

## Discussion

Murphy et al. ((2024)) also highlights the sensitivity of the naive CKA estimator to the number of features, with a focus on the ratio of the number of inputs and features. They motivate the problem by showing that the naive CKA ( Equation 1 ) between two independent random matrices X X and Y Y (their entries are i.i.d. standard normal) takes a non-zero value which depends on the ratio of the number of stimuli and neurons: 1 ( 1 + P Q a ) ​ ( 1 + P Q b ) \frac{1}{\sqrt{\left(1+\frac{P}{Q_{a}}\right)\left(1+\frac{P}{Q_{b}}\right)}} (13) in the limit of large P P , Q a Q_{a} , and Q b Q_{b} . They observe that CKA ^ S \widehat{\text{CKA}}_{S} takes value 0 0 in the same setup, resolving this issue of having spurious similarities between independent random matrices. While this observation is correct, Murphy et al. ((2024)) interprets this result as CKA ^ S \widehat{\text{CKA}}_{S} resolving the issue of the naive CKA being sensitive to Q a Q_{a} and Q b Q_{b} (denoted by P 1 P_{1} and P 2 P_{2} in their paper), since CKA ^ S \widehat{\text{CKA}}_{S} value is 0 0 regardless of Q a Q_{a} and Q b Q_{b} . Extrapolating this interpretation, they apply CKA ^ S \widehat{\text{CKA}}_{S} to real neural data in an attempt to resolve the issue of CKA sensitivity to the number of neurons.

However, here we explain that this interpretation is only valid when the true CKA is 0, and therefore cannot be generalized to other cases. CKA ^ S \widehat{\text{CKA}}_{S} returns 0 0 in the random matrix setup, not because it corrects neuron sampling, but because it corrects the stimulus sampling. We can see that by taking the number of stimuli to infinity ( P → ∞ P\to\infty ) in Equation 13 , which makes the CKA estimate approach 0 0 . In practical scenarios, where the measurement matrices are not random, CKA ^ S \widehat{\text{CKA}}_{S} is still sensitive to the number of features as we have shown in Equation 8 . CKA ^ C \widehat{\text{CKA}}_{C} presented in this paper resolves this problem.

Finally, our estimator for the ℋ \mathcal{H} -value can be used for de-biasing other HSIC-based CKA measures, such as angular-CKA introduced in Williams et al. ((2021)) and Representational Similarity Analysis (RSA) ( Kriegeskorte et al., (2008) ) , since these measures also ignore the bias coming from finite neuron sampling.

## Conclusion

We have addressed a key limitation in applying CKA to scenarios where only a subset of features—such as neurons or model units—is observable. While CKA is widely used to compare representations in machine learning and neuroscience, existing estimators systematically underestimate similarity when features are undersampled, a common issue in neural data collection or reduced-dimension network analyses.

Theoretically, we showed how the stimulus-corrected ℋ \mathcal{H} estimator still fails to yield unbiased CKA under partial column sampling. We derived a bias-corrected estimator that handles both input and feature sampling, more accurately recovering true representation similarity. Empirically, our method performed well on both synthetic data and electrophysiological recordings from the ventral visual stream (V1, V4, and IT). Conventional approaches underestimated alignment, particularly for high-dimensional or semantically complex representations, whereas our estimator revealed consistent neural similarity patterns and clearer evidence of object-category disentanglement.

## References

Cadena et al. ((2024)) Cadena, S.A., Willeke, K.F., Restivo, K., Denfield, G., Sinz, F.H., Bethge, M.Ecker, A.S. (2024). Diverse task-driven modeling of macaque v4 reveals functional specialization towards semantic tasks. PLOS Computational Biology 20 5 e1012056.

Cloos et al. ((2024)) Cloos, N., Siegel, M., Brincat, S.L., Miller, E.K. & Cueva, C.J. (2024). Differentiable optimization of similarity scores between models and brains. In Iclr 2024 workshop on representational alignment.

Cortes et al. ((2012)) Cortes, C., Mohri, M. & Rostamizadeh, A. (2012). Algorithms for learning kernels based on centered alignment. The Journal of Machine Learning Research 13 795–828.

Davari et al. ((2022)) Davari, M., Horoi, S., Natik, A., Lajoie, G., Wolf, G. & Belilovsky, E. (2022). Reliability of cka as a similarity measure in deep learning. arXiv preprint arXiv:2210.16156 .

Gretton et al. ((2005)) Gretton, A., Bousquet, O., Smola, A. & Schölkopf, B. (2005). Measuring statistical dependence with hilbert-schmidt norms. In International conference on algorithmic learning theory ( 63–77).

Han et al. ((2023)) Han, Y., Poggio, T.A. & Cheung, B. (2023). System identification of neural systems: If we got it right, would we know? In International conference on machine learning ( 12430–12444).

He et al. ((2016)) He, K., Zhang, X., Ren, S. & Sun, J. (2016). Deep residual learning for image recognition. In Proceedings of the ieee conference on computer vision and pattern recognition ( 770–778).

Hebart et al. ((2019)) Hebart, M.N., Dickter, A.H., Kidder, A., Kwok, W.Y., Corriveau, A., Van Wicklin, C. & Baker, C.I. (2019). Things: A database of 1,854 object concepts and more than 26,000 naturalistic object images. PloS one 14 10 e0223792.

Hoeffding ((1948)) Hoeffding, W. (1948). A class of statistics with asymptotically normal distribution. The Annals of Mathematical Statistics 19 3 293–325.

Hubel & Wiesel ((1962)) Hubel, D.H. & Wiesel, T.N. (1962). Receptive fields, binocular interaction and functional architecture in the cat’s visual cortex. The Journal of physiology 160 1 106.

Hung et al. ((2005)) Hung, C.P., Kreiman, G., Poggio, T. & DiCarlo, J.J. (2005). Fast readout of object identity from macaque inferior temporal cortex. Science 310 5749 863–866.

Kornblith et al. ((2019)) Kornblith, S., Norouzi, M., Lee, H. & Hinton, G. (2019). Similarity of neural network representations revisited. In International conference on machine learning ( 3519–3529).

Kriegeskorte & Kievit ((2013)) Kriegeskorte, N. & Kievit, R.A. (2013). Representational geometry: integrating cognition, computation, and the brain. Trends in cognitive sciences 17 8 401–412.

Kriegeskorte et al. ((2008)) Kriegeskorte, N., Mur, M. & Bandettini, P.A. (2008). Representational similarity analysis-connecting the branches of systems neuroscience. Frontiers in systems neuroscience 2 249.

Kubilius et al. ((2018)) Kubilius, J., Schrimpf, M., Nayebi, A., Bear, D., Yamins, D.L. & DiCarlo, J.J. (2018). Cornet: Modeling the neural mechanisms of core object recognition. BioRxiv 408385.

Majaj et al. ((2015)) Majaj, N.J., Hong, H., Solomon, E.A. & DiCarlo, J.J. (2015). Simple learned weighted sums of inferior temporal neuronal firing rates accurately predict human core object recognition performance. Journal of Neuroscience 35 39 13402–13418.

Murphy et al. ((2024)) Murphy, A., Zylberberg, J. & Fyshe, A. (2024). Correcting biased centered kernel alignment measures in biological and artificial neural networks. arXiv preprint arXiv:2405.01012 .

Nguyen et al. ((2021)) Nguyen, T., Raghu, M. & Kornblith, S. (2021). Do wide and deep networks learn the same things? uncovering how neural network representations vary with width and depth. In International conference on learning representations. https://openreview.net/forum?id=KJNcAkY8tY4

Nonaka et al. ((2021)) Nonaka, S., Majima, K., Aoki, S.C. & Kamitani, Y. (2021). Brain hierarchy score: Which deep neural networks are hierarchically brain-like? IScience 24 9 .

Papale et al. ((2025)) Papale, P., Wang, F., Self, M.W. & Roelfsema, P.R. (2025). An extensive dataset of spiking activity to reveal the syntax of the ventral stream. Neuron .

Raghu et al. ((2021)) Raghu, M., Unterthiner, T., Kornblith, S., Zhang, C. & Dosovitskiy, A. (2021). Do vision transformers see like convolutional neural networks? Advances in neural information processing systems 34 12116–12128.

Schrimpf et al. ((2018)) Schrimpf, M., Kubilius, J., Hong, H., Majaj, N.J., Rajalingham, R., Issa, E.B.others (2018). Brain-score: Which artificial neural network for object recognition is most brain-like? BioRxiv 407007.

Song et al. ((2012)) Song, L., Smola, A., Gretton, A., Bedo, J. & Borgwardt, K. (2012). Feature selection via dependence maximization. Journal of Machine Learning Research 13 5 .

Sucholutsky et al. ((2023)) Sucholutsky, I., Muttenthaler, L., Weller, A., Peng, A., Bobu, A., Kim, B.others (2023). Getting aligned on representational alignment. arXiv preprint arXiv:2310.13018 .

Williams et al. ((2021)) Williams, A.H., Kunz, E., Kornblith, S. & Linderman, S. (2021). Generalized shape metrics on neural representations. Advances in Neural Information Processing Systems 34 4738–4750.

Yamins & DiCarlo ((2016)) Yamins, D.L. & DiCarlo, J.J. (2016). Using goal-driven deep learning models to understand sensory cortex. Nature neuroscience 19 3 356–365.

Yamins et al. ((2014)) Yamins, D.L., Hong, H., Cadieu, C.F., Solomon, E.A., Seibert, D. & DiCarlo, J.J. (2014). Performance-optimized hierarchical models predict neural responses in higher visual cortex. Proceedings of the national academy of sciences 111 23 8619–8624.

## Supplementary Information

## Appendix A Full derivation of the ℋ \mathcal{H} estimators

### A.1 Expression of ℋ \mathcal{H} in terms of the uncentered kernel

Recall that the kernel is defined as k ⁡ ( x , y ) = ⟨ ϕ ⁡ ( x ) , ϕ ⁡ ( y ) ⟩ k(x,y)=\left\langle\phi(x),\phi(y)\right\rangle and the centered kernel is k ′ ​ ( x , y ) = ⟨ ϕ ⁡ ( x ) − 𝔼 x ​ [ ϕ ⁡ ( x ) ] , ϕ ⁡ ( y ) − 𝔼 x ​ [ ϕ ⁡ ( x ) ] ⟩ k^{\prime}(x,y)=\left\langle\phi(x)-\mathbb{E}_{x}\left[\phi(x)\right],\phi(y)-\mathbb{E}_{x}\left[\phi(x)\right]\right\rangle . Then ℋ \mathcal{H} is defined as

ℋ ( k ( a ) , k ( b ) ) = 𝔼 x , y [ k ′ ( x , y ) ( a ) k ′ ( x , y ) ( b ) ] \mathcal{H}(k^{(a)},k^{(b)})=\mathbb{E}_{x,y}\left[k^{\prime}{}^{(a)}(x,y)k^{\prime}{}^{(b)}(x,y)\right]

Note that k ′ k^{\prime} can be written in terms of k k :

k ′ ​ ( x , y ) = ⟨ ϕ ⁡ ( x ) − 𝔼 x ​ [ ϕ ⁡ ( x ) ] , ϕ ⁡ ( y ) − 𝔼 x ​ [ ϕ ⁡ ( x ) ] ⟩ k^{\prime}(x,y)=\left\langle\phi(x)-\mathbb{E}_{x}\left[\phi(x)\right],\phi(y)-\mathbb{E}_{x}\left[\phi(x)\right]\right\rangle

= k ⁡ ( x , y ) − 𝔼 z ​ [ k ⁡ ( x , z ) ] − 𝔼 z ​ [ k ⁡ ( z , y ) ] + 𝔼 z , w ​ [ k ⁡ ( z , w ) ] . =k(x,y)-\mathbb{E}_{z}\left[k(x,z)\right]-\mathbb{E}_{z}\left[k(z,y)\right]+\mathbb{E}_{z,w}\left[k(z,w)\right]. Plugging this into our definition of ℋ \mathcal{H} , we arrive at the expression of ℋ \mathcal{H} explicitly in terms of the original uncentered kernels, k ( a ) k^{(a)} and k ( b ) k^{(b)} :

ℋ ⁡ ( k ( a ) , k ( b ) ) = 𝔼 x , y ​ [ k ( a ) ​ ( x , y ) ​ k ( b ) ​ ( x , y ) ] − 2 ​ 𝔼 x , y , z ​ [ k ( a ) ​ ( x , y ) ​ k ( b ) ​ ( x , z ) ] + 𝔼 x , y ​ [ k ( a ) ​ ( x , y ) ] ​ 𝔼 x , y ​ [ k ( b ) ​ ( x , y ) ] \mathcal{H}(k^{(a)},k^{(b)})=\mathbb{E}_{x,y}\left[k^{(a)}(x,y)k^{(b)}(x,y)\right]-2\mathbb{E}_{x,y,z}\left[k^{(a)}(x,y)k^{(b)}(x,z)\right]+\mathbb{E}_{x,y}\left[k^{(a)}(x,y)\right]\mathbb{E}_{x,y}\left[k^{(b)}(x,y)\right]

We want to derive an estimator where each term has an expected value that is equal to an individual term above.

### A.2 Identifying the source of bias in the naive estimator

We begin our analysis by considering the naive ℋ \mathcal{H} estimator:

ℋ ^ 0 ​ ( Φ ( a ) , Φ ( b ) ) = 1 P 2 ​ tr ​ ( H ​ K ( a ) ​ H ​ K ( b ) ) \widehat{\mathcal{H}}_{0}(\Phi^{(a)},\Phi^{(b)})=\frac{1}{P^{2}}\text{tr}\left(HK^{(a)}HK^{(b)}\right) where H = I − 1 P ​ 11 ⊤ H=I-\frac{1}{P}11^{\top} and K ( a ) = 1 Q ​ Φ ( a ) ​ Φ ( a ) ⊤ K^{(a)}=\frac{1}{Q}\Phi^{(a)}\Phi^{(a)\top} . When we expand this, we get

ℋ ^ 0 ​ ( Φ ( a ) , Φ ( b ) ) = 1 P 2 ​ ∑ i ​ j = 1 P K i ​ j ( a ) ​ K i ​ j ( b ) − 2 P 3 ​ ∑ i ​ j ​ l = 1 P K i ​ j ( a ) ​ K j ​ l ( b ) + 1 P 4 ​ ∑ i ​ j ​ l ​ m = 1 P K i ​ j ( a ) ​ K l ​ m ( b ) . \widehat{\mathcal{H}}_{0}(\Phi^{(a)},\Phi^{(b)})=\frac{1}{P^{2}}\sum_{ij=1}^{P}K_{ij}^{(a)}K_{ij}^{(b)}-\frac{2}{P^{3}}\sum_{ijl=1}^{P}K_{ij}^{(a)}K_{jl}^{(b)}+\frac{1}{P^{4}}\sum_{ijlm=1}^{P}K_{ij}^{(a)}K_{lm}^{(b)}.

Note that K ( a ) K^{(a)} and K ( b ) K^{(b)} are dependent on the stimuli { x i } i = 1 P \left\{x_{i}\right\}_{i=1}^{P} such that, for example K i ​ j ( a ) = k ( a ) ​ ( x i , x j ) K_{ij}^{(a)}=k^{(a)}(x_{i},x_{j}) . Now, let us take the average of ℋ ^ 0 \widehat{\mathcal{H}}_{0} over all possible stimuli sets { x i } i = 1 P \left\{x_{i}\right\}_{i=1}^{P} sampled from 𝒳 \mathcal{X} . We can compute the expected value separately for each term and add them together for the final result. Let us consider the first term:

𝔼 { x i } i = 1 P ​ [ 1 P 2 ​ ∑ i ​ j = 1 P K i ​ j ( a ) ​ K i ​ j ( b ) ] = 1 P 2 ​ 𝔼 { x i } i = 1 P ​ [ ∑ i ​ j = 1 P k ( a ) ​ ( x i , x j ) ​ k ( b ) ​ ( x i , x j ) ] \mathbb{E}_{\left\{x_{i}\right\}_{i=1}^{P}}\left[\frac{1}{P^{2}}\sum_{ij=1}^{P}K_{ij}^{(a)}K_{ij}^{(b)}\right]=\frac{1}{P^{2}}\mathbb{E}_{\left\{x_{i}\right\}_{i=1}^{P}}\left[\sum_{ij=1}^{P}k^{(a)}(x_{i},x_{j})k^{(b)}(x_{i},x_{j})\right]

= 1 P 2 ​ ∑ i ≠ j P 𝔼 x , x ′ ​ [ k ( a ) ​ ( x , x ′ ) ​ k ( b ) ​ ( x , x ′ ) ] + 1 P 2 ​ ∑ i P 𝔼 x ​ [ k ( a ) ​ ( x , x ) ​ k ( b ) ​ ( x , x ) ] =\frac{1}{P^{2}}\sum_{i\neq j}^{P}\mathbb{E}_{x,x^{\prime}}\left[k^{(a)}(x,x^{\prime})k^{(b)}(x,x^{\prime})\right]+\frac{1}{P^{2}}\sum_{i}^{P}\mathbb{E}_{x}\left[k^{(a)}(x,x)k^{(b)}(x,x)\right]

= ( P − 1 P ) ​ 𝔼 x , x ′ ​ [ k ( a ) ​ ( x , x ′ ) ​ k ( b ) ​ ( x , x ′ ) ] + 1 P ​ 𝔼 x ​ [ k ( a ) ​ ( x , x ) ​ k ( b ) ​ ( x , x ) ] =\left(\frac{P-1}{P}\right)\mathbb{E}_{x,x^{\prime}}\left[k^{(a)}(x,x^{\prime})k^{(b)}(x,x^{\prime})\right]+\frac{1}{P}\mathbb{E}_{x}\left[k^{(a)}(x,x)k^{(b)}(x,x)\right]

= 𝔼 x , x ′ ​ [ k ( a ) ​ ( x , x ′ ) ​ k ( b ) ​ ( x , x ′ ) ] + 1 P ​ ( 𝔼 x ​ [ k ( a ) ​ ( x , x ) ​ k ( b ) ​ ( x , x ) ] − 𝔼 x , x ′ ​ [ k ( a ) ​ ( x , x ′ ) ​ k ( b ) ​ ( x , x ′ ) ] ) =\mathbb{E}_{x,x^{\prime}}\left[k^{(a)}(x,x^{\prime})k^{(b)}(x,x^{\prime})\right]+\frac{1}{P}\left(\mathbb{E}_{x}\left[k^{(a)}(x,x)k^{(b)}(x,x)\right]-\mathbb{E}_{x,x^{\prime}}\left[k^{(a)}(x,x^{\prime})k^{(b)}(x,x^{\prime})\right]\right) Notice that in the second equality, we separate the sum into two parts, based on whether the indices overlap or not: i = j i=j vs i ≠ j i\neq j . This is essential since the expected values are different for these two cases. Note that in the second term here is 𝒪 ⁡ ( 1 / P ) \mathcal{O}(1/P) , contributing as bias. We would not have gotten this bias if the term 1 P 2 ​ ∑ i ​ j = 1 P K i ​ j ( a ) ​ K i ​ j ( b ) \frac{1}{P^{2}}\sum_{ij=1}^{P}K_{ij}^{(a)}K_{ij}^{(b)} was instead defined as 1 P ⁡ ( P − 1 ) ​ ∑ i ≠ j P K i ​ j ( a ) ​ K i ​ j ( b ) \frac{1}{P\left(P-1\right)}\sum_{i\neq j}^{P}K_{ij}^{(a)}K_{ij}^{(b)} . In that case, the expected value does not have a bias:

𝔼 { x i } i = 1 P ​ [ 1 P ⁡ ( P − 1 ) ​ ∑ i ≠ j P K i ​ j ( a ) ​ K i ​ j ( b ) ] = 𝔼 x , x ′ ​ [ k ( a ) ​ ( x , x ′ ) ​ k ( b ) ​ ( x , x ′ ) ] . \mathbb{E}_{\left\{x_{i}\right\}_{i=1}^{P}}\left[\frac{1}{P\left(P-1\right)}\sum_{i\neq j}^{P}K_{ij}^{(a)}K_{ij}^{(b)}\right]=\mathbb{E}_{x,x^{\prime}}\left[k^{(a)}(x,x^{\prime})k^{(b)}(x,x^{\prime})\right]. Note that 1 P ⁡ ( P − 1 ) \frac{1}{P\left(P-1\right)} is introduced as a scaling factor since P ⁡ ( P − 1 ) P\left(P-1\right) is the number of summand in ∑ i ≠ j P \sum_{i\neq j}^{P} .

### A.3 Derivation of the stimulus-corrected ℋ \mathcal{H} estimator

Applying the same logic to the rest of the terms in ℋ ^ 0 ​ ( Φ ( a ) , Φ ( b ) ) \widehat{\mathcal{H}}_{0}(\Phi^{(a)},\Phi^{(b)}) , we arrive at an estimator that removes the bias from input sampling:

ℋ ^ S ​ ( Φ ( a ) , Φ ( b ) ) = 1 P ⁡ ( P − 1 ) ​ ∑ i ≠ j P K i ​ j ( a ) ​ K i ​ j ( b ) − 2 P ​ ( P − 1 ) ​ ( P − 2 ) ​ ∑ i ≠ j ≠ l P K i ​ j ( a ) ​ K j ​ l ( b ) + 1 P ⁡ ( P − 1 ) ​ ( P − 2 ) ​ ( P − 3 ) ​ ∑ i ≠ j ≠ l ≠ m P K i ​ j ( a ) ​ K l ​ m ( b ) . \widehat{\mathcal{H}}_{S}(\Phi^{(a)},\Phi^{(b)})=\frac{1}{P\left(P-1\right)}\sum_{i\neq j}^{P}K_{ij}^{(a)}K_{ij}^{(b)}-\frac{2}{P\left(P-1\right)\left(P-2\right)}\sum_{i\neq j\neq l}^{P}K_{ij}^{(a)}K_{jl}^{(b)}+\frac{1}{P\left(P-1\right)\left(P-2\right)\left(P-3\right)}\sum_{i\neq j\neq l\neq m}^{P}K_{ij}^{(a)}K_{lm}^{(b)}. This is equivalent to the estimator by Song et al. Note that these sums over disjoint indices, e.g. ∑ i ≠ j ≠ l P \sum_{i\neq j\neq l}^{P} , are practically difficult to compute, so they are decomposed into a linear combination of regular sums. Denoting K i ​ j ( a ) ​ K l ​ m ( b ) K_{ij}^{(a)}K_{lm}^{(b)} as v i ​ j ​ l ​ m v_{ijlm} , the first term can be rewritten as:

∑ i ≠ j P K i ​ j ( a ) ​ K i ​ j ( b ) = ∑ i ​ j v i ​ j ​ i ​ j − ∑ i v i ​ i ​ i ​ i , \sum_{i\neq j}^{P}K_{ij}^{(a)}K_{ij}^{(b)}=\sum_{ij}v_{ijij}-\sum_{i}v_{iiii}, whereas the second term can be:

∑ i ≠ j ≠ l P K i ​ j ( a ) ​ K j ​ l ( b ) = ∑ i ​ j ​ l v i ​ j ​ j ​ l − ∑ i ​ j v i ​ i ​ i ​ j − ∑ i ​ j v i ​ j ​ j ​ j − ∑ i ​ j v i ​ j ​ i ​ j + ∑ i 2 ​ v i ​ i ​ i ​ i , \sum_{i\neq j\neq l}^{P}K_{ij}^{(a)}K_{jl}^{(b)}=\sum_{ijl}v_{ijjl}-\sum_{ij}v_{iiij}-\sum_{ij}v_{ijjj}-\sum_{ij}v_{ijij}+\sum_{i}2v_{iiii}, and the third term:

∑ i ≠ j ≠ l ≠ m P K i ​ j ( a ) ​ K l ​ m ( b ) = ∑ i ​ j ​ l ​ m v i ​ j ​ l ​ m − ∑ i ​ j ​ l ( v i ​ i ​ j ​ l + v j ​ l ​ i ​ i + 4 ​ v i ​ j ​ j ​ l ) + ∑ i ​ j ( v i ​ i ​ j ​ j + 4 ​ v i ​ i ​ i ​ j + 4 ​ v i ​ j ​ j ​ j + 2 ​ v i ​ j ​ i ​ j ) − ∑ i 6 ​ v i ​ i ​ i ​ i \sum_{i\neq j\neq l\neq m}^{P}K_{ij}^{(a)}K_{lm}^{(b)}=\sum_{ijlm}v_{ijlm}-\sum_{ijl}\left(v_{iijl}+v_{jlii}+4v_{ijjl}\right)+\sum_{ij}\left(v_{iijj}+4v_{iiij}+4v_{ijjj}+2v_{ijij}\right)-\sum_{i}6v_{iiii}

Replacing the terms in ℋ ^ S ​ ( Φ ( a ) , Φ ( b ) ) \widehat{\mathcal{H}}_{S}(\Phi^{(a)},\Phi^{(b)}) with these new notations, we arrive at the following expression:

ℋ ^ S ( Φ ( a ) , Φ ( b ) ) = 1 P ⁡ ( P − 3 ) × ( ( ∑ i ​ j v i ​ j ​ i ​ j − ∑ i v i ​ i ​ i ​ i ) − 2 ( P − 2 ) ​ ( ∑ i ​ j ​ l v i ​ j ​ j ​ l − ∑ i ​ j v i ​ i ​ i ​ j − ∑ i ​ j v i ​ j ​ j ​ j + ∑ i v i ​ i ​ i ​ i ) + 1 ( P − 1 ) ​ ( P − 2 ) ​ ( ∑ i ​ j ​ l ​ m v i ​ j ​ l ​ m − ∑ i ​ j ​ l v i ​ i ​ j ​ l − ∑ i ​ j ​ l v j ​ l ​ i ​ i + ∑ i ​ j v i ​ i ​ j ​ j ) ) \widehat{\mathcal{H}}_{S}(\Phi^{(a)},\Phi^{(b)})=\frac{1}{P\left(P-3\right)}\times\\ \left(\left(\sum_{ij}v_{ijij}-\sum_{i}v_{iiii}\right)-\frac{2}{\left(P-2\right)}\left(\sum_{ijl}v_{ijjl}-\sum_{ij}v_{iiij}-\sum_{ij}v_{ijjj}+\sum_{i}v_{iiii}\right)+\frac{1}{\left(P-1\right)\left(P-2\right)}\left(\sum_{ijlm}v_{ijlm}-\sum_{ijl}v_{iijl}-\sum_{ijl}v_{jlii}+\sum_{ij}v_{iijj}\right)\right) (S1)

ℋ ^ S ( Φ ( a ) , Φ ( b ) ) = 1 P 3 ​ ( P − 3 ) × ∑ i ​ j ​ l ​ m ( ( v i ​ j ​ i ​ j − v i ​ i ​ i ​ i P ) − 2 ​ P P − 2 ​ ( v i ​ j ​ j ​ l − v i ​ i ​ i ​ j P − v i ​ j ​ j ​ j P + v i ​ i ​ i ​ i P 2 ) + P 2 ( P − 1 ) ​ ( P − 2 ) ​ ( v i ​ j ​ l ​ m − v i ​ i ​ j ​ l P − v j ​ l ​ i ​ i P + v i ​ i ​ j ​ j P 2 ) ) \widehat{\mathcal{H}}_{S}(\Phi^{(a)},\Phi^{(b)})=\frac{1}{P^{3}\left(P-3\right)}\times\\ \sum_{ijlm}\left(\left(v_{ijij}-\frac{v_{iiii}}{P}\right)-\frac{2P}{P-2}\left(v_{ijjl}-\frac{v_{iiij}}{P}-\frac{v_{ijjj}}{P}+\frac{v_{iiii}}{P^{2}}\right)+\frac{P^{2}}{\left(P-1\right)\left(P-2\right)}\left(v_{ijlm}-\frac{v_{iijl}}{P}-\frac{v_{jlii}}{P}+\frac{v_{iijj}}{P^{2}}\right)\right) (S2)

Suppose K ′ K^{\prime} is a version of K K whose diagonal elements are 0 0 , and v ′ i ​ j ​ l ​ m ≔ K ′ i ​ j ( a ) K ′ l ​ m ( b ) v^{\prime}_{ijlm}\coloneqq K^{\prime}{}_{ij}^{(a)}K^{\prime}{}_{lm}^{(b)} . Then, the above expression simplifies to

ℋ ^ S ​ ( Φ ( a ) , Φ ( b ) ) = 1 P ⁡ ( P − 3 ) ​ ( ∑ i ​ j v i ​ j ​ i ​ j ′ − 2 ( P − 2 ) ​ ∑ i ​ j ​ l v i ​ j ​ j ​ l ′ + 1 ( P − 1 ) ​ ( P − 2 ) ​ ∑ i ​ j ​ l ​ m v i ​ j ​ l ​ m ′ ) \widehat{\mathcal{H}}_{S}(\Phi^{(a)},\Phi^{(b)})=\frac{1}{P\left(P-3\right)}\left(\sum_{ij}v^{\prime}_{ijij}-\frac{2}{\left(P-2\right)}\sum_{ijl}v^{\prime}_{ijjl}+\frac{1}{\left(P-1\right)\left(P-2\right)}\sum_{ijlm}v^{\prime}_{ijlm}\right)

= 1 P ⁡ ( P − 3 ) ( tr ( K ′ K ′ ( a ) ) ( b ) − 2 ( P − 2 ) 1 ⊤ K ′ K ′ ( a ) 1 ( b ) + 1 ( P − 1 ) ​ ( P − 2 ) 1 ⊤ K ′ 11 ⊤ ( a ) K ′ 1 ( b ) ) =\frac{1}{P\left(P-3\right)}\left(\text{tr}\left(K^{\prime}{}^{(a)}K^{\prime}{}^{(b)}\right)-\frac{2}{\left(P-2\right)}1^{\top}K^{\prime}{}^{(a)}K^{\prime}{}^{(b)}1+\frac{1}{\left(P-1\right)\left(P-2\right)}1^{\top}K^{\prime}{}^{(a)}11^{\top}K^{\prime}{}^{(b)}1\right) which is the exact expression found in Song et al. ((2012)) .

### A.4 Derivation of the stimulus-neuron-corrected ℋ \mathcal{H} estimator

Our estimator assumes that the features are also sampled, and the features are correlated or identical in (a) and (b). Let us consider a single term v i ​ j ​ l ​ m v_{ijlm} (redefined here as K i ​ j ( a ) ​ K l ​ m ( a ) K_{ij}^{(a)}K_{lm}^{(a)} reflecting a = b a=b ) and see how the sampling of features contribute to the bias.

v i ​ j ​ l ​ m = K i ​ j ( a ) ​ K l ​ m ( a ) = 1 Q 2 ​ ∑ α ​ β Φ i ​ α ( a ) ​ Φ j ​ α ( a ) ​ Φ l ​ β ( a ) ​ Φ m ​ β ( a ) v_{ijlm}=K_{ij}^{(a)}K_{lm}^{(a)}=\frac{1}{Q^{2}}\sum_{\alpha\beta}\Phi_{i\alpha}^{(a)}\Phi_{j\alpha}^{(a)}\Phi_{l\beta}^{(a)}\Phi_{m\beta}^{(a)} If we were to take the expected value of v i ​ j ​ l ​ m v_{ijlm} over the feature sampling, we need to consider the overlapping of the feature indices α \alpha and β \beta . Therefore, the following should be an unbiased estimator of v i ​ j ​ l ​ m v_{ijlm} :

c i ​ j ​ l ​ m ′ = 1 Q ⁡ ( Q − 1 ) ​ ∑ α ≠ β Φ i ​ α ( a ) ​ Φ j ​ α ( a ) ​ Φ l ​ β ( a ) ​ Φ m ​ β ( a ) = 1 Q ⁡ ( Q − 1 ) ​ ( ∑ α ​ β Φ i ​ α ( a ) ​ Φ j ​ α ( a ) ​ Φ l ​ β ( a ) ​ Φ m ​ β ( a ) − ∑ α Φ i ​ α ( a ) ​ Φ j ​ α ( a ) ​ Φ l ​ α ( a ) ​ Φ m ​ α ( a ) ) . c^{\prime}_{ijlm}=\frac{1}{Q\left(Q-1\right)}\sum_{\alpha\neq\beta}\Phi_{i\alpha}^{(a)}\Phi_{j\alpha}^{(a)}\Phi_{l\beta}^{(a)}\Phi_{m\beta}^{(a)}=\frac{1}{Q\left(Q-1\right)}\left(\sum_{\alpha\beta}\Phi_{i\alpha}^{(a)}\Phi_{j\alpha}^{(a)}\Phi_{l\beta}^{(a)}\Phi_{m\beta}^{(a)}-\sum_{\alpha}\Phi_{i\alpha}^{(a)}\Phi_{j\alpha}^{(a)}\Phi_{l\alpha}^{(a)}\Phi_{m\alpha}^{(a)}\right). In the main text we define c i ​ j ​ l ​ m c_{ijlm} , which is simply c i ​ j ​ l ​ m ′ c^{\prime}_{ijlm} without the 1 Q ⁡ ( Q − 1 ) \frac{1}{Q\left(Q-1\right)} factor, i.e. c i ​ j ​ l ​ m = Q ⁡ ( Q − 1 ) ​ c i ​ j ​ l ​ m ′ c_{ijlm}=Q\left(Q-1\right)c^{\prime}_{ijlm} . Therefore, finally, our estimator can be expressed as

ℋ ^ C ( Φ ( a ) , Φ ( a ) ) = 1 P 3 ​ ( P − 3 ) ​ Q ​ ( Q − 1 ) × ∑ i ​ j ​ l ​ m ( ( c i ​ j ​ i ​ j − c i ​ i ​ i ​ i P ) − 2 ​ P P − 2 ​ ( c i ​ j ​ j ​ l − c i ​ i ​ i ​ j P − c i ​ j ​ j ​ j P + c i ​ i ​ i ​ i P 2 ) + P 2 ( P − 1 ) ​ ( P − 2 ) ​ ( c i ​ j ​ l ​ m − c i ​ i ​ j ​ l P − c j ​ l ​ i ​ i P + c i ​ i ​ j ​ j P 2 ) ) \widehat{\mathcal{H}}_{C}(\Phi^{(a)},\Phi^{(a)})=\frac{1}{P^{3}\left(P-3\right)Q\left(Q-1\right)}\times\\ \sum_{ijlm}\left(\left(c_{ijij}-\frac{c_{iiii}}{P}\right)-\frac{2P}{P-2}\left(c_{ijjl}-\frac{c_{iiij}}{P}-\frac{c_{ijjj}}{P}+\frac{c_{iiii}}{P^{2}}\right)+\frac{P^{2}}{\left(P-1\right)\left(P-2\right)}\left(c_{ijlm}-\frac{c_{iijl}}{P}-\frac{c_{jlii}}{P}+\frac{c_{iijj}}{P^{2}}\right)\right) (S3)

## Appendix B Bias analysis

### B.1 Generative process framework

Here we formalize the problem setup by defining a general formulation of the process that generates the measurement matrices. Consider a pair of measurement matrices Φ ( a ) \Phi_{(a)} and Φ ( b ) \Phi_{(b)} of two systems a a and b b . Let x i ∈ 𝒳 x_{i}\in\mathcal{X} be latent variables for the i i th rows of Φ ( a ) \Phi_{(a)} and Φ ( b ) \Phi_{(b)} , and u α ∈ 𝒰 u_{\alpha}\in\mathcal{U} and v α ∈ 𝒱 v_{\alpha}\in\mathcal{V} be column-latent variables for the α \alpha th column of Φ ( a ) \Phi_{(a)} , and Φ ( b ) \Phi_{(b)} , respectively. Let P P be the number of row-latent variables sampled, and Q a Q_{a} and Q b Q_{b} be the numbers of column-latent variables sampled. Let ϕ a : 𝒳 × 𝒰 → ℝ \phi_{a}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\times\mathcal{U}\to\mathbb{R} be a map that defines the measurement value. Then we assume the entries of Φ ( a ) ∈ ℝ P × Q a \Phi_{(a)}\in\mathbb{R}^{P\times Q_{a}} and Φ ( b ) ∈ ℝ P × Q b \Phi_{(b)}\in\mathbb{R}^{P\times Q_{b}} are defined as Φ i ​ α ( a ) = ϕ a ​ ( x i , u α ) , and ​ Φ i ​ β ( b ) = ϕ b ​ ( x i , v β ) . \Phi_{i\alpha}^{(a)}=\phi_{a}(x_{i},u_{\alpha}),\text{and }\Phi_{i\beta}^{(b)}=\phi_{b}(x_{i},v_{\beta}). In some cases, a latent space might be fully observed in the measurement. For example, N N dimensional layer of neural network activation is given by a feature map ψ ( a ) : 𝒳 → ℝ N \psi_{(a)}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\to\mathbb{R}^{N} , where 𝒳 \mathcal{X} is the input space, and the column latent set 𝒰 \mathcal{U} of cardinality N N is fully observed. Here, 𝒰 \mathcal{U} would be a set of trained neural network weights. The uncentered kernel would be simply k ′ ​ ( x , x ′ ) = 1 N ​ ψ ​ ( x ) ​ ψ ​ ( x ′ ) ⊤ k^{\prime}(x,x^{\prime})=\frac{1}{N}\psi(x)\psi(x^{\prime})^{\top} . If the latent variables cannot be fully observed, we assume there are probability measures over the latent spaces: ρ 𝒳 \rho_{\mathcal{X}} is the probability measures over 𝒳 \mathcal{X} and, and similarly for ρ 𝒰 \rho_{\mathcal{U}} and ρ 𝒱 \rho_{\mathcal{V}} . Assume ϕ a \phi_{a} is square integrable w.r.t. ρ 𝒳 \rho_{\mathcal{X}} and ρ 𝒰 \rho_{\mathcal{U}} , and similarly ϕ b \phi_{b} is also square integrable. The associated uncentered kernels are defined as k a ′ ​ ( x , x ′ ) = ∫ d ​ ρ 𝒰 ​ ( u ) ​ ϕ a ​ ( x , u ) ​ ϕ a ​ ( x ′ , u ) k^{\prime}_{a}(x,x^{\prime})=\int d\rho_{\mathcal{U}}(u)\>\phi_{a}(x,u)\phi_{a}(x^{\prime},u) and k b ′ ​ ( y , y ′ ) = ∫ d ​ ρ 𝒱 ​ ( v ) ​ ϕ b ​ ( y , v ) ​ ϕ b ​ ( y ′ , v ) k^{\prime}_{b}(y,y^{\prime})=\int d\rho_{\mathcal{V}}(v)\>\phi_{b}(y,v)\phi_{b}(y^{\prime},v) . We may also define the associated kernel integral operator:

T k ​ f = ∫ d ​ ρ 𝒳 ​ ( x ) ​ k ′ ​ ( ⋅ , x ) ​ f ​ ( x ) . T_{k}f=\int d\rho_{\mathcal{X}}(x)\,k^{\prime}(\cdot,x)f(x). (S4) Later, the eigenvalues of the operator T k T_{k} will be relevant.

We also define associated covariance kernels: k ~ a ′ ​ ( u , u ′ ) = ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ a ​ ( x , u ) ​ ϕ a ​ ( x , u ′ ) \tilde{k}^{\prime}_{a}(u,u^{\prime})=\int d\rho_{\mathcal{X}}(x)\>\phi_{a}(x,u)\phi_{a}(x,u^{\prime}) and k ~ b ′ ​ ( v , v ′ ) = ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ b ​ ( x , v ) ​ ϕ b ​ ( x , v ′ ) \tilde{k}^{\prime}_{b}(v,v^{\prime})=\int d\rho_{\mathcal{X}}(x)\>\phi_{b}(x,v)\phi_{b}(x,v^{\prime}) .

### B.2 Biases in CKA estimators contributed by the biases in HSIC estimators

Here we derive the analytical expression of the biases of the CKA estimators. Here we assume that each feature is already centered, i.e. ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ a ​ ( x , ⋅ ) = 0 \int d\rho_{\mathcal{X}}(x)\>\phi_{a}(x,\cdot)=0 , and ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ b ​ ( x , ⋅ ) = 0 \int d\rho_{\mathcal{X}}(x)\>\phi_{b}(x,\cdot)=0 . We first compute the expected values of the HSIC estimators. Let 𝒮 𝒳 ≔ { x i } i = 1 P \mathcal{S}_{\mathcal{X}}\coloneqq\left\{x_{i}\right\}_{i=1}^{P} , ℱ 𝒰 ≔ { u i } i = 1 Q a \mathcal{F}_{\mathcal{U}}\coloneqq\left\{u_{i}\right\}_{i=1}^{Q_{a}} , and ℱ 𝒱 ≔ { v i } i = 1 Q b \mathcal{F}_{\mathcal{V}}\coloneqq\left\{v_{i}\right\}_{i=1}^{Q_{b}} be the sets of latent variables sampled independently from ρ 𝒳 , \rho_{\mathcal{X}}, ρ 𝒰 \rho_{\mathcal{U}} , and ρ 𝒱 \rho_{\mathcal{V}} respectively. The naive HSIC estimator is then given by

ℋ 0 ^ ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) = 1 P 2 ​ Q a ​ Q b ​ ∑ i ​ j ∑ α ∑ β ϕ a ​ ( x i , u α ) ​ ϕ a ​ ( x j , u α ) ​ ϕ b ​ ( x i , v β ) ​ ϕ b ​ ( x j , v β ) \widehat{\mathcal{H}_{0}}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}})=\frac{1}{P^{2}Q_{a}Q_{b}}\sum_{ij}\sum_{\alpha}\sum_{\beta}\phi_{a}(x_{i},u_{\alpha})\phi_{a}(x_{j},u_{\alpha})\phi_{b}(x_{i},v_{\beta})\phi_{b}(x_{j},v_{\beta}) if u α u_{\alpha} and v β v_{\beta} are independently sampled for all α \alpha and β \beta combinations, which corresponds to the numerator of CKA. If the column latent variables are identical, and ϕ a = ϕ b \phi_{a}=\phi_{b} , this corresponds to the HSICs in the denominator of naive CKA: ℋ 0 ^ ​ ( ϕ a , 𝒮 𝒳 , ℱ 𝒰 ) = 1 P 2 ​ Q a 2 ​ ∑ i ​ j ∑ α ​ β ϕ a ​ ( x i , u α ) ​ ϕ a ​ ( x j , u α ) ​ ϕ a ​ ( x i , u β ) ​ ϕ a ​ ( x j , u β ) , \widehat{\mathcal{H}_{0}}(\phi_{a},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}})=\frac{1}{P^{2}Q_{a}^{2}}\sum_{ij}\sum_{\alpha\beta}\phi_{a}(x_{i},u_{\alpha})\phi_{a}(x_{j},u_{\alpha})\phi_{a}(x_{i},u_{\beta})\phi_{a}(x_{j},u_{\beta}), ℋ 0 ^ ​ ( ϕ b , 𝒮 𝒳 , ℱ 𝒱 ) = 1 P 2 ​ Q b 2 ​ ∑ i ​ j ∑ α ​ β ϕ b ​ ( x i , v α ) ​ ϕ b ​ ( x j , v α ) ​ ϕ b ​ ( x i , v β ) ​ ϕ b ​ ( x j , v β ) \widehat{\mathcal{H}_{0}}(\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{V}})=\frac{1}{P^{2}Q_{b}^{2}}\sum_{ij}\sum_{\alpha\beta}\phi_{b}(x_{i},v_{\alpha})\phi_{b}(x_{j},v_{\alpha})\phi_{b}(x_{i},v_{\beta})\phi_{b}(x_{j},v_{\beta})

Here, we assume the naive CKA is computed after N N trials, across which the empirical means of the HSIC are computed::

CKA ^ 0 ​ ( ϕ a , ϕ b , ρ 𝒳 , ρ 𝒰 , ρ 𝒱 ) = ∑ l = 1 N ℋ 0 ^ ​ ( ϕ a , ϕ b , 𝒮 𝒳 ( l ) , ℱ 𝒰 ( l ) , ℱ 𝒱 ( l ) ) ∑ l = 1 N ℋ 0 ^ ​ ( ϕ a , 𝒮 𝒳 ( l ) , ℱ 𝒰 ( l ) ) ​ ∑ l = 1 N ℋ 0 ^ ​ ( ϕ b , 𝒮 𝒳 ( l ) , ℱ 𝒱 ( l ) ) \widehat{\text{CKA}}_{0}(\phi_{a},\phi_{b},\rho_{\mathcal{X}},\rho_{\mathcal{U}},\rho_{\mathcal{V}})=\frac{\sum_{l=1}^{N}\widehat{\mathcal{H}_{0}}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}}^{(l)},\mathcal{F}_{\mathcal{U}}^{(l)},\mathcal{F}_{\mathcal{V}}^{(l)})}{\sqrt{\sum_{l=1}^{N}\widehat{\mathcal{H}_{0}}(\phi_{a},\mathcal{S}_{\mathcal{X}}^{(l)},\mathcal{F}_{\mathcal{U}}^{(l)})\sum_{l=1}^{N}\widehat{\mathcal{H}_{0}}(\phi_{b},\mathcal{S}_{\mathcal{X}}^{(l)},\mathcal{F}_{\mathcal{V}}^{(l)})}} We assume N N is large such that all three ℋ ^ \widehat{\mathcal{H}} empirical averages have small variance of order 𝒪 ⁡ ( 1 N ) \mathcal{O}\left(\frac{1}{N}\right) . In this limit, the following approximation is valid:

⟨ CKA ^ 0 ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) ⟩ ≈ ⟨ ℋ 0 ^ ​ ( ϕ a , ϕ b , 𝒮 𝒳 ( l ) , ℱ 𝒰 ( l ) , ℱ 𝒱 ( l ) ) ⟩ ⟨ ℋ 0 ^ ​ ( ϕ a , 𝒮 𝒳 ( l ) , ℱ 𝒰 ( l ) ) ⟩ ​ ⟨ ℋ 0 ^ ​ ( ϕ b , 𝒮 𝒳 ( l ) , ℱ 𝒱 ( l ) ) ⟩ . \left\langle\widehat{\text{CKA}}_{0}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}})\right\rangle\approx\frac{\left\langle\widehat{\mathcal{H}_{0}}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}}^{(l)},\mathcal{F}_{\mathcal{U}}^{(l)},\mathcal{F}_{\mathcal{V}}^{(l)})\right\rangle}{\sqrt{\left\langle\widehat{\mathcal{H}_{0}}(\phi_{a},\mathcal{S}_{\mathcal{X}}^{(l)},\mathcal{F}_{\mathcal{U}}^{(l)})\right\rangle\left\langle\widehat{\mathcal{H}_{0}}(\phi_{b},\mathcal{S}_{\mathcal{X}}^{(l)},\mathcal{F}_{\mathcal{V}}^{(l)})\right\rangle}}. This approximation is also valid when P P , Q a Q_{a} , and Q b Q_{b} are all large. We empirically observe that the expected value of CKA ^ 0 \widehat{\text{CKA}}_{0} obtained via this approximation still accurately predict the expected values of CKA ^ 0 \widehat{\text{CKA}}_{0} empirically computed even with N = 1 N=1 and often even in addition to small P P , Q a Q_{a} , and Q b Q_{b} . This allows us to understand the bias of CKA ^ 0 \widehat{\text{CKA}}_{0} contributed from the biases of ℋ ^ 0 \widehat{\mathcal{H}}_{0} , isolated from the bias contributed by taking products, inverse, and square root of the ℋ ^ 0 \widehat{\mathcal{H}}_{0} estimates, which unnecessarily complicates the analysis.

First, let us take the expected value of the ℋ ^ 0 \widehat{\mathcal{H}}_{0} in the numerator: ⟨ ℋ 0 ^ ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) ⟩ = \displaystyle\left\langle\widehat{\mathcal{H}_{0}}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}})\right\rangle= 1 P 2 ​ Q a ​ Q b ​ ∑ i ​ j ∑ α ∑ β ⟨ ϕ a ​ ( x i , u α ) ​ ϕ a ​ ( x j , u α ) ​ ϕ b ​ ( x i , v β ) ​ ϕ b ​ ( x j , v β ) ⟩ \displaystyle\frac{1}{P^{2}Q_{a}Q_{b}}\sum_{ij}\sum_{\alpha}\sum_{\beta}\left\langle\phi_{a}(x_{i},u_{\alpha})\phi_{a}(x_{j},u_{\alpha})\phi_{b}(x_{i},v_{\beta})\phi_{b}(x_{j},v_{\beta})\right\rangle (S5) = \displaystyle= ⟨ k a ​ ( x , y ) ​ k b ​ ( x , y ) ⟩ x , y + 1 P ​ ( ⟨ k a ​ ( x , x ) ​ k b ​ ( x , x ) ⟩ x − ⟨ k a ​ ( x , y ) ​ k b ​ ( x , y ) ⟩ x , y ) \displaystyle\left\langle k_{a}(x,y)k_{b}(x,y)\right\rangle_{x,y}+\frac{1}{P}\left(\left\langle k_{a}(x,x)k_{b}(x,x)\right\rangle_{x}-\left\langle k_{a}(x,y)k_{b}(x,y)\right\rangle_{x,y}\right) (S6) = \displaystyle= ⟨ k a ​ ( x , y ) ​ k b ​ ( x , y ) ⟩ x , y ​ ( 1 + 1 P ​ ( ζ a ​ b − 1 ) ) \displaystyle\left\langle k_{a}(x,y)k_{b}(x,y)\right\rangle_{x,y}\left(1+\frac{1}{P}\left(\zeta_{ab}-1\right)\right) (S7)

where we have introduced a new variable ζ a ​ b ≔ ⟨ k a ​ ( x , x ) ​ k b ​ ( x , x ) ⟩ x ⟨ k a ​ ( x , y ) ​ k b ​ ( x , y ) ⟩ x , y \zeta_{ab}\coloneqq\frac{\left\langle k_{a}(x,x)k_{b}(x,x)\right\rangle_{x}}{\left\langle k_{a}(x,y)k_{b}(x,y)\right\rangle_{x,y}} that is sensitive to the alignment of the representation. Note that the expected value of the stimulus-corrected estimator ℋ S ^ ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) \widehat{\mathcal{H}_{S}}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}}) can be obtained by taking P → ∞ P\to\infty limit in Equation S7 , which simply ⟨ k a ​ ( x , y ) ​ k b ​ ( x , y ) ⟩ \left\langle k_{a}(x,y)k_{b}(x,y)\right\rangle , which means that the numerator ℋ S ^ \widehat{\mathcal{H}_{S}} is unbiased in this problem setup. However, in an alternative problem setup, such as the trial-to-trial CKA, the numerator ℋ S ^ \widehat{\mathcal{H}_{S}} is still biased, since the neurons are identical across (a) and (b), i.e. u α = v α u_{\alpha}=v_{\alpha} . The bias of the numerator ℋ S ^ \widehat{\mathcal{H}_{S}} in this alternative problem setup is similar to the biases in the denominator ℋ S ^ \widehat{\mathcal{H}_{S}} ’s of both problem setup.

Next, let us take the expected value of one of the ℋ ^ 0 \widehat{\mathcal{H}}_{0} ’s in the denominator:

⟨ ℋ 0 ^ ​ ( ϕ a , 𝒮 𝒳 , ℱ 𝒰 ) ⟩ = 1 P 2 ​ Q a 2 ​ ∑ i ​ j ∑ α ​ β ⟨ ϕ a ​ ( x i , u α ) ​ ϕ a ​ ( x j , u α ) ​ ϕ a ​ ( x i , u β ) ​ ϕ a ​ ( x j , u β ) ⟩ \left\langle\widehat{\mathcal{H}_{0}}(\phi_{a},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}})\right\rangle=\frac{1}{P^{2}Q_{a}^{2}}\sum_{ij}\sum_{\alpha\beta}\left\langle\phi_{a}(x_{i},u_{\alpha})\phi_{a}(x_{j},u_{\alpha})\phi_{a}(x_{i},u_{\beta})\phi_{a}(x_{j},u_{\beta})\right\rangle

= ⟨ k a ​ ( x , y ) 2 ⟩ − 1 P ​ ( ⟨ k a ​ ( x , y ) 2 ⟩ − ⟨ k a ​ ( x , x ) 2 ⟩ ) − 1 Q a ​ ( ⟨ k a ​ ( x , y ) 2 ⟩ − ⟨ k ~ a ​ ( w , w ) 2 ⟩ ) + 1 P ​ Q a ​ ( ⟨ k a ​ ( x , y ) 2 ⟩ − ⟨ k a ​ ( x , x ) 2 ⟩ − ⟨ k ~ a ​ ( w , w ) 2 ⟩ + ⟨ ϕ a ​ ( x , w ) 4 ⟩ ) =\left\langle k_{a}(x,y)^{2}\right\rangle-\frac{1}{P}\left(\left\langle k_{a}(x,y)^{2}\right\rangle-\left\langle k_{a}(x,x)^{2}\right\rangle\right)-\frac{1}{Q_{a}}\left(\left\langle k_{a}(x,y)^{2}\right\rangle-\left\langle\tilde{k}_{a}(w,w)^{2}\right\rangle\right)+\\ \frac{1}{PQ_{a}}\left(\left\langle k_{a}(x,y)^{2}\right\rangle-\left\langle k_{a}(x,x)^{2}\right\rangle-\left\langle\tilde{k}_{a}(w,w)^{2}\right\rangle+\left\langle\phi_{a}(x,w)^{4}\right\rangle\right) (S8)

= ⟨ k a ​ ( x , y ) 2 ⟩ ​ [ 1 + 1 P ​ ( γ a ψ a − 1 ) + 1 Q a ​ ( γ a ψ ~ a − 1 ) − 1 P ​ Q a ​ ( γ a ψ a + γ a ψ ~ a − γ a ρ a − 1 ) ] =\left\langle k_{a}(x,y)^{2}\right\rangle\left[1+\frac{1}{P}\left(\frac{\gamma_{a}}{\psi_{a}}-1\right)+\frac{1}{Q_{a}}\left(\frac{\gamma_{a}}{\tilde{\psi}_{a}}-1\right)-\frac{1}{PQ_{a}}\left(\frac{\gamma_{a}}{\psi_{a}}+\frac{\gamma_{a}}{\tilde{\psi}_{a}}-\frac{\gamma_{a}}{\rho_{a}}-1\right)\right] where we have introduced new variables γ a = ⟨ k a ​ ( x , x ) ⟩ 2 ⟨ k a ​ ( x , y ) 2 ⟩ \gamma_{a}=\frac{\left\langle k_{a}(x,x)\right\rangle^{2}}{\left\langle k_{a}(x,y)^{2}\right\rangle} , ψ a ≔ ⟨ k a ​ ( x , x ) ⟩ 2 ⟨ k a ​ ( x , x ) 2 ⟩ \psi_{a}\coloneqq\frac{\left\langle k_{a}(x,x)\right\rangle^{2}}{\left\langle k_{a}(x,x)^{2}\right\rangle} , ψ ~ a ≔ ⟨ k ~ a ​ ( w , w ) ⟩ 2 ⟨ k ~ a ​ ( w , w ) 2 ⟩ \tilde{\psi}_{a}\coloneqq\frac{\left\langle\tilde{k}_{a}(w,w)\right\rangle^{2}}{\left\langle\tilde{k}_{a}(w,w)^{2}\right\rangle} , and ρ a = ⟨ ϕ a ​ ( x , w ) 2 ⟩ 2 ⟨ ϕ a ​ ( x , w ) 4 ⟩ \rho_{a}=\frac{\left\langle\phi_{a}(x,w)^{2}\right\rangle^{2}}{\left\langle\phi_{a}(x,w)^{4}\right\rangle} . Each of them is participation ratio (PR), i.e. effective/soft count, of some quantities. For discrete quantities, let us define PR as ( ∑ i a i ) 2 ∑ i a i 2 . \frac{\left(\sum_{i}a_{i}\right)^{2}}{\sum_{i}a_{i}^{2}}. It is easy to see that, if N N number of a i a_{i} ’s take value 1 1 and the rest 0 0 , then PR is N N , indicating it is a soft count of non-zero a i a_{i} . The continuous version is

( ∫ d ​ μ ​ ( t ) ​ f ​ ( t ) ) 2 ∫ d ​ μ ​ ( t ) ​ f ​ ( t ) 2 , \frac{\left(\int d\mu(t)\,f(t)\right)^{2}}{\int d\mu(t)\,f(t)^{2}}, which we call PR of f f w.r.t. to μ \mu . With these definitions, we can interpret γ a \gamma_{a} as the PR of the eigenvalues of T k a T_{k_{a}} , i.e. intrinsic dimensionality of T k T_{k} , ψ a \psi_{a} as the PR of ∫ d ​ ρ 𝒰 ​ ( w ) ​ ϕ a ​ ( ⋅ , w ) 2 \int d\rho_{\mathcal{U}}(w)\>\phi_{a}(\cdot,w)^{2} w.r.t. ρ 𝒳 \rho_{\mathcal{X}} , ψ ~ a \tilde{\psi}_{a} as the PR of ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ a ​ ( x , ⋅ ) 2 \int d\rho_{\mathcal{X}}(x)\>\phi_{a}(x,\cdot)^{2} w.r.t. ρ 𝒰 \rho_{\mathcal{U}} , and ρ a \rho_{a} as the PR of ϕ a 2 \phi_{a}^{2} w.r.t. ρ 𝒳 ⊗ ρ 𝒰 \rho_{\mathcal{X}}\otimes\rho_{\mathcal{U}} . Suppose Φ ∞ ( a ) \Phi_{\infty}^{(a)} is the measurement matrix in the limit of infinite stimulus and neuron samples, and 𝐊 ∞ ( a ) \mathbf{K}_{\infty}^{(a)} is the corresponding Gram matrix, equivalent to T k a T_{k_{a}} . Then we can loosely say γ a \gamma_{a} is the intrinsic dimensionality of 𝐊 ∞ \mathbf{K}_{\infty} , ψ a \psi_{a} is the effective number of rows of Φ ∞ ( a ) \Phi_{\infty}^{(a)} with non-zero lengths, ψ a ~ \tilde{\psi_{a}} is the effective number of columns of Φ ∞ ( a ) \Phi_{\infty}^{(a)} with non-zero lengths, and ρ a \rho_{a} is the effective number of non-zero entries in Φ ∞ ( a ) \Phi_{\infty}^{(a)} .

Then, with the expected values of ℋ ^ 0 \widehat{\mathcal{H}}_{0} derived above the expected value of CKA ^ 0 \widehat{\text{CKA}}_{0} can be approximated as

⟨ CKA ^ 0 ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) ⟩ ≈ ( 1 + 1 P ​ ( ζ a ​ b − 1 ) ) ​ CKA ( 1 + γ a ψ a − 1 P + γ a ψ ~ a − 1 Q a − γ a ψ a + γ a ψ ~ a − γ a ρ a − 1 P ​ Q a ) ​ ( 1 + γ b ψ b − 1 P + γ b ψ ~ b − 1 Q b − γ b ψ b + γ b ψ ~ b − γ b ρ b − 1 P ​ Q b ) \left\langle\widehat{\text{CKA}}_{0}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}})\right\rangle\approx\frac{\left(1+\frac{1}{P}\left(\zeta_{ab}-1\right)\right)\text{CKA}}{\sqrt{\left(1+\frac{\frac{\gamma_{a}}{\psi_{a}}-1}{P}+\frac{\frac{\gamma_{a}}{\tilde{\psi}_{a}}-1}{Q_{a}}-\frac{\frac{\gamma_{a}}{\psi_{a}}+\frac{\gamma_{a}}{\tilde{\psi}_{a}}-\frac{\gamma_{a}}{\rho_{a}}-1}{PQ_{a}}\right)\left(1+\frac{\frac{\gamma_{b}}{\psi_{b}}-1}{P}+\frac{\frac{\gamma_{b}}{\tilde{\psi}_{b}}-1}{Q_{b}}-\frac{\frac{\gamma_{b}}{\psi_{b}}+\frac{\gamma_{b}}{\tilde{\psi}_{b}}-\frac{\gamma_{b}}{\rho_{b}}-1}{PQ_{b}}\right)}} (S9)

where CKA ≔ ⟨ k a ​ ( x , y ) ​ k b ​ ( x , y ) ⟩ ⟨ k a ​ ( x , y ) 2 ⟩ ​ ⟨ k b ​ ( x , y ) 2 ⟩ \text{CKA}\coloneqq\frac{\left\langle k_{a}(x,y)k_{b}(x,y)\right\rangle}{\sqrt{\left\langle k_{a}(x,y)^{2}\right\rangle\left\langle k_{b}(x,y)^{2}\right\rangle}} is the true CKA. Taking P → ∞ P\to\infty , we obtain the expected value of CKA ^ S ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) \widehat{\text{CKA}}_{S}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}}) estimator:

⟨ CKA ^ S ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) ⟩ ≈ CKA ( 1 + γ a ψ ~ a − 1 Q a ) ​ ( 1 + γ b ψ ~ b − 1 Q b ) \left\langle\widehat{\text{CKA}}_{S}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}})\right\rangle\approx\frac{\text{CKA}}{\sqrt{\left(1+\frac{\frac{\gamma_{a}}{\tilde{\psi}_{a}}-1}{Q_{a}}\right)\left(1+\frac{\frac{\gamma_{b}}{\tilde{\psi}_{b}}-1}{Q_{b}}\right)}} If ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ a ​ ( x , ⋅ ) 2 \int d\rho_{\mathcal{X}}(x)\>\phi_{a}(x,\cdot)^{2} and ∫ d ​ ρ 𝒳 ​ ( x ) ​ ϕ b ​ ( x , ⋅ ) 2 \int d\rho_{\mathcal{X}}(x)\>\phi_{b}(x,\cdot)^{2} are constant, i.e. the norm of the activation is normalized for each neuron, then ψ ~ a = ψ ~ b = 1 \tilde{\psi}_{a}=\tilde{\psi}_{b}=1 , which gives

⟨ CKA ^ S ​ ( ϕ a , ϕ b , 𝒮 𝒳 , ℱ 𝒰 , ℱ 𝒱 ) ⟩ ≈ CKA ( 1 + γ a − 1 Q a ) ​ ( 1 + γ b − 1 Q b ) \left\langle\widehat{\text{CKA}}_{S}(\phi_{a},\phi_{b},\mathcal{S}_{\mathcal{X}},\mathcal{F}_{\mathcal{U}},\mathcal{F}_{\mathcal{V}})\right\rangle\approx\frac{\text{CKA}}{\sqrt{\left(1+\frac{\gamma_{a}-1}{Q_{a}}\right)\left(1+\frac{\gamma_{b}-1}{Q_{b}}\right)}}

### B.3 Analyzing the overall ℋ ^ \widehat{\mathcal{H}} -based CKA estimators

All CKA estimators based on ℋ ^ \widehat{\mathcal{H}} , including ours, have bias coming from taking the product of ℋ ^ \widehat{\mathcal{H}} ’s that are correlated to each other, taking inverse, and taking the square root. As a simple analog of CKA estimator, consider the following problem. We have fixed quantities X X , A A , and B B , and corresponding potentially biased estimators x x , a a , and b b , respectively. Here, x x , a a , and b b are correlated. We wish to study the properties of the estimator T = x a ​ b , T=\frac{x}{\sqrt{ab}}, as an estimator for θ = X A ​ B . \theta=\frac{X}{\sqrt{AB}}. Here, the fixed quantities X X , A A , and B B correspond to the ground truth ℋ \mathcal{H} values, and x x , a a , and b b correspond to the ℋ ^ \widehat{\mathcal{H}} estimators.

We denote δ x \displaystyle\delta_{x} = E [ x ] − X , δ a = E [ a ] − A , δ b = E [ b ] − B , \displaystyle=E[x]-X,\quad\delta_{a}=E[a]-A,\quad\delta_{b}=E[b]-B, σ x 2 \displaystyle\sigma_{x}^{2} = Var ( x ) , σ a 2 = Var ( a ) , σ b 2 = Var ( b ) , \displaystyle=\text{Var}(x),\quad\sigma_{a}^{2}=\text{Var}(a),\quad\sigma_{b}^{2}=\text{Var}(b), σ x ​ a \displaystyle\sigma_{xa} = Cov ( x , a ) , σ x ​ b = Cov ( x , b ) , σ a ​ b = Cov ( a , b ) . \displaystyle=\text{Cov}(x,a),\quad\sigma_{xb}=\text{Cov}(x,b),\quad\sigma_{ab}=\text{Cov}(a,b). We also define the deviations: Δ ​ x = x − X , Δ ​ a = a − A , Δ ​ b = b − B , \Delta x=x-X,\quad\Delta a=a-A,\quad\Delta b=b-B, with E ⁡ [ Δ ​ x ] = δ x , E ⁡ [ Δ ​ a ] = δ a , E ⁡ [ Δ ​ b ] = δ b . E[\Delta x]=\delta_{x},\quad E[\Delta a]=\delta_{a},\quad E[\Delta b]=\delta_{b}.

We consider the function f ⁡ ( x , a , b ) = x a ​ b , f(x,a,b)=\frac{x}{\sqrt{ab}}, and expand it about the point ( X , A , B ) (X,A,B) to second order:

f ⁡ ( x , a , b ) ≈ f ⁡ ( X , A , B ) + f x ​ Δ ​ x + f a ​ Δ ​ a + f b ​ Δ ​ b + 1 2 ​ [ f x ​ x ​ ( Δ ​ x ) 2 + f a ​ a ​ ( Δ ​ a ) 2 + f b ​ b ​ ( Δ ​ b ) 2 + 2 ​ f x ​ a ​ Δ ​ x ​ Δ ​ a + 2 ​ f x ​ b ​ Δ ​ x ​ Δ ​ b + 2 ​ f a ​ b ​ Δ ​ a ​ Δ ​ b ] f(x,a,b)\approx f(X,A,B)+f_{x}\Delta x+f_{a}\Delta a+f_{b}\Delta b+\frac{1}{2}\left[f_{xx}(\Delta x)^{2}+f_{aa}(\Delta a)^{2}+f_{bb}(\Delta b)^{2}+2f_{xa}\Delta x\Delta a+2f_{xb}\Delta x\Delta b+2f_{ab}\Delta a\Delta b\right] with all derivatives evaluated at ( X , A , B ) (X,A,B) .

The zeroth-order term:

f ⁡ ( X , A , B ) = X A ​ B . f(X,A,B)=\frac{X}{\sqrt{AB}}.

The first order derivatives: f x ​ ( X , A , B ) = 1 A ​ B , f a ​ ( X , A , B ) = − X 2 ​ A ​ A ​ B , f b ​ ( X , A , B ) = − X 2 ​ B ​ A ​ B . f_{x}(X,A,B)=\frac{1}{\sqrt{AB}},\quad f_{a}(X,A,B)=-\frac{X}{2A\sqrt{AB}},\quad f_{b}(X,A,B)=-\frac{X}{2B\sqrt{AB}}.

Next, the second-order derivatives:

f x ​ x ​ ( X , A , B ) = 0 , f x ​ a ​ ( X , A , B ) = − 1 2 ​ A ​ A ​ B , f x ​ b ​ ( X , A , B ) = − 1 2 ​ B ​ A ​ B , f_{xx}(X,A,B)=0,\quad f_{xa}(X,A,B)=-\frac{1}{2A\sqrt{AB}},\quad f_{xb}(X,A,B)=-\frac{1}{2B\sqrt{AB}},\quad

f a ​ a ​ ( X , A , B ) = 3 ​ X 4 ​ A 2 ​ A ​ B , f b ​ b ​ ( X , A , B ) = 3 ​ X 4 ​ B 2 ​ A ​ B , f a ​ b ​ ( X , A , B ) = X 4 ​ A ​ B ​ A ​ B . f_{aa}(X,A,B)=\frac{3X}{4A^{2}\sqrt{AB}},\quad f_{bb}(X,A,B)=\frac{3X}{4B^{2}\sqrt{AB}},\quad f_{ab}(X,A,B)=\frac{X}{4AB\sqrt{AB}}. Plugging them into our expansion, we get

x a ​ b ≈ X A ​ B ​ ( 1 + Δ ​ x X − 1 2 ​ Δ ​ a A − 1 2 ​ Δ ​ b B + 3 8 ​ ( Δ ​ a ) 2 A 2 + 3 8 ​ ( Δ ​ b ) 2 B 2 − 1 2 ​ Δ ​ x ​ Δ ​ a X ​ A − 1 2 ​ Δ ​ x ​ Δ ​ b X ​ B + 1 4 ​ Δ ​ a ​ Δ ​ b A ​ B ) \frac{x}{\sqrt{ab}}\approx\frac{X}{\sqrt{AB}}\left(1+\frac{\Delta x}{X}-\frac{1}{2}\frac{\Delta a}{A}-\frac{1}{2}\frac{\Delta b}{B}+\frac{3}{8}\frac{(\Delta a)^{2}}{A^{2}}+\frac{3}{8}\frac{(\Delta b)^{2}}{B^{2}}-\frac{1}{2}\frac{\Delta x\Delta a}{XA}-\frac{1}{2}\frac{\Delta x\Delta b}{XB}+\frac{1}{4}\frac{\Delta a\Delta b}{AB}\right)

Now let us take the expected value of the above. We use E ⁡ [ Δ ​ x ] \displaystyle E[\Delta x] = δ x , E [ Δ a ] = δ a , E [ Δ b ] = δ b , \displaystyle=\delta_{x},\quad E[\Delta a]=\delta_{a},\quad E[\Delta b]=\delta_{b}, E ⁡ [ ( Δ ​ x ) 2 ] \displaystyle E[(\Delta x)^{2}] = σ x 2 + δ x 2 , E [ ( Δ a ) 2 ] = σ a 2 + δ a 2 , E [ ( Δ b ) 2 ] = σ b 2 + δ b 2 , \displaystyle=\sigma_{x}^{2}+\delta_{x}^{2},\quad E[(\Delta a)^{2}]=\sigma_{a}^{2}+\delta_{a}^{2},\quad E[(\Delta b)^{2}]=\sigma_{b}^{2}+\delta_{b}^{2}, E ⁡ [ Δ ​ x ​ Δ ​ a ] \displaystyle E[\Delta x\Delta a] = σ x ​ a + δ x δ a , E [ Δ x Δ b ] = σ x ​ b + δ x δ b , E [ Δ a Δ b ] = σ a ​ b + δ a δ b . \displaystyle=\sigma_{xa}+\delta_{x}\delta_{a},\quad E[\Delta x\Delta b]=\sigma_{xb}+\delta_{x}\delta_{b},\quad E[\Delta a\,\Delta b]=\sigma_{ab}+\delta_{a}\delta_{b}. With these definitions, we get the following as the expected value up to the second order approximation: E ⁡ [ x a ​ b ] ≈ X A ​ B ​ ( 1 + δ x X − 1 2 ​ δ a A − 1 2 ​ δ b B + 3 8 ​ σ a 2 + δ a 2 A 2 + 3 8 ​ σ b 2 + δ b 2 B 2 − 1 2 ​ σ x ​ a + δ x ​ δ a X ​ A − 1 2 ​ σ x ​ b + δ x ​ δ b X ​ B + 1 4 ​ σ a ​ b + δ a ​ δ b A ​ B ) . E\left[\frac{x}{\sqrt{ab}}\right]\approx\frac{X}{\sqrt{AB}}\left(1+\frac{\delta_{x}}{X}-\frac{1}{2}\frac{\delta_{a}}{A}-\frac{1}{2}\frac{\delta_{b}}{B}+\frac{3}{8}\frac{\sigma_{a}^{2}+\delta_{a}^{2}}{A^{2}}+\frac{3}{8}\frac{\sigma_{b}^{2}+\delta_{b}^{2}}{B^{2}}-\frac{1}{2}\frac{\sigma_{xa}+\delta_{x}\delta_{a}}{XA}-\frac{1}{2}\frac{\sigma_{xb}+\delta_{x}\delta_{b}}{XB}+\frac{1}{4}\frac{\sigma_{ab}+\delta_{a}\delta_{b}}{AB}\right). Therefore, the bias is given by

E ⁡ [ x a ​ b ] − X A ​ B ≈ X A ​ B ​ ( δ x X − 1 2 ​ δ a A − 1 2 ​ δ b B + 3 8 ​ σ a 2 + δ a 2 A 2 + 3 8 ​ σ b 2 + δ b 2 B 2 − 1 2 ​ σ x ​ a + δ x ​ δ a X ​ A − 1 2 ​ σ x ​ b + δ x ​ δ b X ​ B + 1 4 ​ σ a ​ b + δ a ​ δ b A ​ B ) . E\left[\frac{x}{\sqrt{ab}}\right]-\frac{X}{\sqrt{AB}}\approx\frac{X}{\sqrt{AB}}\left(\frac{\delta_{x}}{X}-\frac{1}{2}\frac{\delta_{a}}{A}-\frac{1}{2}\frac{\delta_{b}}{B}+\frac{3}{8}\frac{\sigma_{a}^{2}+\delta_{a}^{2}}{A^{2}}+\frac{3}{8}\frac{\sigma_{b}^{2}+\delta_{b}^{2}}{B^{2}}-\frac{1}{2}\frac{\sigma_{xa}+\delta_{x}\delta_{a}}{XA}-\frac{1}{2}\frac{\sigma_{xb}+\delta_{x}\delta_{b}}{XB}+\frac{1}{4}\frac{\sigma_{ab}+\delta_{a}\delta_{b}}{AB}\right).

Now let us put this back into the perspective of the CKA estimation. In that context, δ x \delta_{x} is the bias of an ℋ \mathcal{H} -estimator for the numerator, whereas δ a \delta_{a} and δ b \delta_{b} are the bias of the ℋ \mathcal{H} -estimators in the numerator of CKA estimator. The variances σ x 2 \sigma_{x}^{2} , σ a 2 \sigma_{a}^{2} , and σ b 2 \sigma_{b}^{2} are the variances of the corresponding ℋ \mathcal{H} -estimators. Again, X X , A A , B B corresponds to the true ℋ \mathcal{H} -values in the true CKA. In our estimator CKA ^ C \widehat{\text{CKA}}_{C} , the ℋ ^ C \widehat{\mathcal{H}}_{C} biases δ x \delta_{x} , δ a \delta_{a} , and δ b \delta_{b} are zero, substantially reducing the bias in CKA ^ C \widehat{\text{CKA}}_{C} . However, in CKA ^ 0 \widehat{\text{CKA}}_{0} and CKA ^ S \widehat{\text{CKA}}_{S} , their ℋ \mathcal{H} estimators have non-zero biases δ x , δ a , δ b > 0 \delta_{x},\delta_{a},\delta_{b}>0 of order 𝒪 ⁡ ( 1 P + 1 Q ) \mathcal{O}\left(\frac{1}{P}+\frac{1}{Q}\right) (except, for some setups, δ x \delta_{x} is zero in CKA ^ S \widehat{\text{CKA}}_{S} ), adding even more bias to the CKA estimations. Note that all ℋ \mathcal{H} estimators have variance ( σ x 2 \sigma_{x}^{2} , σ a 2 \sigma_{a}^{2} , σ b 2 \sigma_{b}^{2} , σ x ​ a \sigma_{xa} , σ x ​ b \sigma_{xb} , and σ a ​ b \sigma_{ab} ) of order 𝒪 ⁡ ( 1 P + 1 Q ) \mathcal{O}\left(\frac{1}{P}+\frac{1}{Q}\right) . Therefore, as long as CKA estimation is built based on the ℋ \mathcal{H} estimators, the overall bias of CKA is of order 𝒪 ⁡ ( 1 P + 1 Q ) \mathcal{O}\left(\frac{1}{P}+\frac{1}{Q}\right) .

In a scenario where we have two systems (a) and (b) and we aim to compare the representations in (a) and (b) for a single input distribution. However, for each trial, one uses a unique set of inputs and observes a unique set of neurons independently drawn from distributions. Then, instead of taking the empirical average of CKA ^ \widehat{\text{CKA}} , one can take the empirical averages of ℋ ^ \widehat{\mathcal{H}} and use that average to compute CKA. If there are N N trials, these empirical averages of ℋ ^ \widehat{\mathcal{H}} have the variance of order 𝒪 ​ ( 1 N ​ ( 1 P + 1 Q ) ) \mathcal{O}\left(\frac{1}{N}\left(\frac{1}{P}+\frac{1}{Q}\right)\right) . This situation corresponds to when N N number of labs use distinct individual animals to study V1 and they use distinct sets of natural images. Alternatively, if one aims to compute trial-to-trial CKA and there are M M trials, then by performing all pairwise comparisons, one can reduce the bias of CKA estimators significantly since N = ( M 2 ) N=\binom{M}{2} . In summary, having multiple trials allows our estimator CKA ^ C \widehat{\text{CKA}}_{C} to have bias of 𝒪 ​ ( 1 N ​ ( 1 P + 1 Q ) ) \mathcal{O}\left(\frac{1}{N}\left(\frac{1}{P}+\frac{1}{Q}\right)\right) , while the other estimators CKA ^ 0 \widehat{\text{CKA}}_{0} and CKA ^ S \widehat{\text{CKA}}_{S} still has bias of 𝒪 ⁡ ( 1 P + 1 Q ) \mathcal{O}\left(\frac{1}{P}+\frac{1}{Q}\right) .

## Appendix C Additional results on neural data

In this section, we present the results on the second monkey subject. Overall, our observations are consistent across the two monkeys.

### C.1 Brain-to-model alignment

The CKA estimates between CNNs and the second monkey is shown in Figure S1 .

### C.2 Object disentanglement

The CKA estimates object disentanglement in the second monkey is shown in Figure S2 .

## Appendix D Code availability

The code for the estimators and generating all figures is available in https://github.com/badooki/CKA/ .

## Instructions for reporting errors

We are continuing to improve HTML versions of papers, and your feedback helps enhance accessibility and mobile support. To report errors in the HTML that will help us improve conversion and rendering, choose any of the methods listed below:

Click the "Report Issue" ( ) button, located in the page header.

Tip: You can select the relevant text first, to include it in your report.

Our team has already identified the following issues . We appreciate your time reviewing and reporting rendering errors we may not have found yet. Your efforts will help us improve the HTML versions for all readers, because disability should not be a barrier to accessing research. Thank you for your continued support in championing open access for all.

Have a free development cycle? Help support accessibility at arXiv! Our collaborators at LaTeXML maintain a list of packages that need conversion , and welcome developer contributions .