Original Article Analysis Department of Biomedical Engineering, City
University of New York, New York, NY 10031, USA<EMAIL_ADDRESS>This work
was supported through grants from NSF (DRL-1660548) and ARL
(W911NF-10-2-0022).

# Correlated Components Analysis – Extracting Reliable Dimensions in
Multivariate Data

Lucas C. Parra  Department of Biomedical Engineering, City University of New
York, New York, NY 10031, USA  Stefan Haufe  Berlin Center for Advanced
Neuroimaging, Charité - Universitätsmedizin Berlin, 10117 Berlin, Germany
Jacek P. Dmochowski  Department of Biomedical Engineering, City University of
New York, New York, NY 10031, USA

###### Abstract

How does one find dimensions in multivariate data that are reliably expressed
across repetitions? For example, in a brain imaging study one may want to
identify combinations of neural signals that are reliably expressed across
multiple trials or subjects. For a behavioral assessment with multiple
ratings, one may want to identify an aggregate score that is reliably
reproduced across raters. Correlated Components Analysis (CorrCA) addresses
this problem by identifying components that are maximally correlated between
repetitions (e.g. trials, subjects, raters). Here we formalize this as the
maximization of the ratio of between-repetition to within-repetition
covariance. We show that this criterion maximizes repeat-reliability, defined
as mean over variance across repeats, and that it leads to CorrCA or to multi-
set Canonical Correlation Analysis, depending on the constraints.
Surprisingly, we also find that CorrCA is equivalent to Linear Discriminant
Analysis for zero-mean signals, which provides an unexpected link between
classic concepts of multivariate analysis. We present an exact parametric test
of statistical significance based on the F-statistic for normally distributed
independent samples, and present and validate shuffle statistics for the case
of dependent samples. Regularization and extension to non-linear mappings
using kernels are also presented. The algorithms are demonstrated on a series
of data analysis applications, and we provide all code and data required to
reproduce the results.

###### keywords:

Inter-subject correlation, electro-encephalography, behavioral ratings, inter-
rater correlation, canonical correlation analysis

##  1 Introduction

Consider the following scenario: A group of people are shown a movie while
their neural activity is recorded. Guided by the assumption that the movie
evokes similar brain responses in different subjects, the goal is to identify
movie-related brain activity that is common across individuals [1]. Because
activity is distributed across multiple sensors, we would like to find linear
combinations of sensors that have a similar time course in all subjects. In
other words, we would like to identify underlying _components_ or factors that
have high _inter-subject correlation_.

In a different setting, the task may be to assess motor skills. A clinician
observes a group of individuals performing various tasks and provides a rating
for each person and task [2]. Typically, ratings are combined across tasks so
that each person obtains a single aggregate performance score. But how should
individual task ratings be combined? Traditionally, this is done by simple
averaging of related scores. We propose instead a data-driven approach that
aims to make the aggregate scores as consistent as possible across different
raters. This means that we want to combine scores so as to maximize _inter-
rater reliability_.

What is common in these two scenarios is the goal of identifying directions in
multivariate data sets that are reliably reproduced across repetitions.
Correlated Components Analysis (CorrCA) accomplishes this goal. CorrCA was
originally developed in the context of neuroimaging studies to extract similar
activity in multiple subjects [1], and was re-developed independently to
extract reliable brain responses across repeated trials in a single subject
[3]. Here we will show that the technique can also be used to identify
aggregate ratings with high inter-rater reliability.

More generally, CorrCA is applicable whenever a data volume is available with
dimensions $T\times D\times N$, where $N$ are repeated measures with the same
sensors. The method identifies directions in the $D$-dimensional space along
which the data maximally correlate between $N$ repeated measurements, with
correlation measured across $T$ samples. This is similar to Principal
Components Analysis [PCA, 4] in that CorrCA finds a set of $D$-dimensional
projection vectors that linearly decompose a data set (of size $T\times D$).
Instead of capturing dimensions with maximum variance within a single data set
as in PCA, CorrCA captures the dimensions with maximum correlation between the
$N$ data sets. As we will show here, the optimality criterion of CorrCA can
also be used to derive Canonical Correlation Analysis [CCA, 5], including
multi-set CCA [MCCA, 6, 7]. The methods differs in that CorrCA uses a single
set of linear projections for all data sets assuming they repeated measures
are all taken in the same space, whereas CCA and MCCA yield a different set of
projection vectors for each data-set. MCCA is well established for the purpose
of identifying similar brain activity across subjects [8, 9, 10, 11, 12, 13].
Here we will show for brain signals and rating data that CorrCA can achieve
better performance than MCCA by virtue of the smaller number of free
parameters.

The main contribution of this paper is to provide a didactic yet formal
characterization of CorrCA and with this, establish a few novel theoretical
results. In particular, we show that maximizing the correlation between
repeated measurements is indeed equivalent to maximizing the reliability of
the projected data. By “reliability” we mean the average over repetitions,
divided by the variance over repetitions, which is also a sensible definition
for signal-to-noise ratio. We establish a direct link to MCCA, as well as
Linear Discriminant Analysis (LDA). CorrCA and MCCA both maximize inter-
subject correlation, defined here the ratio of between- to within-subjects
covariance, while LDA maximizes the between-class over the within-class
variance. Surprisingly, LDA and CorrCA yield the same result for zero-mean
signals. This direct equivalence allows us to formulate a strict statistical
test for significance based on Fisher’s variance ratio (F-statistic). The
equivalence to LDA suggests a straightforward extension of CorrCA to non-
linear mappings using the well-known “kernel” approach. We also provide a
validation of CorrCA for identifying shared dimensions using simulated data
and validate the proposed tests of statistical significance. Another important
contribution of this paper is to formalize inter-subject correlation (ISC) as
a correlation metric which applies to more than two signals. ISC differs from
conventional definitions of correlation, yet seamlessly integrates into the
frameworks of Linear Discriminant Analysis, Canonical Correlation Analysis,
and the F-statistic.

Before we present the mathematical treatment, we begin with a simple two-
dimensional example to illustrate the basic concept of identifying reliable
signal components that are contained in a shared linear dimension.

Figure 1: Extraction of a correlated component from 2D time courses. A common
electrical current source, $s(t)$, is captured by voltage sensors $x_{1l}(t)$
and $x_{2l}(t)$ for subject $l$. This common source contributes to both
subjects with identical scaling $a$ and $b$: $x_{1l}(t)=as(t)+c\xi_{l}(t)$,
$x_{2l}(t)=bs(t)+d\xi_{l}(t)$, where $\xi_{l}(t)$ is a source that contributes
a different value to each subject $l$ but with the same mixing coefficients
$c,d$. Note that voltages from multiple current sources combine additively in
a resistive medium. (A) Time courses $x_{1l}(t)$ for both subjects. (B) Time
courses $x_{2l}(t)$ (C) Same data, but now individual time points are
connected (by gray line) between the subjects. The source directions $[a,b]$
and $[c,d]$ are indicated with a red and blue arrow respectively and are the
same for both subjects, i.e. they are “shared” dimensions. (D) Projection
$y_{1l}(t)$ with maximal correlation. The time courses of the two subjects are
perfectly correlated with inter-subject correlation ISC=1. (E) Projection
$y_{2l}(t)$ that is not common between subjects. (F) Along dimension $y_{1}$,
both subjects’ data are expressed identically, i.e. this is the (“reliable”)
source dimension that is shared between the two subjects. Values along
dimension $y_{2}$ differ randomly between subjects. This dimension is not
“reliable” across subjects.

### Extracting shared dimensions — an illustrative example

For this illustrative example we consider a simulation of electrical brain
signals recorded from multiple subjects. Assume that we have two sensors,
$x_{1}(t)$ and $x_{2}(t)$, which record temporal signals from each of two
subjects, as shown in Figure 1 A & B. The figure caption details the
generation of the data; briefly, there are two “sources” of signal: one has a
common waveform for both subjects, while a second source has a subject-
dependent waveform. The common waveform indroduces a correlation between
subjects at both sensors, although the exact shared waveform is not
immediately apparent in the individual time courses shown in panels A & B).
The source directions become clear from a scatter plot of the same data (panel
C, blue and red arrows). Here, for each time point $t$ we have connected the
samples of both subjects. Evidently, at any point in time, the signals of the
two subjects lie along a line (in the direction of the source not shared by
the subjects, blue arrow), but the lines are shifted in parallel (in direction
of the common source, red arrow). The locations along the gray lines are
determined by the subject-dependent waveform, while the shift of the parallel
line from one time instance to the next is determined by the common waveform.
It is important to note that while both directions are shared by the two
subjects, only one of these dimensions expresses a common waveform. Now we
select a linear transformation, ${\bf V}$, to map this data onto a component
space, ${\bf y}(t)={\bf V}^{\top}{\bf x}(t)$, such that the component signals
are maximally correlated between subjects. With this choice, the time course
for both subjects in the first component dimension, $y_{1}(t)$, are now
perfectly correlated (panel D). Thus, we have extracted the common source of
variation in these two subjects, whereas the second component dimension,
$y_{2}(t)$, now captures a waveform that is not necessarily shared between
subjects (panel E). As a result, the gray lines in the scatter plot are now
vertical (panel F), separating the waveforms into reproducible signal
($y_{1}$) and non-reproducible noise dimensions ($y_{2}$). Evidently, ${\bf
V}$ captured the source directions in these data (blue and red arrows). The
objective of Correlated Components Analysis is to identify the unknown
projection matrix ${\bf V}$ by analyzing time courses ${\bf x}(t)$. The
fundamental assumption of CorrCA is that there is a reproducible signal and a
non-reproducible noise, and that the directions of the reproducible signal are
shared between subjects.

Readers wishing to omit the theoretical treatment below may proceed directly
to Section 3, which presents several applications of CorrCA to real-world data
sets.

##  2 Theory

###  2.1 Problem statement: Maximal inter-subject correlation

The objective of CorrCA is to find dimensions in multivariate data that
maximize correlation between repetitions. Denote the $D$-dimensional data by
${\bf x}_{i}^{l}\in\mathbb{R}^{D}$, where $i=1,\dots,T$ enumerates exemplars,
and $l=1,\dots,N$ represents repeated renditions of these data. In the case of
brain signals, the exemplars represent time points, so that $i$ indexes a
sequence in time of length $T$, and $l$ indexes $N$ repeated recordings
(subjects or trials). In the case of a behavioral assessment instrument, $i$
may index individuals receiving a rating, while $l$ indexes multiple raters,
or repeated ratings given to the same individuals by a single rater. In both
instances, the observations are $D$-dimensional ($D$ sensors that record brain
signals, or $D$ items rated in the assessment instrument). For simplicity, in
the following we will employ the terms “time” and “subjects” in reference to
$i$, such that correlation is measured _between subjects_ and _across time_.
When appropriate, we will rephrase the resulting expressions for the case of
assessing inter-rater correlation or test-retest correlation.

The goal is now to identify a linear combination of these $D$ measures,
defined by a projection vector ${\bf v}\in{\mathbb{R}}^{D}$:

$y^{l}_{i}={\bf v}^{\top}{\bf x}^{l}_{i},$ (1)

such that the correlation between $N$ subjects (repetitions) is maximized. We
define inter-subject correlation (ISC) as the ratio of the between-subject
covariance, $r_{B}$, and the within-subject covariance, $r_{W}$:

$\rho=\frac{1}{N-1}\frac{r_{B}}{r_{W}}.$ (2)

The between-subject covariance is a sum over all pairs of subjects, and
within-subject covariance is a sum over all subjects:

$\displaystyle r_{B}$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\sum_{k=1,k\neq
l}^{N}\sum_{i=1}^{T}(y_{i}^{l}-{\bar{y}}_{*}^{l})(y_{i}^{k}-{\bar{y}}_{*}^{k}),$
(3) $\displaystyle r_{W}$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\sum_{i=1}^{T}(y_{i}^{l}-{\bar{y}}_{*}^{l})(y_{i}^{l}-{\bar{y}}_{*}^{l}),$
(4)

where, $\bar{y}^{l}_{*}=\frac{1}{T}\sum_{i=1}^{T}y_{i}^{l}$, is the sample-
mean of the projected data for subject $l$. The factor $(N-1)^{-1}$ in
definition (2) is required so that correlation is normalized to $\rho\leq 1$
(see Appendix G.2). To simplify the notation, in our “variance” and
“covariance” definitions we generally omit common scaling factors, in this
case we omitted $(T-1)^{-1}N^{-1}$.

While we refer to $\rho$ as the inter-subject correlation (ISC), it could also
be termed the inter-rater correlation or inter-repeat correlation (IRC)
depending on the context. In Appendix G.1 we discuss this definition of
correlation and its relation to other measures of correlation, such as intra-
class correlation or Pearson’s correlation coefficient. Note that $r_{W}$ is
also the variance around the subject-mean.

By inserting Eq. (1) into definition (2) it follows readily that

$\rho=\frac{1}{N-1}\frac{{\bf v}^{\top}{\bf R}_{B}{\bf v}}{{\bf v}^{\top}{\bf
R}_{W}{\bf v}},$ (5)

where ${\bf R}_{B}$ and ${\bf R}_{W}$ are the between-subject and within-
subject covariance matrices of ${\bf x}_{i}^{k}$ defined analogously to
(3)-(4):

$\displaystyle{\bf R}_{B}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{T}\sum_{l=1}^{N}\sum_{k=1,k\neq l}^{N}({\bf
x}_{i}^{l}-{\bar{\bf x}}_{*}^{l})({\bf x}_{i}^{k}-{\bar{\bf
x}}_{*}^{k})^{\top},$ (6) $\displaystyle{\bf R}_{W}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{T}\sum_{l=1}^{N}({\bf x}_{i}^{l}-{\bar{\bf
x}}_{*}^{l})({\bf x}_{i}^{l}-{\bar{\bf x}}_{*}^{l})^{\top}.$ (7)

Here ${\bar{\bf x}}_{*}^{l}=\frac{1}{T}\sum_{i=1}^{T}{\bf x}_{i}^{l}$ is the
sample-mean vector for subject $l$. The projection vector ${\bf v}$ that
maximizes $\rho$ can be found as the solution of $\partial\rho/\partial{\bf
v}^{\top}=0$, which yields:

${\bf R}_{B}{\bf v}=\rho\,{\bf R}_{W}{\bf v}.$ (8)

Assuming that ${\bf R}_{W}$ is invertible, ${\bf v}$ is an eigenvector of
${\bf R}_{W}^{-1}{\bf R}_{B}$ with eigenvalue $\rho$. Maximal ISC is achieved
by projecting the data onto the eigenvector of ${\bf R}_{W}^{-1}{\bf R}_{B}$
with largest eigenvalue. We refer to this principal eigenvector as ${\bf
v}_{1}$. In addition to this projection with maximal ISC, there may be
additional dimensions that capture significant ISC in the data. Our goal is to
find all such projection vectors, ${\bf v}_{d},d=1\ldots D$. Thus, a natural
objective is to maximize ISC summed over all components. However, we want each
component to capture a different, independent aspect of the data. We therefore
require components to be uncorrelated, similar to the constraints enforced by
Independent Components Analysis or Canonical Correlation Analysis. (The
alternative of requiring orthogonal projection vectors, as in PCA is discussed
in Appendix A.5.) Thus, in total the goal is to maximize the following cost
function with respect to ${\bf V}=[{\bf v}_{1},\ldots{\bf v}_{K}]$:

$J({\bf V})=\sum_{d=1}^{D}\rho_{d}=\sum_{d=1}^{D}\frac{{\bf v}_{d}^{\top}{\bf
R}_{B}{\bf v}_{d}}{{\bf v}_{d}^{\top}{\bf R}_{W}{\bf v}_{d}}$ (9)

subject to the constraints that components are mutually uncorrelated, i.e.,
${\bf v}_{d}^{\top}{\bf R}_{W}{\bf v}_{c}=0$, for all $c\neq d$. The solution
to this problem is given by the eigenvector matrix ${\bf V}$ that satisfies
the eigenvalue equation (see Appendix A.4):

${\bf R}_{B}{\bf V}={\bf R}_{W}{\bf V}{\bm{\Lambda}},$ (10)

where ${\bm{\Lambda}}$ is a diagonal matrix with diagonal terms measuring the
ISC, $\rho_{d}$, of the corresponding projections of the data, $y_{d}={\bf
v}_{d}^{\top}{\bf x}$. Note that these additional projections of the data are
uncorrelated, as desired, because their correlation matrix, ${\bf R}^{y}_{W}$
— defined as in Eq. (7) but for ${\bf y}$ — is diagonal:

${\bf R}^{y}_{W}={\bf V}^{\top}{\bf R}_{W}{\bf V}={\bm{\Lambda}}_{W}.$ (11)

The fact that ${\bm{\Lambda}}_{W}$ is a diagonal matrix is a well-know
property of the generalized eigenvalue equation (10) (see Appendix A.1).

To summarize, the principal eigenvector of the between-subject over the
within-subject covariance matrices provides the projection of the data that
maximizes ISC. Subsequent smaller eigenvectors provide additional projections
of the data that are uncorrelated with one another. The corresponding
eigenvalues are equal to the ISC that each projection achieves. It follows
that the last projection captures the direction with the smallest ISC. Note
that the corresponding eigenvectors are uniquely defined, provided that the
ISC are distinct (see Appendix A.3).

In the illustrative example of Figure 1 with $D=2$, the eigenvalues were
$\rho_{1}=1$ and $\rho_{2}=0.04$, for the red and blue directions respectively
(see Figure 1 D and 1 F). Evidently the shared dimension where the two
subjects vary in unison (red) has been identified. A concrete application of
CorrCA to the problem of identifying brain activity that is shared between
subjects is described in section 3.1. A thorough evaluation of the conditions
under which the algorithm identifies the correct projection vectors ${\bf
v}_{d}$ will be presented in section 4.1 (Fig. 7). The main finding is that
the shared waveforms (signal) need to project in the same direction across
subjects. Performance degrades if there is additive noise, or equivalently, if
the directions of the non-reproducible waveforms (noise) differ between
subjects.

Given a set of extracted correlated components, the question arises as to how
many of these components exhibit statistically significant correlations. For
zero mean, independent test data we can use the F-statistic to determine
statistical significance (see section 2.4, Eq. 26). Statistical testing is
described in Appendix D, where we also treat the case of non-independent data
using shuffle statistics. The methods for statistical testing will be
validated with numerical simulations in section 4.

###  2.2 Fast computation of inter-subject correlation

In Eq. (6), the term $k=l$ is excluded from the sum over $k$. Note that the
excluded term is precisely ${\bf R}_{W}$, so that adding the two covariances
together completes a sum with all pairs, including $k=l$. We will refer to
this, therefore, as the total covariance:

$\displaystyle{\bf R}_{T}$ $\displaystyle=$ $\displaystyle{\bf R}_{B}+{\bf
R}_{W}$ (12) $\displaystyle=$
$\displaystyle\sum_{i=1}^{T}\sum_{l=1}^{N}\sum_{k=1}^{N}({\bf
x}_{i}^{l}-{\bar{\bf x}}_{*}^{l})({\bf x}_{i}^{k}-{\bar{\bf
x}}_{*}^{k})^{\top}.$ (13)

This relationship is useful because ${\bf R}_{T}$ can be simplified to:

${\bf R}_{T}=N^{2}\sum_{i=1}^{T}({\bar{\bf x}}_{i}^{*}-{\bar{\bf
x}}_{*}^{*})({\bar{\bf x}}_{i}^{*}-{\bar{\bf x}}_{*}^{*})^{\top},$ (14)

where ${\bar{\bf x}}_{i}^{*}=\frac{1}{N}\sum_{l=1}^{N}{\bf x}_{i}^{l}$ and
${\bar{\bf x}}_{*}^{*}=\frac{1}{T}\sum_{i=1}^{T}{\bar{\bf x}}_{i}^{*}$ are the
mean across subjects and the grand mean, respectively. Thus, surprisingly, to
compute inter-subject correlation, one never has to explicitly correlate pairs
of subjects, because ${\bf R}_{B}$ can be computed from ${\bf R}_{T}$ and
${\bf R}_{W}$, and neither one requires a sum over pairs of subjects. With an
order of $N$ operations, this makes the computation more efficient than the
direct implementation of Eqs. (6)-(7), which requires an order of $N^{2}$
operations. The direct implementation was used in all previous work [e.g., 1,
3, 14], but becomes problematic when one uses populations sizes reaching
hundreds of subjects or more [e.g. 15, 16].

###  2.3 Relationship to Linear Discriminant Analysis

Perhaps the naming and the additivity of the covariance matrices introduced
here — Eqs. (6), (7), (12) — reminds some readers of the between-class and
within-class scatter matrices used in Linear Discriminant Analysis (LDA). In
the current notation these are:

$\displaystyle{\bf S}_{T}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{T}\sum_{l=1}^{N}({\bf x}_{i}^{l}-{\bar{\bf
x}}_{*}^{*})({\bf x}_{i}^{l}-{\bar{\bf x}}_{*}^{*})^{\top}$ (15)
$\displaystyle{\bf S}_{B}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{T}N({\bar{\bf x}}_{i}^{*}-{\bar{\bf
x}}_{*}^{*})({\bar{\bf x}}_{i}^{*}-{\bar{\bf x}}_{*}^{*})^{\top}$ (16)
$\displaystyle{\bf S}_{W}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{T}\sum_{l=1}^{N}({\bf x}_{i}^{l}-{\bar{\bf
x}}_{i}^{*})({\bf x}_{i}^{l}-{\bar{\bf x}}_{i}^{*})^{\top}.$ (17)

Here, index $i$ enumerates the different classes and index $l$ enumerates the
exemplars in each class (see Table 2 in Appendix G.1 for an overview on the
use of indices in different algorithms). Therefore $\bar{{\bf
x}}_{i}^{*}=\frac{1}{N}\sum_{l=1}^{N}\bar{{\bf x}}_{i}^{l}$ and $\bar{{\bf
x}}_{*}^{*}=\frac{1}{T}\sum_{i=1}^{T}\bar{\bf{x}}_{i}^{*}$ are the class mean
and grand mean vectors, respectively. Scatter matrices also satisfy an
additivity rule: ${\bf S}_{T}={\bf S}_{B}+{\bf S}_{W}$ [17]. The goal of LDA
is to maximize the separation between classes, which is defined as the ratio
of between- to within-class variance:

$S=\frac{{\bf v}^{\top}{\bf S}_{B}{\bf v}}{{\bf v}^{\top}{\bf S}_{W}{\bf v}}.$
(18)

As before, the maximal separation $S$ is obtained with the eigenvector of
${\bf S}_{B}^{-1}{\bf S}_{W}$ corresponding to the largest eigenvalue $S$. As
before, we can extract additional projections that achieve separation of
classes. If we require that these projections are uncorrelated, then they are
given by the additional eigenvectors of ${\bf S}_{B}^{-1}{\bf S}_{W}$. This
concept was first proposed by [18] and generalizes the case of two classes
introduced by [19]. We refer to this as classical LDA, and note that it
maximizes the sum of variance ratios $S$ computed separately for each
projection (see A.4). An alternative approach is needed if instead the goal is
to maximize the ratio of summed variances [20], or if one prefers orthogonal
projection [21].

This optimality criterion and resulting eigenvalue problem look strikingly
similar to that of CorrCA in (5). Despite the similar naming, however, the
scatter matrices differ from the between-subject and within-subject covariance
matrices, and it is not immediately obvious how the two optimality criteria
relate to one another. There is, however, a close relationship between the two
(see Appendix B):

$\displaystyle N{\bf S}_{B}$ $\displaystyle=$ $\displaystyle{\bf R}_{W}+{\bf
R}_{B},$ (19) $\displaystyle N{\bf S}_{W}$ $\displaystyle=$
$\displaystyle(N-1){\bf R}_{W}-{\bf R}_{B}+NT{\bf S}_{M},$ (20)

where ${\bf S}_{M}$ is the covariance of $\bar{\bf x}^{l}_{*}$ across
subjects:

${\bf S}_{M}=\sum_{l=1}^{N}({\bar{\bf x}}_{*}^{l}-{\bar{\bf
x}}_{*}^{*})({\bar{\bf x}}_{*}^{l}-{\bar{\bf x}}_{*}^{*})^{\top}.$ (21)

When these sample-mean vectors are equal across subjects then the scatter
matrices are a linear combination of the between- and within-subject
covariance matrices. Note that equal sample-means, ${\bar{\bf
x}}_{*}^{l}={\bar{\bf x}}_{*}^{k}$, does not imply equal class-means,
${\bar{\bf x}}_{i}^{*}={\bar{\bf x}}_{j}^{*}$. Zero-mean signals will
naturally satisfy this condition and will typically still have different class
means. Under this assumption, we can also write a simple functional
relationship between class separation $S$ and inter-subject correlation
$\rho$:

$S=\frac{r_{W}+r_{B}}{(N-1)r_{W}-r_{B}}=\frac{\rho+(N-1)^{-1}}{1-\rho}.$ (22)

This relationship will be useful in section 2.4 where we establish a link
between ISC and the classic F-statistic (Eq. 26), which yields a parametric
test for statistical significance of ISC. Note that $S$ is monotonically
increasing with $\rho$ (for $\rho\leq 1$), because the slope of this
relationship is strictly positive:

$\frac{\partial S}{\partial\rho}=\frac{N}{N-1}\frac{1}{(1-\rho)^{2}}>0.$ (23)

This means that maximizing class separation $S$ in (18) also results in
maximal ISC $\rho$, provided that we equate the means. Thus, finding vectors
${\bf v}$ that maximize $S$ of Eq. (18) gives the same set of solutions as
maximizing $\rho$ in Eq. (5).

To see this, note that the goal of maximizing a ratio of two matrices, such as
in Eq. (5), can also be achieved by joint-diagonalization of the two matrices
(see Appendix A.1):

$\begin{split}{\bf V}^{\top}{\bf R}_{B}{\bf V}&={\bm{\Lambda}}_{B},\\\ {\bf
V}^{\top}{\bf R}_{W}{\bf V}&={\bm{\Lambda}}_{W},\end{split}$ (24)

where ${\bm{\Lambda}}_{B}$ and ${\bm{\Lambda}}_{W}$ are diagonal matrices.
Note that if ${\bf V}$ jointly diagonalizes two matrices, it also diagonalizes
any linear combinations of the two [22, briefly recapitulated in Appendix
A.2]. Because ${\bf S}_{B}$ and ${\bf S}_{W}$ are linear combination of ${\bf
R}_{B}$ and ${\bf R}_{W}$, then the ${\bf V}$ that satisfies (24) also
diagonalizes ${\bf S}_{B}$ and ${\bf S}_{W}$. Thus, this ${\bf V}$ captures
the eigenvectors to ${\bf R}_{W}^{-1}{\bf R}_{B}$ as well as ${\bf
S}_{W}^{-1}{\bf S}_{B}$, and satisfies in both cases the constraint of
mutually uncorrelated components. Because the eigenvalues of the two problems
are monotonically related, the eigenvectors will be sorted in the same order
for both problems. Thus, the solutions to both optimization problems are the
same, including the constraint of uncorrelated components. To summarize,
CorrCA and classical LDA yield the same result, provided that the sample mean
vectors are equal across subjects.

The equivalence between CorrCA and LDA is perhaps surprising. Note that LDA
attempts to separate classes so that exemplars of one class do not overlap
with the exemplars of other classes. In the illustrative example of Figure 1,
time samples $i$ take on the role of classes. There are 20 classes, each with
two exemplars (connected with gray lines). Classes are well separated in the
direction of component $y_{1}$ (red arrow), whereas they are overlapping in
the original dimensions $x_{1}$ and $x_{2}$. Note that the signals are zero-
mean across samples but have differing class-means. In this toy example, there
is no variation between the two exemplars in each class (vertical lines in
panel F indicate that the $y_{1}$ values are identical within a class). Thus,
CorrCA has found a dimension of the data ($y_{1}$) in which the classes are
perfectly separated. Zero within-class variance leads to infinite separation
(Eq. 18) and this in turn to unit correlation (Eq. 22). More generally, we
have demonstrated with Eq. (22) that increasing separation between samples is
equivalent to increasing correlation across samples. To understand this, note
that separation is high when multiple repeats have similar values, yet these
values have to be different for different samples. This means that multiple
repeats have to vary together across samples. In other words, the variations
across samples have to be highly correlated across repeats.

Notice that the equivalence between CorrCA and LDA requires the same number of
exemplars in each class – something that is not usually the case in
conventional LDA. In fact, there is a one-to-one link between exemplars across
classes (gray lines in panels C and F of Figure 1).

Table 1: Relationship between the optimality criteria of different algorithms.

CorrCA |  LDA |  JD
---|---|---
$\rho=\frac{1}{N-1}\frac{{\bf v}^{\top}{\bf R}_{B}{\bf v}}{{\bf v}^{\top}{\bf R}_{W}{\bf v}}$ | $S=\frac{{\bf v}^{\top}{\bf S}_{B}{\bf v}}{{\bf v}^{\top}{\bf S}_{W}{\bf v}}\stackrel{{\scriptstyle*}}{{=}}\frac{\rho+(N-1)^{-1}}{1-\rho}$ | $\frac{{\bf v}^{\top}{\bf S}_{B}{\bf v}}{{\bf v}^{\top}{\bf S}_{T}{\bf v}}\stackrel{{\scriptstyle*}}{{=}}(N-1)\rho+1$

  * •

$\rho$: inter-subject correlation, $S$: class-separation. * assuming zero-mean
signals, or unbiased raters (${\bf S}_{M}=0$).

###  2.4 Maximum reliability and the F-statistic

In the present context, the scatter matrices ${\bf S}_{B}$ and ${\bf S}_{W}$
are interesting for another, perhaps more important reason. Consider the case
where $l$ represents repeated measurements of raters, or repeated measures on
the same subjects, so that $\rho$ now measures inter-rater correlation, or
inter-repeat correlation (IRC). According to Eq. (16), ${\bf S}_{B}$ measures
(except for a scaling factor) the sample-covariance of ${\bar{\bf
x}}_{i}^{*}$, which is the mean across repeats. On the other hand, according
to Eq. (17), ${\bf S}_{W}$ measures the sample-covariance averaged over
repeats. When projected onto $y$ with ${\bf v}$ they define the variance of
the mean across repeats, $\sigma_{\bar{y}}^{2}$, and average variance around
these means, ${\bar{\sigma}_{y}}^{2}$, respectively. The ratio of the two
variances is a sensible definition for signal-to-noise ratio (SNR):

$S=\frac{\sigma_{\bar{y}}^{2}}{\bar{\sigma}_{y}^{2}}.$ (25)

We take this SNR also as a metric of repeat-reliability. In this view, the
results of the previous section show that maximizing correlation between
repeats is equivalent to maximizing repeat-reliability. In particular, Eq.
(22) provides the relationship between SNR and ISC. An application of CorrCA
to the problem of identifying ratings that are reliably reproduced between
different raters is described in Sections 3.2 and 3.3. In Section 3.4 we
describe results on the problem of identifying components of brain activity
that is reliably reproduced across repeated trials. In all instances, the goal
is to maximize $\rho$, which is equivalent to maximizing $S$.

A component decomposition method that specifically maximizes SNR was
previously described in de Cheveigné and Parra [23]. This method can be
formulated as a joint diagonalization (JD) problem similar to CorrCA or LDA.
When SNR is defined as in Eq. (25), then the JD approach diagonalizes matrices
${\bf S}_{T}$ and ${\bf S}_{B}$ in the present notation. By the same argument
as in section 2.3, it is clear that this JD approach provides the same
solution as LDA and CorrCA. The optimization criteria for the three techniques
are summarized in Table 1.

We note that $S$ is a ratio of variances, which permits a direct link to the
$F$-statistic, as used in the analysis of variance:

$F=\frac{T(N-1)}{T-1}S=\frac{T(N-1)\rho+T}{(T-1)(1-\rho)}.$ (26)

For normally and independently distributed samples this quantity follows the
$F$-distribution with degrees of freedom $d_{1}=T(N-1)$ and $d_{2}=T-1$ [24,
Chapter 6]. We can therefore use the $F$-distribution to assess statistical
significance for $\rho$. This gives us a strict parametric test of
significance for each CorrCA component (on unseen data and provided the means
are equalized, i.e. $\bar{\bf x}^{l}_{*}=\bar{\bf x}^{*}_{*}$). We will
leverage this observation when we perform tests for statistical significance
of $\rho$ evaluated on independently and identically distributed test data.
For non-independent samples, which are more typical of the temporal signals of
interest, we will rely instead on shuffle-statistics, which can also be
applied directly to training data (see Appendix D). Both methods will be
evaluated for accuracy in section 4.

###  2.5 Forward model

An import aspect of multivariate models is parameter interpretation. CorrCA is
an example of a “backward model” in the sense that the observed data is
transformed into components that capture the source of covariation. It shares
this property with a host of other methods such as PCA and ICA. An alternative
strategy is “forward modeling” such as Factor Analysis or General Linear
Models, which transforms known variables to explain the observed data [25].
Such models often capture the physical processes underlying data generation.
For instance, in the case of electromagnetic brain signals, a forward model
refers to the “image” that a current source in the brain generates at the
scalp sensors.

For CorrCA, the backward model is given by the matrix ${\bf V}$ which maps
observations ${\bf x}$ to components ${\bf y}={\bf V}^{\top}{\bf x}$. A
corresponding forward model can be defined as the projection ${\bf A}$ that
best recovers measurements ${\bf x}$ from the components ${\bf y}$ [26, 27]:

$\hat{\bf x}={\bf A}{\bf y}.$ (27)

The least-squares estimate of this forward model projection is

$\hat{\bf A}={\bf R}_{W}{\bf V}({\bf V}^{\top}{\bf R}_{W}{\bf V})^{-1}.$ (28)

Note that the columns ${\bf a}_{k}$ of this matrix capture how correlated a
putative source ${\bf y}_{k}$ is with the observed sensor recordings ${\bf
x}$. This approach of recovering the forward model from the backward model is
explained in more detail in [25], along with a discussion on how forward and
backward model parameters need to be interpreted. In the illustrative example
of Figure 1, the forward model is shown in panel C as red and blue arrows for
first and second component respectively. Evidently, the directions recovered
by CorrCA match the underlying source orientations described in that example.

The utility of the forward model will be demonstrated for the case of brain
activity, where we are interested in how the activity of a recovered neural
source manifests at the sensor level (Figure 2 and Figure 5). These forward
models can be used directly to perform source localization to identify the
spatial origin of the corresponding current sources [25]. In the case of
questionnaire ratings (section 3.3) we will interpret the meaning of a
component by inspecting the forward model, as it measures the correlation of a
component with answers of raters on different topics. For the case of ratings,
we may also be interested in which items lead to a reliable aggregate score
(sections 3.2). In these cases, an analysis of backward model weights may
provide valuable insight on what data features are most important to construct
a reliable aggregate measure.

Note that if the original data is uncorrelated, then ${\bf R}_{W}$ is diagonal
and the projections ${\bf V}$ are orthogonal. In that case the forward and
backward models are identical, except of an overall scale for each component.
In that case, we can directly inspect the backward model, as is customarily
done in PCA, which has orthogonal projections and one inspects the component
“loadings”. In the example of Figure 4D the covariance ${\bf R}_{W}$ is
approximately diagonal and one can inspect the results with either the forward
of backward model.

###  2.6 Relationship to multi-set Canonical Correlation Analysis

Thus far we have assumed that all repetitions should receive a common
projection vector. This is sensible in the case of multiple trials if we
expect that the sources of interest behave similarly across repetitions.
However, in the case of diverse subjects or raters, it may be more appropriate
to assume that each repetition should be uniquely weighted. For instance,
brain activity recorded at a given location may differ between subjects due to
anatomical or functional differences. This brings us to the case treated by
Canonical Correlation Analysis (CCA) for $M=2$, or more generally, multi-set
CCA (MCCA) for $M>2$. There are multiple implementations of MCCA [7, 28]. We
will derive one of these instantiations as a maximization of ISC as defined in
Eq. (2). Assign to every subject $l$ a unique projection vector ${\bf v}^{l}$:

$y^{l}_{i}={\bf v}^{l\top}{\bf x}^{l}_{i}.$ (29)

By inserting this into definition (2), it follows readily that

$\rho=\frac{\sum_{l=1}^{N}\sum_{k=1,k\neq l}^{N}{\bf v}^{l\top}{\bf
R}^{lk}{\bf v}^{k}}{(N-1)\sum_{l=1}^{N}{\bf v}^{l\top}{\bf R}^{ll}{\bf
v}^{l}},$ (30)

where ${\bf R}^{lk}$ is the cross-covariance matrix between ${\bf x}_{i}^{l}$
and ${\bf x}_{i}^{k}$:

${\bf R}^{lk}=\sum_{i=1}^{T}({\bf x}_{i}^{l}-{\bar{\bf x}}_{*}^{l})({\bf
x}_{i}^{k}-{\bar{\bf x}}_{*}^{k})^{\top}.$ (31)

The projection vectors ${\bf v}^{l}$ that maximize $\rho$ can be found as the
solution of $\partial\rho/\partial{\bf v}^{l\top}=0$, which yields for each
$l$ the following equation:

$\frac{1}{N-1}\sum_{k=1,k\neq l}^{N}{\bf R}^{lk}{\bf v}^{k}=\rho\,{\bf
R}^{ll}{\bf v}^{l}.$ (32)

This set of equations for all ${\bf v}^{l}$ can be written as a single
equation via concatenation to a single $ND$-dimensional column vector ${\bf
v}=\left[{\bf v}^{1\top}\ldots{\bf v}^{N\top}\right]^{\top}$:

${\bf R}{\bf v}=\lambda\,{\bf D}{\bf v},$ (33)

where, $\lambda=(N-1)\rho+1$, ${\bf R}$ is a matrix combining all ${\bf
R}^{lk}$, and where the covariance matrices ${\bf R}^{ll}$ are arranged in a
block-diagonal matrix ${\bf D}$:

${\bf R}=\left[\begin{array}[]{cccc}{\bf R}^{11}&{\bf R}^{12}&\ldots&{\bf
R}^{1N}\\\ {\bf R}^{21}&{\bf R}^{22}&\ldots&{\bf R}^{2N}\\\
\vdots&\vdots&\ddots&\vdots\\\ {\bf R}^{N1}&{\bf R}^{N2}&\cdots&{\bf
R}^{NN}\\\ \end{array}\right],{\bf D}=\left[\begin{array}[]{cccc}{\bf
R}^{11}&0&\ldots&0\\\ 0&{\bf R}^{22}&\ldots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&{\bf R}^{NN}\\\ \end{array}\right].$
(34)

If ${\bf D}$ is invertible, the stationary points of the ISC are now the
eigenvectors of ${\bf D}^{-1}{\bf R}$, or more generally, the generalized
eigenvectors of Eq. (33). We can arrange all such eigenvectors as the columns
of a single matrix ${\bf V}=[{\bf v}_{1}\ldots{\bf v}_{n}\ldots{\bf v}_{D}]$.
The eigenvector with the largest eigenvalue $\lambda$ maximizes the ISC
because $\rho$ increases linearly with $\lambda$. The subsequent eigenvectors
maximize ISC subject to the constraint that the component signals are
uncorrelated from the previous components (see Appendix A.4). In other words,
the solution ${\bf V}$ diagonalizes ${\bf D}$:

${\bf V}^{T}{\bf D}{\bf V}={\bm{\Lambda}}.$ (35)

Formulating the MCCA problem as an eigenvalue problem can be derived in
multiple ways [29]. In Kettenring [7] it is presented as the maximization of
the eigenvalue of the correlation matrix between subjects (MAXVAR criterion),
but it can also be derived by maximizing the sum of correlations between
subjects (SUMCORR criterion) [28]. Here we have derived MCCA as a maximization
of ISC. In this sense CorrCA can be seen as a version of MCCA constrained to
have identical projection vectors for all $N$ data sets (repeats/subjects).
MCCA has $ND$ parameters per component, whereas CorrCA has only $D$. When
there is an abundance of data ($T\gg ND$) and there is reason to believe that
the $N$ datasets should be weighted individually, then MCCA may be appropriate
[e.g., 30]. However, when $N$ is large compared to the data size ($N>T/D$)
CorrCA may be preferable due to the reduced number of free parameters. This
will be demonstrated in section 3.1.

When there are only two data sets ($M=2$) then the present MCCA algorithm
reduces to conventional CCA [31, see Table 4.1]. In the illustrative example
of Figure 1, CorrCA by definition gives the same result for both subjects (red
and blue arrows), but we find that the directions derived with CCA differ
between subjects (not shown). In section 3.3, we will apply CCA to ratings
data from two distinct raters (pairs of parent and child), in which case
distinct projections may be warranted.

##  3 Applications

In the following sections we demonstrate existing and new applications of
CorrCA, while also examining issues related to statistical significance
testing and regularization of the model. We also compare CorrCA with MCCA/CCA
and kernel-CorrCA.

###  3.1 Extracting brain activity that correlates between subjects

We have previously used CorrCA to analyze brain activity recorded from
multiple subjects [1, 32, 33, 14, 34, 35, 15, 36]. Briefly, subjects in these
experiments watched movies or listened to narrated stories while their neural
responses were recorded with multiple sensors (electroencephalography (EEG) —
electrical potentials measured at multiple locations on the scalp surface).
CorrCA was then used to extract components that were most reliable across
subjects. These experiments showed that brain signals are significantly
correlated between subjects on a rapid time-scale of less than a second
(signals were high-pass filtered at 0.5 Hz). While the ISC values are small
($\rho\approx 0.05$), these values are statistically significant and very
reliable across sessions [1] even in realistic environments [34].
Interestingly, the level of correlation between subjects as measured by ISC is
indicative of how attentive subjects are towards the stimulus [33] and,
therefore, predictive of how well individuals remember the content of the
videos [14, see Appendix G.3 for a definition of ISC for individual subjects].
ISC of the EEG also predicts the retention of TV and online audiences [32,
35].

Figure 2:  Correlated components of the EEG evoked by natural video stimuli.
(A) The ISC achieved by both CorrCA (blue bars) and MCCA (green bars) on EEG
data from 18 subjects viewing a 197-second video clip (“Sunday at Rocco’s”,
2007, Storycorp). Here MCCA was regularized with TSVD ($K=12$; see Appendix
C), while CCA was not regularized. Darker shades indicate that the learned
components were found to be statistically significant using circular shuffle
statistics at $\alpha<0.05$ (see Appendix D): 3 significant components were
found with CorrCA, while 1 was identified with MCCA. Light shades indicate
non-significant ISC. (B) ISC measured on independent test data. Test-set
performance was computed by training on 4/5 of the time samples and testing on
the 1/4 of left-out samples.. Here CorrCA produces higher ISCs than MCCA
across 4 of the top 5 components (only the first is significant) (C) The
spatial distributions of the first three correlated components displayed here
as “forward models” (see section 2.5). The three components exhibit smooth
topographies with unique spatial profiles, suggesting that each component
recovers genuine and distinct neural activity. (D) The forward models of the
subject-specific components learned by MCCA. Each column displays the forward
models of the components from an individual subject (i.e., subjects 1, 5, 9,
13, and 17). The components are markedly less smooth, consistent with over-
fitting, as MCCA learns separate parameters for each of 18 subjects. For
details on this data set, please refer to Cohen and Parra [14].

Here we reanalyzed the data of [14] to compare CorrCA with MCCA (see section
2.6) on their ability to learn reliable projections of the data, while also
using the opportunity to demonstrate forward model computation (section 2.5)
and statistical testing (Appendix D). The data set consisted of the EEG from
$N=18$ subjects as their brain activity was captured with $D=64$ scalp
electrodes across $T=50344$ time samples corresponding to the 197-second video
stimulus. We first compared CorrCA with regularized MCCA in terms of the
number of reliable dimensions extracted by each, computed here using a non-
parametric test using circular shuffle statistics (see Appendix D). The ISC
values found on the entire data set are shown for the first 5 components in
Figure 2 A , where blue bars correspond to CorrCA and green to MCCA. Moreover,
darker shading indicates that the ISC of the component was found to be
statistically significant: 3 significant dimensions were found by CorrCA,
while only 1 was found by MCCA, despite the (training set) ISC values being
higher for the latter. Indeed, when evaluating ISC on unseen test data the ISC
drops substantially for MCCA as compared to the training data, and CorrCA
yielded higher ISC values than MCCA for 4 of the first 5 dimensions (Figure 2
B).

Both of these findings are indicative of over-fitting with MCCA, despite the
use of regularization. We used the Truncated Singular Value Decomposition
(TSVD, similar to what is described in Appendix C) and selected a subspace of
$K=12$ dimensions as this gave the largest test-set ISC (summed over the first
3 components). Note that MCCA fits a model for each subject. In this example
it uses 216 parameters per component ($N=18$ subjects, $K=12$ dimensions),
whereas CorrCA uses $D=64$ parameters per component on these data. While
stronger regularization of MCCA is possible, the regularization is agnostic to
the shared dimensions, which explains why ISC dropped with stronger
regularization on the test data. CorrCA can be seen as a form of
regularization that addresses over-fitting, without removing dimensions that
significantly correlate between subjects.

Both CorrCA and MCCA learned spatial filters that were applied to the multi-
electrode EEG data to maximize the ISC. The activity captured by each spatial
filter (component) is best depicted by “forward models” which are illustrated
row-wise for the first three correlated components in Figure 2 C & D. The
components extracted with CorrCA are marked by a smooth topography consistent
with activation of the cerebral cortex (Figure 2 C). The components produced
by MCCA are less smooth (Figure 2D shows subject-specific forward models for
the first three components of 5 of the 18 subjects, i.e., subjects 1, 5, 9,
13, and 17). Note that only the first component replicates, for some subjects,
the bilateral occipital activation indicative of visual cortex activity seen
in panel C (top). Other topographies are unlikely to be generated by a genuine
neural source, and instead further support the assertion that MCCA overfit the
data.

Here and in the following examples we computed model parameters ${\bf v}$ by
leaving out samples $i$, and measured test-set ISC on the left-out samples.
The alternative could have been to leave out subjects. However, the number of
entries in the between-subject correlation matrix ${\bf R}_{B}$ defined in Eq.
(6) scales quadratically with the number of subjects $N$, but linearly with
number of samples $T$. Thus, we expect the accuracy of ${\bf R}_{B}$ and the
resulting eigenvectors ${\bf v}$ in Eq. (8) to decrease faster when removing
subjects instead of samples.

###  3.2 Maximizing inter-rater reliability

The theory laid out in section 2.4 shows that one can use CorrCA to identify
components with maximal inter-rater reliability. In this following example,
clinicians assessed the motor skill of children by measuring performance on
various standard tasks (Zurich Neuromotor Assessment). The rating obtained for
each task here is the time taken by the child to perform each task. This time
is normalized relative to the values of a standard group of children of the
same age [37]. Although this is an objective measure, there is nonetheless
variability due to variable performance of the child when repeating the task,
and also subjective variability in the observer operating the stop-watch. In
this specific data set, we have ratings from $T=30$ children on $D=12$ tasks
from $N=2$ raters. The 12 measures are mostly positively correlated across the
children (Figure 3 A, see panel G for descriptive task labels). Reliability
between raters is measured here as the IRC, which differs across tasks (Figure
3 B).

Figure 3: Ratings from the Zurich Neuromotor Assessment for $T=30$ children,
$D=12$ tasks that were rated, and $N=2$ clinical raters. The ratings are
continuous numbers and normed by age. The goal of CorrCA is to identify
aggregate scores that have high correlation between raters (IRC). (A, B, C)
Original rating data. (D, E, F) Components extracted by CorrCA. One of the two
raters provided two repeated ratings for each child so that one can evaluate
correlation between the repeated tests (test-retest). (A, E) show IRC between
the two different raters. (B, F) show IRC between the two repeated ratings of
the same rater. Dark-shaded bars indicated statistically significant IRC (in
panel E using circular shuffle statistics for training data; in panels B, C, F
using F-statistics with Bonferroni correction; $\alpha<0.05$; see Appendix D)
(G) Projection matrix ${\bf V}$ for the first 5 significant components. Data
were provided by Valentin Rousson [37]. See text for details on each panel.

For many assessment instruments, individual ratings are summed up to obtain a
total score. It is also natural to compute sub-scores by combining ratings on
related items. In these sums, ratings effectively obtain equal weights, which
is a sensible but nonetheless arbitrary choice. CorrCA can provide weightings
that will give the most reliable response across raters, as reproducibility of
results is a desirable goal. In this specific example, the goal is to find the
weightings of the 12 motor tasks that provide reliable aggregate (sub-)
scores, and which capture independent aspects of these ratings (independent to
second order). The resulting weighting vectors ${\bf V}$ are shown in Figure 3
G. Each column corresponds to a different component and has been normalized to
have unit norm. Components are sorted in decreasing order of IRC (Figure 3 E)
and therefore decreasing order of inter-rater reliability (see section 2.4,
Eq. 25). Of these components, the first five reach statistical significance at
an $\alpha$-level of 0.05 (highlighted in Figure 3 E). Inspecting the
“backward model” ${\bf V}$ it is clear that first component relies heavily on
a single task (task 2, which has the highest IRC) but has small positive and
negative contributions from other tasks. Components 2 and 3 are mostly a
contrast between tasks 3 and 8, and 2 and 4, respectively. The detailed
results can be obtained by running the code accompanying this article. In
total, the method extracts different aggregate scores of high reliability,
relying in some instances on a small subset of tasks. This has the potential
of simplifying the assessment of motor performance while increasing its
reliability. An important aspect of these various components is that, by
design, they are uncorrelated from one another (uncorrelated across the 30
children; Figure 3 D), whereas the original task ratings all measure similar
aspects (they are predominantly positively correlated; Figure 3 A).

In this data set, children were also evaluated a second time by the same rater
in order to assess test-retest reliability (Figure 3 C). Comparing panels B
and C, it is evident that the test-retest reliability is somewhat higher than
the inter-rater reliability (i.e., the evaluation does indeed depend to some
degree on the subjective judgment of the rater). The components extracted with
the highest IRC also exhibit the strongest test-retest correlations (at least
for the first few components). This confirms that these aggregate measures
generalize to data that was not used in the CorrCA optimization, and that
minimizing subjective variability of the observer also minimizes test-retest
variability.

###  3.3 Identifying agreement between child and parent

We also provide an additional example of CorrCA on subjective ratings data
related to the style of parenting. In this example, the goal is to identify
questions on which parent and child agree. Specifically, we analyzed data
collected with the Alabama Parenting Questionnaire [APQ, 38] as part of the
Healthy Brain Network initiative [16]. In the APQ, parents and children
answered $D=42$ questions pertaining to their relationship by providing a
numerical value from 1 to 5 for each question. From this, the APQ averages
various questions to provide aggregate scores related to different aspects:
“involvement score”, “positive parenting score”, “poor supervision score”,
“inconsistent discipline score”, and “corporal punishment score”. At the time
of this analysis, complete data were available for $T=616$ children/parent
pairs (ages 6-17, 40% female), i.e. $N=2$ raters. The first observation from
the data is that the answers from children and parents are only weakly
correlated ($\rho=0.17\pm 0.1$, Figure 4 A). The aggregate scores fare even
worse ($\rho=0.064\pm 0.13$). A two-sample t-test suggests that these
aggregate scores are less reliable than individual questions
($t(47)=2.4,p=0.022$) including several scores with negative correlations.
Evidently, children do not much agree with their parents when it comes to
parenting issues, despite answering nearly identical questions on the topic.

We applied CorrCA to find a combination of questions on which parents and
children agree. To this end, we divided families at random into a training and
a test set of equal size ($T=308$) and computed the projections ${\bf V}$ that
maximized IRC on the training set. We then evaluated the resulting projections
in terms of the correlation between child and parent for both train and test
sets (Figure 4 B & C). By design, IRC on the training set is largest for the
first component. On the test set, the first component also has a relatively
large IRC, whereas the other components do not differ substantially in IRC
from the original set of questions (IRC hovers around 0.2 in the original
ratings, ${\bf X}$, and the aggregate component scores, ${\bf Y}$). The lower
IRC values on the test set as compared to the training set indicate the
presence of over-fitting, although the discrepancy has been reduced here using
shrinkage regularization (using $\gamma=0.5$; see Appendix C). To determine
how many significantly correlated dimensions are actually present in these
data, we used circular shuffle statistics (see Appendix D) on all the data
($T=616$) and found $K=8$ components with $p<0.05$. On the test data alone,
computation of the F-statistic revealed $K=8$ significant components (Figure 4
B & C, dark shaded).

Figure 4: Ratings from the Alabama Parenting Questionnaire for $T=616$
families and $D=42$ questions, with $N=2$ ratings from a child and one of its
parents. Here, the goal of CorrCA is to maximize IRC, i.e., to identify
agreement between child and parent. Data were collected by the Healthy Brain
Networks Initiative [16]. Data were divided evenly and at random into a
training and a test set. Second and third columns show results for training
and test sets, respectively. Red error bars indicate the standard deviation
across 100 random partitions of the data into training and test sets, and
sampling with replacement within that (bootstrap estimates). Dark-shaded bars
indicate significant IRC at $\alpha<0.05$ (F-statistic, Bonferroni corrected
with $T=308$ for panels A and C; and circular shuffle statistic with $T=616$
for panels B and E). (A) IRC between parent and child on all original
questions. Dark-shaded bars indicate significant IRC (tested on Pearson’s
correlations, $\alpha=0.05,T=616$. (B) IRC for the components extracted with
CorrCA on the training set. Dark-shaded bars indicate significant IRC
(circular shuffle statistic, $\alpha<0.05,T=616$). (C) IRC on the test set.
Only the first component achieves strong IRC. In this case, we used shrinkage
regularization ($\gamma=0.5$, see Appendix C). Dark-shaded bars indicate
significant IRC (F-statistic, Bonferroni corrected, $\alpha<0.05,T=308$) (D)
Correlation between questions across all $616$ families. (E) IRC for
components extracted with MCCA (section 2.6) and kernel-CorrCA (Appendix F).
Dark-shaded bars indicate significant IRC (circular shuffle statistic,
$\alpha<0.05,T=308$) (F) Comparison of IRC on the test set for the 4 strongest
components extracted with CorrCA, MCCA and kernel-CorrCA. (G) Component
weights for the strongest component extracted with CorrCA.

The component with the highest IRC has a combination of questions contributing
to reasonably reliable aggregate scores ($\rho=0.61$ on test set, Figure 4 G).
To interpret this component we inspect its forward model weights (see section
2.5.) The following six questions correlate most strongly to the first
component (weightings of ${\bf V}$ are given in parenthesis): Question 32:
“You are at home without an adult being with you” (-0.38); Question 11: “Your
parent helps with your homework” (+0.35); Question 19: “Your child goes out
with a set time to be home” (-0.33); Question 21: “You go out after dark
without an adult with you” (-0.33); Question 41: “You use time-out as
punishment” (+0.33). Question 6: “You fail to leave a note or let your parents
know where you are going” (-0.28). The sign of the forward weights indicate
here whether the answer to that question correlate positively or negatively
with the aggregate measure [25]. Therefore, this measure seems to capture the
independence of the child versus parental involvement. The positive
correlation of Question 41 with this independence-vs-involvement dimension is
interesting as it suggests that this form of punishment is associated with a
child’s independence. In total, while parents and children largely disagreed,
there is moderate agreement on how independent their lives are. Given this
result, we suspected that this metric increases with the child’s age. Indeed,
there is a strong correlation of this component with the child’s age across
the cohort (Pearson’s $r=0.7$ and $r=0.67$ for parent and child scores
respectively, $p<10^{-80},T=616$; with no significant difference with the
child’s sex). We did not analyze subsequent components due to the proximity of
residual eigenvalues (Figure 4, panel C) and noise level may not permit
accurate identification (see section 4.1).

Given that parents and children may genuinely differ in their judgments, it
would make sense to allow for different weighting to the questions. Therefore,
we also applied MCCA to this data as described in section 2.6. We found only
$K=4$ significant components at $p<0.05$ (as opposed to $K=8$ for CorrCA;
Figure 4E, dark shaded). ISC summed over these first four components was
numerically smaller for MCCA vs CorrCA (0.69 vs 1.10, values obtained on the
test data without regularization; Figure 4 F). Therefore, it appears that on
this data the sample size of $T=616$ was not sufficiently large to support a
doubling in the number of free parameters (from 42 for CorrCA to 84 for MCCA).

We have also tested the non-linear kernel-CorrCA method (see Appendix F) on
this ratings dataset (Figure 4 E). There was a small gain in terms of inter-
rater correlation or inter-trial correlation respectively in the test data
(Figure 4 F). However, the projections of the strongest components were not
substantially different from the linear projections, suggesting that there
were no strong non-linear relationships in the data. We find a similar result
with kernel-CorrCA for the the data presented in the next section 3.4 (results
not shown).

###  3.4 Identifying stimulus-related components in EEG with high SNR

PCA is often used to extract a low-dimensional representation of high-
dimensional data. To this end, PCA selects orthogonal projections that capture
most of the variance of the projected samples. When multiple repetitions of
the data are available, CorrCA can be used for the same purpose of reducing
dimensionality of the data, but with the objective of capturing dimensions
that are reliably reproduced across repetitions. Here we apply this to a
conventional scenario in neuroimaging: neural activity is collected over time
with multiple sensors, while the same experimental stimulus is repeated in
multiple trials. We propose to use CorrCA to identify spatio-temporal patterns
of activity that are well preserved across trials, similar to what was
proposed in [3, 23, 39, 40]. Ideally, the corresponding time-courses of
multiple components are uncorrelated in time so that they capture different
activity. CorrCA guarantees this by design as it diagonalizes ${\bf R}_{W}$.
CorrCA is therefore similar to other source separation methods such as
Independent Components Analysis [41], Denoising Source Separation [42], or
Joint-Decorrelation [23]. In fact, CorrCA can be viewed as a form of blind-
source separation, as no labels are required for the analysis.

Here we analyze EEG data that were collected during an Eriksen Flanker Task
[26]. Subjects were presented with a collection of arrows pointing left or
right. They were asked to decide as quickly as possible between two
alternative choices (the direction of the center arrow) by pushing
corresponding buttons with either their left or right hand. We analyze data
from a single participant who responded correctly in 550 trials, and made an
error in 46 trials (Figure 5). This error, which participants notice as soon
as they make their incorrect response, leads to a negative EEG potential 0-100
ms after the button press (this period is indicated between black and red
vertical lines in panel B & C for electrode #7, see electrode location in
panel H). This phenomenon is known as the “error-related negativity” [ERN,
43]. Button pushes also elicit the “readiness potential” — a gradually
increasing potential preceding the button push. Aside of the readiness
potential and error-related negativity there are a series of other potential
deflections following the button push. In this example, we have $T=176$ time
samples, $D=64$ electrodes, and $N=46$ trials. Time is aligned across trials
so that the button press is at time $t=0$ ms. The traditional approach to
analyzing such data is to average the electric potentials over trials,
yielding what are called event-related potentials (ERPs, panel B). The
reliability of these trial-averages is relatively low (panels B& D show mean
and standard-error across trials, and panels C & E show single trial data). We
selected electrodes #7 and #22 for display as these capture distinct time
courses and locations identified in prior literature. The SNR (as defined in
section 2.4) is maximal for the same cluster of frontal electrodes (#12, #17,
#44, #43) and the mean SNR across all electrodes is $S\sqrt{N}=0.72$ (Figure 5
A). We are multiplying SNR with $\sqrt{N}$ because, typically, the goal is to
determine if the mean value differs significantly from zero, and thus
standard-error of the mean (SEM) is the relevant quantity.

Figure 5: Event-related potentials. EEG recordings from a single participant
collected during a speed two-alternative forced-choice task (Eriksen Flanker
Task). (A-E) Original EEG data $x_{j}(t)$, highlighting two channels, $j=7$
(Fz electrode, SNR=0.97) where the “error-related negativity” is expected to
be most pronounced (between 0 and 100ms; black and red vertical lines), and
$j=22$ (PCz electrode, SNR=1.3) where a “readiness potential” is typically
observed. (F-L) First and second correlated components with the strongest
inter-trial correlation. (A & F) SNR computed from the inter-trial correlation
(Eq. 25) for original data and extracted components respectively. Dark shaded
bars indicate significant components (Circular shuffle statistic on the
training set; $\alpha<0.05$; see Appendix D) (B, D, G, F) Average across
$N=46$ repeated trials. Blue-shaded area indicates standard-error of the mean.
Evidently the extracted components have a smaller SEM relative to the mean
values. The SNR measures are multiplied with $\sqrt{N}$ to parallel the SEM
shown here. SNR of 2, for instance, indicates that the mean is 2 SEM away from
zero baseline. (C, I, E, F) EEG resolved over all trials to visualize the
inter-trial variability. (H & K) Forward model of the two strongest components
with distributions consistent with the ERN and readiness potential
respectively. Data from [26].

In our earlier work we specifically looked for a projection that distinguishes
the EEG evoked response between error and correct responses [26]. Here the
goal is to identify reliable error responses. When we apply CorrCA to these
data, we obtain a series of components with high inter-trial correlation
(Figure 5 F). CorrCA components were extracted using shrinkage regularization
($\gamma=0.4$, see Appendix C) on a training set of $N=45$ trials (leaving 1
trial out for testing). Performance was evaluated on the leave-one-out test
set. The first three components (highlighted in dark blue) are significant at
the $\alpha=0.05$ level according to a non-parametric test based on circularly
shifted surrogate data (see Appendix D). The first two of these components
have an SNR of $S\sqrt{N}=3.3$ and $1.6$. Evidently, component 1 is very
reliable (panel G shows that the mean is several SEM away from the baseline).
The deflections are now obvious even in single trials without the need to
average (panel I). The spatial distribution of this component activity (panel
H) is indicative of the well-known ERN. Component 2 (panel J) is indicative of
the readiness potential [44]. It rises leading up to the button push at time
$t=0$ ms, and has a spatial distribution over parietal and motor areas (panel
K).

Note that CorrCA has reduced the $D=64$ dimensional data to just a few
components that capture the reliable portion of the data (high SNR; panel F).
By compactly representing neural activity in a small number of components, it
becomes easier to probe for changes in activity due to manipulation of
experimental variables, whose effects on the original $D=64$ dimensional data
may be obscured due to the lower SNR and the multiple comparison problem.

Like most data decomposition algorithms, CorrCA does not require that the
underlying data have any time structure. For time series data, this implies
that the correlated signal of interest is observed with the same delay in each
channel and repetition. If this assumption is violated, drops in performance
will occur.

###  3.5 Effect of regularization on the reliability of EEG responses

A common numerical difficulty with generalized eigenvalue problems is the
required matrix inverse. In the case of CorrCA, it is necessary to invert
${\bf R}_{W}$ in equation (8). When some of the eigenvalues are small (or even
zero), then random fluctuations due to noise will dominate the inverse. It is
thus important to regularize the estimate of ${\bf R}_{W}$. In the past, we
have used the Truncated Singular Value Decomposition (TSVD) [45] and shrinkage
[46] for this purpose [1, 32, 33, 14]. These methods are described in detail
in Appendix C.

To evaluate the effects of these regularization methods on the results of
CorrCA we use data from our study on neural responses to advertisements [32,
we select one commercial among several aired during the 2012 Super Bowl]. For
both the TSVD and shrinkage approaches, we spanned the space of regularization
strength (a single parameter, $K$ and $\gamma$, in each case) and measured
ISC. The training set was the EEG data from 12 subjects during the first
viewing of the ad, while the EEG data collected during the second viewing
served as the test set.

The sum of the first 3 ISC values for both regularization methods followed a
similar trajectory with increasing regularization strength (decreasing $K$ and
increasing $\gamma$; Figure 6 B & C). On the training data, ISC drops
monotonically, whereas on the test data a maximum is achieved at intermediate
regularization values. On these data, shrinkage regularization was most
effective at $\gamma=0.4$ and TSVD when retaining $K=20$ out of a total of
$D=64$ dimensions in the within-subject covariance ${\bf R}_{W}$. The forward
models of the first three correlated components for each regularization
technique are shown in Figure 6 A. Notice that regularization tends to smooth
out the topographies of the second and third components compared to no
regularization.

Figure 6: Regularization increases reliability of EEG test data. The inter-
subject correlations of EEG responses from 12 subjects viewing a popular
television advertisement [32]. We took the data of the highest-scoring
advertisement maximized for ISC while varying levels of TSVD and shrinkage
regularization. (A) The spatial distributions of the three most correlated
components. Note that the distributions of the second and third components
were distinct for unregularized, TSVD-regularized and shrinkage-regularized
CorrCA. (B) Performance with TSVD as a function of $K$ (decreasing $K$
corresponds to stronger regularization), and (C) performance with shrinkage as
a function of $\gamma$ (increasing $\gamma$ corresponds to stronger
regularization). Components were extracted on the data from a first viewing of
the advertisement (training) and evaluated on data from a second viewing
(test). Reliability peaks at $\gamma=0.4$ for shrinkage, and when truncating
${\bf R}_{W}$ to $K=20$ dimensions, with very similar peak values of $\approx
0.16$ (sum of first three inter-subject correlations) for the two
regularization approaches.

##  4 Validation on simulated data

We evaluated the conditions under which CorrCA identifies the true source
directions. We also determined the performance of CorrCA at identifying
components with high ISC, and compared this against a simple benchmark based
on PCA of the subject-average data. Additionally, we compared the accuracy of
the statistical tests described in Appendix D at estimating the correct number
of components. To this end, we generated artificial data with known
properties. The data are generated assuming a known number of correlated
components (signals shared across subjects) and additive spatially correlated
noise not shared across subjects. The following factors were varied
parametrically: the number of subjects (raters, repeats), $N$, the number of
samples per subject, $T$, the number of dimensions, $D$, the number of
correlated components, $K$, and the signal-to-noise ratio (SNR) as defined in
the space of multivariate measurements. Each simulated correlated component
was identically reproduced in each subject, so that the ISC of all components
was equal to one. We further varied the temporal dependency structure of the
signal and noise components (IID or pink noise), the distribution of the
signal and noise components (Gaussian or $\chi^{2}$ distributed), and the
heterogeneity of the signal and noise subspaces (either the same or different
for all subjects). Signal components were generated to be correlated across
subjects, but uncorrelated with each other. Noise components were uncorrelated
with each other as well as across subjects. The number of noise components was
fixed to $D$. The following default parameters were used unless otherwise
noted: $N=5$, $T=200$, $D=30$, $K=10$, SNR = 0 dB, no regularization, IID
Gaussian signal and noise components, ISC = 1 for all $K$ components. Details
on data generation and evaluation methods can be found in Appendix E.

###  4.1 Accuracy as a function of Signal-to-Noise Ratio

Figure 7 depicts CorrCA performance as a function of measurement SNR for fixed
$T=200,D=30,N=5,K=10$. SNR values range from -40 dB to 40 dB. It is apparent
that average training and test ISC reach the optimal value of 1 for very high
SNR (panel A). For very low SNR ranges, training ISC significantly exceeds
test ISC, indicating over-fitting. The identification of the correlated
subspace converges to the optimal value for high SNR (panel B). However,
individual components are not perfectly identifiable even for high SNR because
they share the same ISC values (see Appendix A.3 for a discussion). The
identifiability of individual components becomes possible if the ISC of the
simulated signal components are adjusted to linearly decrease from 1 down to
1/$K$, and if the measurement SNR is high (data not shown). In practice, one
should keep in mind that the uniqueness of the ISC (and thereby, their
individual identifiability) depends on the amount of noise that is collinear
to each original correlated dimension. In the IID case, all estimates of $K$
converge to the true number of 10 for high SNR (panel C). Here, the parametric
test and the surrogate data approach using random circular shifts perform
similarly, while the surrogate data approach using phase scrambling performs
considerably worse at low SNR, presumably due to the absence of a well-defined
notion of frequency and phase in the sense of the Fourier transform. In the
case of dependent (pink noise) samples, both surrogate data approaches perform
identically, converging to the true value of $K$ for high SNR (panel D). In
contrast, the parametric test, which incorrectly assumes IID data,
overestimates the number of correlated components across the entire SNR range,
as expected.

Note that the significance test is applied here to multiple components
($D=30$). Under appropriate conditions, all three methods converge to the
correct number of dimensions in these simulate data ($K=10$). This indicates
that all three methods (described in Appendix D) properly account for these
multiple comparisons.

Figure 7: CorrCA performance as a function of signal-to-noise ratio (SNR).
Other parameters were fixed to $T=200$, $D=30$, $N=5$, $K=10$. Plotted are
averages across 100 experimental repetitions, while corresponding standard
deviations are shown as transparent areas. (A) The average ISC of the
reconstructed components converges to the maximal value of 1 for very high
SNR. For very low SNR, the expected value of 0 is exceeded on training data,
indicating overfitting. CorrCA outperforms a simple PCA on the subject-
averaged data. (B) CorrCA finds to correct subspace of correlated components
if the SNR is sufficiently high. Individual components may, however, not be
identifiable if their ISC are too similar. (C) In case of IID data, all
statistical tests to estimate $K$ converge to the true number of 10 for high
SNR. (D) In case of dependent (pink noise) samples, only surrogate data
approaches perform well, while the parametric test overestimates $K$.

###  4.2 Accuracy of estimating dimensionality of correlated components

The performance of CorrCA as a function of the true number of correlated
components, $K$, is depicted in Figure 8. Here, $T=200$, $D=30$, $N=5$, SNR =
0 dB, where $K$ was varied from 3 to 28. As the overlap between signal and
noise subspaces increases with $K$, the average ISC of the reconstructed
correlated components drops (panel A). A similar decline is seen in the
reconstruction of simulated correlated components and their forward models
(data not shown). As a result, the percentage of correlated components that
can be recovered (i.e., reach statistical significance) also decreases with
increasing $K$ (panels B & C).

We also compared the ISC values obtained with CorrCA against a simple
benchmark, namely, PCA of the subject-averaged data. While this simple method
captures the correct directions for high SNR case, CorrCA significantly
outperforms at medium and low SNR values (Figure 7 A). The same is true for
simulations with varying number of underlying dimensions (Figure 8 A), where
we see that CorrCA outperforms this simple method in all instances tested.

Figure 8:  CorrCA performance as a function of the true number of correlated
components, K. Other parameters were fixed to $T=200$, $D=30$, $N=5$, SNR = 0
dB. (A) The average ISC of the reconstructed correlated components drops with
increasing $K$. (B) Consequently, the number of correlated components can be
less well estimated, the higher $K$. For IID data, this number is bound from
above by the true number for all statistical approaches. (C) Overestimation of
$K$ may occur for pink noise data when the parametric F test is used.

Figure 9 shows the estimated number of correlated components as a function of
the number of data points, $T$, and the number of repetitions
(subjects/raters), $N$, with fixed $D=30$, $K=10$, SNR = 0 dB. In case of IID
data, all estimates converge to the true value of 10 for increasing $T$ or
$N$. This is also true for pink noise data in combination with non-parametric
statistical tests. For the parametric test, however, the number of seemingly
correlated components diverges as $T$ increases, and is overestimated
regardless of $N$.

Figure 9:  Estimated number of correlated components as a function of the
number of data points, T, (left panels), and the number of repetitions
(subjects/raters), N, (right panels). Other parameters were fixed to $D=30$,
$K=10$, SNR = 0 dB. (A, B) In case of IID data, all estimates converge to the
true value of 10 for increasing $T$ or $N$. (C, D) For pink noise data, only
non-parametric statistical tests converge to the true value.

###  4.3 Dependence on signal and noise distributions

Figure 10 A & B shows the performance of CorrCA for different signal and noise
distributions. Here, $T=200$, $D=30$, $N=2$, $K=10$, SNR = 0 dB. Signal and
noise components were generated either as Gaussian or squared Gaussian
($\chi^{2}$-distributed). For Gaussian distributed components, we considered
an additional case, where observations ${\bf x}^{l}(t)$ were dichotomized
(thus, Bernoulli distributed) by taking their sign. This simulation scenario
resembled the typical setting in which two raters provide a binary assessment
of the same issue. While non-Gaussianity of the underlying components does not
substantially affect the extraction of correlated components, the same is
considerably more difficult in the presence of dichotomous observations. An
additional analysis was conducted for IID Gaussian data and $N=5$ to assess
the performance of CorrCA when signal and noise covariance structures (as
defined by the mixing matrices ${\bf A}_{\text{s}}^{l}$ and ${\bf
A}_{\text{n}}^{l}$, $n=1,\ldots,N$) differ across subjects (Figure 10 panel C
& D). Significant drops in the number of estimated correlated components as
well as the average ISC are observed if the signal directions, the noise
directions, or both, differ across subjects. Note that varying directions of
the non-reproducible component ${\bf A}_{\text{n}}^{l}$ across subjects is
equivalent to adding noise, so that the dependence on the noise directions
shown in Figure 10 (panel C & D) is a restatement of the dependence on SNR of
Figure 7. The result illustrates that identical directions for the signal
components across subjects are a fundamental assumption of CorrCA. If this
assumption is suspected to be violated, the use of MCCA (section 2.6) may be
beneficial, or approaches that can gradually move from CorrCA to MCCA [47].

Figure 10:  CorrCA performance as a function of signal and noise
distributions. Other parameters were fixed to $T=200$, $D=30$, $K=10$, SNR = 0
dB. (A, B) Performance for Gaussian and squared Gaussian
($\chi^{2}$-distributed) signal and noise components, as well as for
dichotomized (Bernoulli-distributed) observations. Data points were either
drawn IID or sampled from a pink noise process. While non-Gaussianity of the
underlying components does not substantially affect the extraction of
correlated components, the same is considerably more difficult in the presence
of dichotomous observations. (C, D) performance on IID Gaussian data when
signal and noise covariance structure differs across repetitions (e.g.,
subjects, raters). Significant drops in the number of estimated correlated
components as well as the average ISC are observed if the signal covariance,
the noise covariance, or both, differ across subjects.

##  5 Conclusion

The goal of this work was to provide a formal, yet didactic and comprehensive
treatment of this new component analysis technique. The analytic development
resulted in a direct link to MCCA and a surprising equivalence to LDA, thus
linking the new methods with classic concepts of multivariate analysis. We
also found an efficient scheme for computing ISC/IRC which surprisingly does
not require computation of pairwise correlations. We also identified the
F-statistic as an exact parametric test for statistical significance, which is
valid in the case of normal and identically distributed samples. For the case
of dependent samples we demonstrated and validated circular-shift shuffle
statistics as an efficient non-parametric test of significance. We also
presented an extension to non-linear mappings using kernels. Finally, perhaps
the most important contribution of this work is the formalization of inter-
subject correlation. Conventional metrics, such as pairwise Pearson’s
correlation are limited to two signals. Inter-subject correlation as defined
here is applicable to more than two signals. ISC differs from conventional
definitions of correlation (see section refapdx:ICC), yet seamlessly
integrates in well established mathematics related to linear discriminants,
canonical correlation, and the F-statistic.

We make all code and data available at http://parralab.org/corrca, and for
long-term storage here [48]. This code is significantly faster than previous
instantiations and has been extensively evaluated. All figures presented in
this paper can be reproduced with this new code base and associated data. We
hope that this analysis and code inspires new uses for correlated component
analysis.

##  6 Acknowledgments

We are grateful to Valentin Rousson who pointed us to work on inter-class
correlation, and Oskar Jenni of the Children University Hospital in Zurich,
Switzerland who collected the data published in [37], which we reanalyzed here
in Figure 3. Similarly, we would like to thank Mike Milhalm for pointing us to
the problem of identifying parent/child agreement and Lindsay Alexander for
providing access to the Healthy Brain Network data [16] for Figure 4. We also
thank Samantha Cohen for providing the data from [14] for Figure 2. We thanks
Samantha Cohen, Jens Madsen Michael X Cohen, and two anonymous reviewers for
useful comments on an earlier version of this manuscript. We want to thank
Flaviano Morone for assisting with the proof in Appendix A.4.

##  Appendix A Generalized eigenvectors for two symmetric matrices

###  A.1 Joint diagonalization

Here we will review the well-know relationship between the joint
diagonalization of two symmetric matrices ${\bf A}$ and ${\bf B}$ and the
eigenvectors of ${\bf B}^{-1}{\bf A}$. Consider eigenvectors arranged as
columns in ${\bf V}$ with the corresponding eigenvalues in the diagonal matrix
${\bm{\Lambda}}$:

${\bf A}{\bf V}={\bf B}{\bf V}{\bm{\Lambda}}.$ (36)

To find the solution ${\bf V}$ for this equation, we can replace the positive
definite matrix ${\bf B}$ with its Cholesky factorization

${\bf B}={\bf L}{\bf L}^{\top}.$ (37)

Insert this into (36) and left-multiply with ${\bf L}^{-1}$ to obtain:

$\displaystyle({\bf L}^{-1}{\bf A}{\bf L}^{-\top})({\bf L}^{\top}{\bf V})$
$\displaystyle=$ $\displaystyle({\bf L}^{\top}{\bf V}){\bm{\Lambda}},$ (38)
$\displaystyle{\bf M}{\bf U}$ $\displaystyle=$ $\displaystyle{\bf
U}{\bm{\Lambda}},$ (39)

with ${\bf M}={\bf L}^{-1}{\bf A}{\bf L}^{-\top}$ and ${\bf U}={\bf
L}^{\top}{\bf V}$. Equation (39) is a conventional eigenvalue equation with
the same eigenvalues as (36). It is now easy to see that the corresponding
eigenvectors ${\bf V}={\bf L}^{-\top}{\bf U}$ diagonalize matrix ${\bf A}$

${\bf V}^{\top}{\bf A}{\bf V}={\bf U}{\bf L}^{-1}{\bf A}{\bf L}^{-\top}{\bf
U}={\bf U}^{\top}{\bf M}{\bf U}={\bm{\Lambda}}.$ (40)

They also diagonalize ${\bf B}$, which we see by left-multiplying (36) with
${\bf V}^{\top}$ and right-multiplying with ${\bm{\Lambda}}^{-1}$. This
yields:

${\bf V}^{\top}{\bf B}{\bf V}={\bf V}^{\top}{\bf A}{\bf
V}{\bm{\Lambda}}^{-1}={\bf I}.$ (41)

Thus, the eigenvectors that solve (36) jointly diagonalize ${\bf A}$ and ${\bf
B}$. The reverse is also true. A matrix ${\bf V}$ that satisfies:

$\displaystyle{\bf V}^{\top}{\bf A}{\bf V}$ $\displaystyle={\bm{\Lambda}}_{A}$
(42) $\displaystyle{\bf V}^{\top}{\bf B}{\bf V}$
$\displaystyle={\bm{\Lambda}}_{B},$ (43)

where ${\bm{\Lambda}}_{A}$ and ${\bm{\Lambda}}_{B}$ are diagonal matrices,
also satisfy (36). To see this, solve (43) for ${\bf V}^{\top}$ and substitute
this into (42), which yields:

${\bm{\Lambda}}_{B}{\bf V}^{-1}{\bf B}^{-1}{\bf A}{\bf V}={\bm{\Lambda}}_{A},$
(44)

and therefore

${\bf B}^{-1}{\bf A}{\bf V}={\bf V}{\bm{\Lambda}}_{B}^{-1}{\bm{\Lambda}}_{A},$
(45)

which is the same as the eigenvalue equation (36) with

${\bm{\Lambda}}={\bm{\Lambda}}_{B}^{-1}{\bm{\Lambda}}_{A}.$ (46)

###  A.2 Linear combinations

A matrix ${\bf V}$ that simultaneously diagonalizes two matrices ${\bf A}$ and
${\bf B}$, as in (42)-(43), will also diagonalize any linear combination of
these two matrices [22], because

$\displaystyle{\bf V}^{\top}(\alpha{\bf A}+\beta{\bf B}){\bf V}$
$\displaystyle=$ $\displaystyle\alpha{\bf V}^{\top}{\bf A}{\bf V}+\beta{\bf
V}^{\top}{\bf B}{\bf V}$ (47) $\displaystyle=$
$\displaystyle\alpha{\bm{\Lambda}}_{A}+\beta{\bm{\Lambda}}_{B},$

which is obviously a diagonal matrix. Therefore, the eigenvectors of ${\bf
A}^{-1}{\bf B}$ are also eigenvectors of $(\alpha{\bf A}+\beta{\bf
B})^{-1}(\gamma{\bf A}+\delta{\bf B})$, assuming vectors $[\alpha,\beta]$ and
$[\gamma,\delta]$ are not collinear.

###  A.3 Identifiability

An additional important observation is that the eigenvectors of (36) are only
uniquely defined for unique eigenvalues. If an eigenvalue is degenerate, say
with multiplicity $K$, then the corresponding eigenvectors span a
$K$-dimensional subspace. Any vector within that subspace is a solution to the
eigenvalue equation with that eigenvalue. In this space any arbitrary rotation
will yield a valid solutions. Thus, when estimating components, it is possible
that solutions with close-by ISC are “mixed”. Such components are not uniquely
identifiable. One should also note that the sign of any eigenvector is
arbitrary. If ${\bf v}_{d}$ is a solution, then so is $-{\bf v}_{d}$. One can
always arbitrarily swap the sign of ${\bf v}_{d}$ for any component $d$, which
also reverses the sign of $y_{d}^{l}$ and ${\bf a}_{d}$ — the $d$th column of
the forward model ${\bf A}$. Thus, keep in mind that in all figures the sign
of individual ${\bf v}_{d}$ and ${\bf a}_{d}$ are arbitrary. The same problem
arises with most signal decomposition methods such as PCA, ICA or CCA. A
common approach used in practice to “fix” the sign is to select one
sensor/dimension with the strongest forward-model coefficient, and set that to
a positive value (if indeed it is significantly larger than other coefficients
with opposite sign).

###  A.4 Optimization of sum ISC

The optimization criterion (9) has the form

$J({\bf V})=\sum_{d}\frac{{\bf v}_{d}^{\top}{\bf A}{\bf v}_{d}}{{\bf
v}_{d}^{\top}{\bf B}{\bf v}_{d}}$ (48)

and is subject to the constraints that components are uncorrelated (relative
to ${\bf B}$):

${\bf v}_{d}^{\top}{\bf B}{\bf v}_{c}=0,\text{for }d\neq c.$ (49)

It will be convenient to convert this optimality criterion into conventional
Rayleigh quotients. To do so use the Cholesky decomposition (37) as before and
replace, ${\bf v}_{d}={\bf L}^{-\top}{\bf w}_{d}$ and ${\bf A}={\bf L}{\bf
M}{\bf L}^{\top}$ to yield:

$J({\bf W})=\sum_{d}\frac{{\bf w}_{d}^{\top}{\bf M}{\bf w}_{d}}{{\bf
w}_{d}^{\top}{\bf w}_{d}}$ (50)

with the constraint now being a conventional orthogonality constraint:

${\bf w}_{d}^{\top}{\bf w}_{c}=0,\text{for }d\neq c.$ (51)

One can now optimize for matrix ${\bf W}$ with the ${\bf w}_{d}$ as its column
vectors. The solution to this constrained optimization problem is well
established [49, Corollary 4.3.39]. It is given by, ${\bf W}={\bf U}$, the
eigenvalue matrix of equation (39). Inserting ${\bf U}={\bf L}{\bf V}$ into
(39) gives the solution to the original problem, which is given by the
generalize value equation (36).

To see that the eigenvector matrix of ${\bf M}$ maximizes $J({\bf W})$, first
note that except for the constraint, the terms in the sum can be optimized
independently, and that they all have the same form. One can therefore
consider the maximum of a single term, written here simply as

$\rho=\frac{{\bf w}^{\top}{\bf M}{\bf w}}{{\bf w}^{\top}{\bf w}},$ (52)

which is to be optimized subject to the orthogonality constraint (51). From
here on the argument proceeds conventionally. Let ${\bf u}_{i},i=1\ldots D$ be
the orthonormal eigenvectors of matrix ${\bf M}$ with eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{D}$. We can express the
solution vectors in this orthogonal basis as

${\bf w}=\sum_{i=1}^{D}\alpha_{i}{\bf u}_{i},$ (53)

where $\alpha_{i}$ is the projection of ${\bf w}$ on that orthonomal basis of
these eigenvectors, $\alpha_{i}={\bf u}_{i}^{T}{\bf w}$. Inserting (53) into
(52) gives

$\rho=\frac{\sum_{i}\sum_{j}\alpha_{i}{\bf u}_{i}^{\top}{\bf M}{\bf
u}_{j}\alpha_{j}}{\sum_{i}\sum_{j}\alpha_{i}{\bf u}_{i}^{\top}{\bf
u}_{j}\alpha_{j}}=\frac{\sum_{i}\lambda_{i}\alpha_{i}^{2}}{\sum_{i}\alpha_{i}^{2}}\leq\frac{\lambda_{1}\sum_{i}\alpha_{i}^{2}}{\sum_{i}\alpha_{i}^{2}}=\lambda_{1}.$
(54)

This upper bound is attained for $\alpha_{1}>0$ and all others equal zero, and
thus ${\bf w}_{1}={\bf u}_{1}$. For the second vector we require that it is
orthogonal to the first, and thus we have to constrain $\alpha_{1}=0$. The sum
in (54) now excludes the first dimension and has therefore an upper bound of
$\lambda_{2}$ which is attained for ${\bf w}_{2}={\bf u}_{2}$. The argument
continues iteratively until for the last component only, ${\bf w}_{D}={\bf
u}_{D}$, is orthogonal to all other projections. It also happens to be the
minimum of $\rho$.

###  A.5 Uncorrelated components vs orthogonal projections

We have required here that the extracted components are uncorrelated, i.e. the
projection vectors diagonalize ${\bf R}_{W}$: ${\bf V}^{\top}{\bf R}_{W}{\bf
V}={\bf I}$. This means that the component dimensions of ${\bf y}_{i}^{l}$ are
uncorrelated, but only in the averaged over all subject/repeats $l$. Note that
for an individual subject, ${\bf y}_{i}^{l}$ may not be uncorrelated. The same
is true for MCCA as derived in section 2.6. Instead of uncorrelated, we could
have required that the projection vectors are orthogonal, ${\bf V}^{\top}{\bf
V}={\bf I}$, i.e. the linear transformation is limited to rotations. To find
the orthogonal (rotations) that maximize ISC, one will have to resort to
constrained optimization algorithms [e.g. 21]. A related issues arises in the
context of Canonical Correlation Analysis which attempts to find a linear
mapping between two datasets such that the components are uncorrelated [50,
51]. Alternatively, the Procrustes problem aims to find a mapping between two
data-sets with orthogonal projections [52]. [The later is what is used for
“hyper-alignment” between subjects in fMRI by 53]. CCA and Procrustes are
related, and in fact one can smoothly titrate between orthogonal and
uncorrelated solutions [54]. In our case we favor uncorrelated components,
because we want the components to capture independent information. At a
minimum, this means the components have to be uncorrelated. In that sense, our
method can be seen as a source separation method, similar to Independent
Component Analysis [41], Denoising Source Separation [42], or others [55]. One
advantage of restricting the transformations to rotations is that euclidean
distances in the new space are preserved. This means that differences found in
the original sensor space are preserved in the new component space. Another
advantage is that all original dimensions contribute equally to the new
dimensions. In contrast, the non-orthogonal shearing operation of CorrCA or
CCA may discount some dimensions and emphasize others. Both these advantages
of orthonormal mappings, however, require that _all_ new dimensions are
considered. This is typically the case in the applications of hyper-alignment
[56]. For the purposes of identifying dimensions that are reliably reproduced,
however, we find in practice that many dimensions are not significantly
correlated between repeats, and one will typically want to remove those.
Similarly, one will want to dismiss sensors in the original data that are not
reproducible across repeats while emphasizing others. Ultimately, however, the
application will dictate whether one should enforce uncorrelated components or
orthogonal projections. In either case, the goal can remain to maximize ISC,
and the solution can be found, either through an eigenvalue question (as we
presented here for uncorrelated components) of with constrained optimization
[for orthogonal projections, along the lines of 21].

##  Appendix B Relationship between Scatter matrices and between- and within-
subject covariance matrices

Comparing Eq. (16) to (14) yields:

${\bf R}_{T}=N{\bf S}_{B}.$ (55)

One can also relate ${\bf R}_{W}$ to ${\bf S}_{T}$ as follows. Expanding the
definition (6) yields (with abbreviated notation for the sums):

$\displaystyle{\bf R}_{W}$ $\displaystyle=$ $\displaystyle\sum_{il}{\bf
x}_{i}^{l}{\bf x}_{i}^{l\top}-\sum_{il}{\bf x}_{i}^{l}{\bar{\bf
x}}_{*}^{l\top}-\sum_{il}{\bar{\bf x}}_{*}^{l}{\bf
x}_{i}^{l\top}+\sum_{il}{\bar{\bf x}}_{*}^{l}{\bar{\bf x}}_{*}^{l\top}$ (56)
$\displaystyle=$ $\displaystyle\sum_{il}{\bf x}_{i}^{l}{\bf
x}_{i}^{l\top}-T\sum_{l}{\bar{\bf x}}_{*}^{l}{\bar{\bf x}}_{*}^{l\top}$
$\displaystyle=$ $\displaystyle\sum_{il}{\bf x}_{i}^{l}{\bf
x}_{i}^{l\top}-T\sum_{l}{\bar{\bf x}}_{*}^{l}{\bar{\bf
x}}_{*}^{l\top}+T\sum_{l}{\bar{\bf x}}_{*}^{*}{\bar{\bf
x}}_{*}^{*\top}-T\sum_{l}{\bar{\bf x}}_{*}^{*}{\bar{\bf x}}_{*}^{*\top}$
$\displaystyle=$ $\displaystyle\left(\sum_{il}{\bf x}_{i}^{l}{\bf
x}_{i}^{l\top}+T\sum_{l}{\bar{\bf x}}_{*}^{*}{\bar{\bf
x}}_{*}^{*\top}\right)-\left(\sum_{il}{\bar{\bf x}}_{*}^{l}{\bar{\bf
x}}_{*}^{l\top}-T\sum_{l}{\bar{\bf x}}_{*}^{*}{\bar{\bf
x}}_{*}^{*\top}\right)$ $\displaystyle=$ $\displaystyle\sum_{il}({\bf
x}_{i}^{l}-{\bar{\bf x}}_{*}^{*})({\bf x}_{i}^{l}-{\bar{\bf
x}}_{*}^{*})^{\top}-\sum_{il}({\bar{\bf x}}_{*}^{l}-{\bar{\bf
x}}_{*}^{*})({\bar{\bf x}}_{*}^{l}-{\bar{\bf x}}_{*}^{*})^{\top}$
$\displaystyle=$ $\displaystyle{\bf S}_{T}-T{\bf S}_{M}.$

Matrix ${\bf S}_{M}$, defined in (21), captures the variance of the mean
across subjects and is zero when all subjects have equal mean: ${\bar{\bf
x}}_{*}^{l}={\bar{\bf x}}_{*}^{*}$ .111 For time domain signals, this is often
the case as it is customary to subtract the mean value (constant offset) to
make the signals zero-mean. For rating data, an unequal mean indicates rater
bias, i.e. a rater gives systematically lower/higher ratings than other
raters, and it may be advantageous to remove this bias by equalizing all mean
values. Note that, when ${\bf S}_{M}=0$, the present definition of ISC is
equivalent to the classic definition of intra-class correlation (ICC, see
Appendix B). But in general, this matrix is not zero, which is important when
generalizing the approach to non-linear transformations. There, enforcing zero
mean in ${\bf x}$ does not guarantee zero mean in ${\bf y}$. With this result,
scatter and covariance matrices can be related as follows:

$\displaystyle{\bf S}_{B}$ $\displaystyle=$ $\displaystyle\frac{1}{N}{\bf
R}_{T}=\frac{1}{N}({\bf R}_{W}+{\bf R}_{B}),$ (57) $\displaystyle{\bf S}_{W}$
$\displaystyle=$ $\displaystyle{\bf S}_{T}-{\bf S}_{B}={\bf R}_{W}+T{\bf
S}_{M}-{\bf S}_{B}$ (58) $\displaystyle=\frac{N-1}{N}{\bf
R}_{W}-\frac{1}{N}{\bf R}_{B}+T{\bf S}_{M}.$

##  Appendix C Regularization and robustness

Note that a high inter-subject correlation corresponds to large eigenvalues of
$\mathbf{R}_{W}^{-1}\mathbf{R}_{B}$. These eigenvalues will be large if the
selected direction exhibits large power in the space of $\mathbf{R}_{B}$, or
_low_ power in the space of $\mathbf{R}_{W}$. The former is meaningful and
desired, but the latter can occur with rank-deficient or impoverished data and
may lead to spuriously high inter-subject correlation. Truncating the
eigenvalue spectrum of $\mathbf{R}_{W}$ has the desired effect of forcing the
extracted data dimensions to have both high power in $\mathbf{R}_{B}$ _and_
$\mathbf{R}_{W}$. This makes it less likely to find spuriously reliable
dimensions. Formally, the TSVD approach is to select $\mathbf{V}$ according
to:

$\displaystyle{\bf R}_{B}{\bf V}=\tilde{{\bf R}}_{W}{\bf V}{\bm{\Lambda}},$
(59)

where the within-subject covariance is now regularized with the following
approximation:

$\displaystyle\tilde{{\bf R}}_{W}=\tilde{{\bf
U}}_{W}\tilde{{\bm{\Lambda}}}_{W}\tilde{{\bf U}}_{W}^{\top},$ (60)

where $\tilde{\bf U}_{W}=\left[{{\bf u}}_{1}\ldots{{\bf u}}_{K}\right]$, and
$\tilde{\bm{\Lambda}}_{W}=\mathtt{diag}\left[\lambda_{1}\ldots\lambda_{K}\right]$
are the $K$ principal eigenvectors of ${\bf R}_{W}$ and associated
eigenvalues. The inverse of ${\bf R}_{W}$ is computed according to:

$\displaystyle\tilde{{\bf R}}_{W}^{-1}=\tilde{{\bf
U}}_{W}\tilde{{\bm{\Lambda}}}_{W}^{-1}\tilde{{\bf U}}_{W}^{\top},$ (61)

which then left-multiplies (59) and yields the regularized solution to the
CorrCA problem.

The parameter $K$ is the number of eigenvectors to retain in the
representation of ${\bf R}_{W}$. Decreasing this number strengthens the level
of regularization, and $K=D$ is equivalent to not regularizing.

Shrinkage regularization operates by a similar principle: dimensions of the
data exhibiting low variance (corresponding to low eigenvalues of ${\bf
R}_{W}$) are enriched by incrementing these eigenvalues:

$\displaystyle\tilde{{\bf R}}_{W}={{\bf U}}_{W}\tilde{\bm{\Lambda}}_{W}{\bf
U}_{W}^{\top},$ (62)

where
$\tilde{\bm{\Lambda}}_{W}=(1-\gamma){\bm{\Lambda}}_{W}+\gamma\bar{\lambda}{\bf
I},$ $\gamma$ is the shrinkage parameter with $0\leq\gamma\leq 1$ , and
$\bar{\lambda}=\mathrm{Tr}({\bf R}_{W})/D$ is the mean eigenvalue of the
(unregularized) within-subjects covariance. The effect of shrinkage, which
retains the full dimensionality of the covariance matrix, is that the smaller
eigenvalues are increased and larger eigenvalues are reduced (hence the term
“shrinkage”).

The advantage of shrinkage regularization is the simplicity of its
implementation, where the regularized covariance can be simply computed as

$\tilde{{\bf R}}_{W}=(1-\gamma){\bf R}_{W}+\gamma\bar{\lambda}{\bf I}$ (63)

without any further modifications to the generalized eigenvalue routine. TSVD
instead requires more careful handling of the null-space when computing the
general eigenvalues (see provided code). On the other hand, TSVD has the
important advantage that the solutions ${\bf V}$ to the regularized eigenvalue
problem (59) diagonalize $\tilde{\bf R}_{W}$ as well as the original ${\bf
R}_{W}$. This means that the extracted components are strictly uncorrelated,
whereas for shrinkage that is only approximately the case.

The performance of TSVD and shrinkage can be evaluated by maximizing ISC on
one data set (training set) and using the resulting projection vectors to
measure ISC on a different data set (test set). We do this on one of the
application examples in section 3.5.

Note that the TSVD approach can also be used to regularize the MCCA algorithm.
In this case, the regularized inverse $\tilde{\bf D}$ has to be computed
separately for each block of the block-diagonal matrix ${\bf D}$ in (34).

##  Appendix D Testing for statistical significance

The assessment of statistical significance of $\rho$ requires knowledge of
their distribution under the null hypothesis. The null hypothesis associated
with CorrCA states that there exists no shared one-dimensional subspace in
which the data of at least two subjects are correlated. Estimating the null
distribution is complicated by two factors. First, CorrCA has free parameters,
the optimization of which inflates the ISC values under the null hypothesis to
an extent governed by the dimensionality of the problem (parameters $T$, $D$,
and $N$). Second, temporal correlations in the data lead to larger variances
than expected under the assumption of independent and identically distributed
(IID) data typical for parametric statistical tests, and, thereby, to false
positive detections of significant ISC values. The IID assumption may be
reasonable for rating data, where rated individuals may be independent from
one another, but it is not fulfilled for time series data such as EEG
recordings. Here, we present three approaches to estimate the significance of
individual correlated components as well as the dimensionality of the
correlated subspace for IID as well as auto-correlated data.

Parametric null distribution—For IID data, we can use the F-statistics given
in Eq. (26) to test for significant ISC (provided the means are equalized).
Evidently, this statistic is not applicable for $\rho$ values obtained on
training data, as the we have adapted parameters to optimize this statistic
itself. However, we can use the optimal projection vectors to assess
statistical significance on separate, previously unseen test data. An
important limitation is that the $T$ exemplars have to be independent. This is
true perhaps for rating data, but typically not true for temporal signals,
where samples in time are often correlated. For such data, the non-parametric
test discussed above are more appropriate. To avoid overestimating the number
of components due to multiple testing (there are always $D$ components to be
tested), we adopt a Bonferroni correction. This limits the family-wise error
rate, i.e. the probability of making one or more false discoveries among the
$D$ tests. Thus, only $S$ values with $p<\alpha/D$ are considered significant,
where $\alpha$ is the desired significance level and $p$ is computed using the
F-statistics.

Phase-scrambled surrogate data— For auto-correlated data we have to use
shuffle statistics. This approach was introduced by [57] as a tool to test for
nonlinearities in time series data. The idea is to create data that are
consistent with the null hypothesis but otherwise resemble the observed
experimental data in terms of temporal and spatial correlations. The original
surrogate data of [57] preserves the amplitude spectrum of the original data
but uses a random phase spectrum, which is achieved through an application of
the Fourier transform and its inverse. This approach was extended to
additionally maintain the spatial covariance structure of multivariate data by
using identical random phase shifts for all variables by [58]. Here, we use
the approach of Prichard and Theiler to test for significant $\rho$ based on
the consideration that no inter-subject correlation can be present after
phase-scrambling the data of each subject. Variants of this approach have
previously been used to test for significant correlations [59, 33, 14, 60].
The p-value for the ISC of a given correlated component is defined as the
frequency with which the ISC computed from the random-phase surrogate data is
exceeded by the maximum ISC of the most strongly correlated component in the
original data. Note that for every shuffle we obtain a set of ISC values,
$\rho_{d}$, for the $d$th component. When we compare the ISC value on the
original data, we measure p-value based on the distribution of maximum $\rho$,
i.e. all $\rho_{d}$ values are compared against the Null distribution of
$\rho_{1}$. With this approach no further correction for multiple comparison
is needed.

Random circular shifts—A simpler approach to create surrogate data for auto-
correlated data is to circularly shift the samples of each subject by a
different offset, where the same offset is used for all variables of the same
subject (i.e. sample index $i=1\ldots T$ are shifted by a random offset $o$:
$i\leftarrow(i+o)\text{mod }T+1$). Just as the phase-scrambling approach, this
procedure ensures that spatial correlations and spectral properties of the
original data are maintained. Additionally, non-linear properties of the
original data are maintained. This approach has previously been used in [1,
32]. P-values are defined analogous to the phase-scrambling approach.

We test the validity of the three methods on simulated data in section 4. The
goal is to determine under which conditions these methods determine the
correct number of underlying correlated components.

##  Appendix E Validation on simulated data: Data generation and evaluation
methods

The generated signal and noise components, ${\bf s}^{l}(t)\in\mathbb{R}^{K}$
and ${\bf n}^{l}(t)\in\mathbb{R}^{D}$, $t=1,\ldots,T$, $l=1,\ldots,N$, were
mapped to the measurement space as ${\bf x}_{\text{s}}^{l}(t)={\bf
A}_{\text{s}}^{l}{\bf s}^{l}(t)$, ${\bf x}_{\text{n}}^{l}(t)={\bf
A}_{\text{n}}^{l}{\bf n}^{l}(t)$ and normalized as ${\bf
x}_{\text{s}}^{l}(t)\leftarrow{\bf x}_{\text{s}}^{l}(t)/||{\bf
x}_{\text{s}}^{l}||_{F}$, ${\bf x}_{\text{n}}^{l}(t)\leftarrow{\bf
x}_{\text{n}}^{l}(t)/||{\bf x}_{\text{n}}^{l}||_{F}$, where $||{\bf x}||_{F}$
denotes the Frobenius norm of the multivariate data set ${\bf x}(t)$. The
mixing matrices ${\bf A}_{\text{s}}^{l}={\bf O}_{\text{s}}^{l}{\bf
D}_{\text{s}}^{l}\in\mathbb{R}^{D\times K}$ and ${\bf A}_{\text{n}}^{l}={\bf
O}_{\text{n}}^{l}{\bf D}_{\text{n}}^{l}\in\mathbb{R}^{D\times D}$ were
generated as products of matrices ${\bf
O}_{\text{s}}^{l}\in\mathbb{R}^{D\times K}$ (${\bf
O}_{\text{n}}^{l}\in\mathbb{R}^{D\times D}$) with random orthonormal columns
and diagonal matrices ${\bf D}_{\text{s}}^{l}\in\mathbb{R}^{K\times K}$ (${\bf
D}_{\text{n}}^{l}\in\mathbb{R}^{D\times D}$), the diagonal entries of which
were drawn randomly as ${{D}}_{ii}=\exp(d_{i})$, $d_{i}\sim{\mathcal{N}}(0,1)$
and normalized to a maximum of 1. This ensured that the non-zero eigenvalues
of the signal and noise covariances matrices decayed exponentially as is the
case for many real-world data sets. Unless otherwise noted, the same mixing
matrices ${\bf A}_{\text{s}}^{l}\equiv{\bf A}_{\text{s}}$ and ${\bf
A}_{\text{n}}^{l}\equiv{\bf A}_{\text{n}}$ were used for all subjects,
reflecting the fundamental assumption of CorrCA that the spatial covariance
structure of signal and noise components is the same in all subjects. The
emulated data were generated as ${\bf x}^{l}(t)=\xi{\bf
x}_{\text{s}}^{l}(t)+(1-\xi){\bf x}_{\text{n}}^{l}(t)$, where the scalar
parameter
$\xi~{}=~{}10^{\nicefrac{{\text{SNR}}}{{20}}}/(1+10^{\nicefrac{{\text{SNR}}}{{20}}})$
was used to adjust the SNR. Note that SNR was defined here in terms of the
contribution of the signal and noise components to the Frobenius norm of the
entire multivariate measurement. It is, therefore, different from the
definition provided in Eqs. (18) and (25), which applies to single channels or
components and relates to the ISC through Eq. (22). Note that in most of our
experiments, the ISC of each simulated signal component was set to 1,
corresponding to infinite SNR according to Eq. (22).

We applied CorrCA to the emulated data and measured its performance in terms
of the achieved ISC (average over the first $K$ components) on the training
set. We also measured the average ISC for a test set, i.e., using the same
CorrCA projections ${\bf V}$ on new emulated data of identical size. Moreover,
we measured the performance of CorrCA to reconstruct the true (emulated)
correlated components ${\bf s}^{l}(t)$ and their corresponding mixing matrices
${\bf A}_{\text{s}}^{l}$ by comparing them to ${\bf y}^{l}(t)$ and the
estimated forward model ${\bf A}$ (section 2.5). The performance metric for
this was the Pearson correlation between individual simulated and
reconstructed components or columns of ${\bf A}_{\text{s}}$. With this, we
measure if each component was correctly identified. However, because
components with equal ISC can be arbitrarily mixed within a subspace (Appendix
A.3), we also assessed identifiability in terms of the angle between the
subspaces spanned by the first $K$ columns of the estimated $\hat{\bf A}$ and
${\bf A}_{\text{n}}$, as well as the angle between the simulated components
${\bf s}^{l}(t)$ and the first $K$ dimensions of ${\bf y}^{l}(t)$
reconstructed by CorrCA. Angles were normalized to the interval $[0,1]$.
Finally, we estimated the number of simulated correlated components, $K$,
using phase-scrambled surrogate data, surrogate data obtained using random
circular shifts, and the parametric F-test. Empirical null distributions were
obtained using 1000 surrogate data sets. The number of components with
p-values smaller than $\alpha=0.05$ was used as an estimate for $K$. For the
parametric test, the $T$ samples were randomly split into training and test
sets of equal size, where the CorrCA projections ${\bf V}$ were obtained on
the training set, and the statistical assessment was conducted on the test
set. For this approach, the median number of components with p-values smaller
than $\alpha=0.05$ across 100 random train/test splits was used as an estimate
of $K$.

All experiments were repeated 100 times. Mean values and standard deviations
of all performance measures are depicted in Figures 7—10.

##  Appendix F Extension to non-linear maps using Kernels

So far we have only considered linear transformations between ${\bf x}$ and
${\bf y}$. Now, given the close relationship between CorrCA and LDA, we can
extend CorrCA to non-linear transformations following the approach of kernel-
LDA [61, 62, 63, 64]. Assume a non-linear mapping ${\bm{\Phi}}(\cdot)$ for the
data of subject $l$:

${\bf Z}^{l}={\bm{\Phi}}({\bf X}^{l}).$ (64)

In this notation, ${\bf X}^{l}$ has dimension $D\times T$ spanning the entire
multivariate time series for subject $l$, and ${\bf Z}^{l}$ has dimensions
$D^{\prime}\times T^{\prime}$. Thus, the non-linear transformation
${\bm{\Phi}}(\cdot)$ is a mapping from space-time to a new space with possibly
quite different dimensions. 222In the example of brain signals, ${\bf Z}^{l}$
could represent, for instance, the power spectrum in source space [65], where
$D^{\prime}$ captures the number of sources and $T^{\prime}$ the number of
frequency bins. In the application examples and in the code we simplify this
by applying the non-linear mapping to each time point separately, so that
$T=T^{\prime}$ and we have a purely spatial nonlinear mapping, similar to
conventional kernel-LDA. However, this more general formulation allows one in
principle to reduce the time dimension as well. This is important because the
computational complexity of the kernel approach scales with $O(T^{3})$, and
thus becomes prohibitive for long time-domain signals. The approach now is to
first extract non-linear features of the data with ${\bm{\Phi}}(\cdot)$, and
then combine these features linearly with projection ${\bf v}$:

${\bf y}^{l}={\bf v}^{\top}{\bf Z}^{l}.$ (65)

Note that the transformation ${\bm{\Phi}}(\cdot)$ does not need to be defined
explicitly. Instead, the “kernel trick” [66] allows one to specify this non-
linearity only implicitly, by defining instead the inner product of vectors in
this new $D^{\prime}$-dimensional feature space with a kernel function ${\bf
K}(\cdot,\cdot)$:

${{\bf Z}^{k}}^{\top}{\bf Z}^{l}={\bf K}({\bf X}^{k},{\bf X}^{l}).$ (66)

Knowing how to compute the inner product is sufficient, if all the expressions
in the algorithm can be formulated in terms of this inner product.333Note that
this inner product is only over dimension $D^{\prime}$. The function ${\bf
K}(\cdot,\cdot)$ is therefore only a ”kernel” in this $D^{\prime}$ dimensional
space. In a slight abuse of notation, we will write ${\bf K}^{kl}={\bf K}({\bf
X}^{k},{\bf X}^{l})$ as shorthand for this matrix defined for each pair of
subjects $lk$. Note that the dimensions of ${\bf K}^{kl}$ are
$T^{\prime}\times T^{\prime}$ and thus this matrix can potentially be quite
large. In cases where the non-linear transformation ${\bm{\Phi}}(\cdot)$ can
be expressed explicitly, and $D^{\prime}\ll T^{\prime}$, one may chose to
operate directly with ${\bf Z}^{l}$. Otherwise, it will be beneficial to
leverage the definition of the non-linear transformation in terms of ${\bf
K}^{kl}$. In the following, we will rewrite the optimization criterion $\rho$
in terms of ${\bf K}^{kl}$ alone.

The crucial step is to define the projection vectors in terms of the samples
of the mean across subjects.

${\bf v}=\sum_{i=1}^{T^{\prime}}\alpha_{i}{\bar{\bf z}}^{*}_{i}={\bar{\bf
Z}}^{*}\bm{\alpha},$ (67)

where ${\bar{\bf z}}^{*}_{i}$ is the across-subject mean vector for exemplar
$i$ in the non-linear mapped space ($i$th column of $\bar{\bf Z}^{*}$). Here,
$\alpha_{i}$ are parameters that indicate how exemplars $\bar{\bf z}_{i}$ are
to be combined to represent the projection vector ${\bf v}$. The approach will
be to find optimal vector $\bm{\alpha}$ instead of vector ${\bf v}$, and thus
we need to rewrite the optimality criterion in term of $\bm{\alpha}$.444As an
alternative to (67), we could have also defined the projection vector in terms
of the samples of all subjects:  ${\bf
v}=\sum_{l=1}^{N}\sum_{i=1}^{T^{\prime}}\alpha_{i}^{l}{\bf z}^{l}_{i}.$ (68)
The derivations can all be repeated analogously and the resulting algorithm is
the same, but with $NT^{\prime}$ degrees for freedom in $\bm{\alpha}$. In the
provided code, this is referred to as the “full model”, whereas equations
derived here for the expansion in terms of the mean (67) are referred to as
the “mean model”. In the example provided in Figure 11, the numerical results
are very similar for both models. The question as to whether the mean model is
sufficient is empirical and comes down to whether there is enough diversity in
the mean response to capture the individual variations across subjects.
Mathematically, this implies that the expansion in terms of the mean is
complete, i.e., the mean response matrix ${\bar{\bf Z}^{*}}$ is full rank.

Combining (65) with (67) one can now express the non-linear components ${\bf
y}$ as a linear combination of ${\bf K}^{lk}$ as follows:

$\displaystyle{\bf y}^{l}$ $\displaystyle=$ $\displaystyle{\bf v}^{\top}{\bf
Z}^{l}=\bm{\alpha}^{\top}\frac{1}{N}\sum_{k=1}^{N}{\bf Z}^{k^{\top}}{\bf
Z}^{l}$ (69) $\displaystyle=$
$\displaystyle\bm{\alpha}^{\top}\frac{1}{N}\sum_{k=1}^{N}{\bf
K}^{kl}=\bm{\alpha}^{\top}{\bar{\bf K}}^{*l},$

where the bar in ${\bar{\bf K}}^{*l}$, indicates that we are taking the mean,
and the asterisk specifies over which index this mean is taken. It analogously
follows that,

${\bf v}^{\top}{\bar{\bf Z}}^{*}=\bm{\alpha}^{\top}{\bar{\bf K}}^{**}.$ (70)

Denote the columns of matrix ${\bf Z}^{l}$ and ${\bar{\bf K}}^{*l}$ as ${\bf
z}_{i}^{l}$ and ${\bar{\bf k}}_{i}^{*l}$. Then we can also write with this
notation for the average:

${\bf v}^{\top}{\bar{\bf z}}_{i}^{*}=\bm{\alpha}^{\top}{\bar{\bf
k}}_{i}^{**}.$ (71)

The within-subject and total variance of the non-linear features ${\bf
z}_{i}^{l}$ can thus be written as:

$\displaystyle r_{W}$ $\displaystyle=$ $\displaystyle{\bf v}^{\top}{\bf
R}_{W}{\bf v}=\sum_{l=1}^{N}\sum_{i=1}^{T^{\prime}}{\bf v}^{\top}({\bf
z}_{i}^{l}-{\bar{\bf z}}_{*}^{l})({\bf z}_{i}^{l}-{\bar{\bf
z}}_{*}^{l})^{\top}{\bf v}$ (72) $\displaystyle=$
$\displaystyle\bm{\alpha}^{\top}{\bf C}_{W}\bm{\alpha},$ $\displaystyle r_{T}$
$\displaystyle=$ $\displaystyle{\bf v}^{\top}{\bf R}_{T}{\bf
v}=N^{2}\sum_{i=1}^{T^{\prime}}{\bf v}^{\top}({\bar{\bf z}}_{i}^{*}-{\bar{\bf
z}}_{*}^{*})({\bar{\bf z}}_{i}^{*}-{\bar{\bf z}}_{*}^{*})^{\top}{\bf v}$ (73)
$\displaystyle=$ $\displaystyle\bm{\alpha}^{\top}{\bf C}_{T}\bm{\alpha},$

where we defined the within-subject and total covariance of ${\bf k}_{i}^{*l}$
as:

$\displaystyle{\bf C}_{W}$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\sum_{i=1}^{T^{\prime}}({\bar{\bf
k}}_{i}^{*l}-{\bar{\bf k}}_{*}^{*l})({\bar{\bf k}}_{i}^{*l}-{\bar{\bf
k}}_{*}^{*l})^{\top},$ (74) $\displaystyle{\bf C}_{T}$ $\displaystyle=$
$\displaystyle N^{2}\sum_{i=1}^{T^{\prime}}({\bar{\bf k}}_{i}^{**}-{\bar{\bf
k}}_{*}^{**})({\bar{\bf k}}_{i}^{**}-{\bar{\bf k}}_{*}^{**})^{\top}.$ (75)

Because of the symmetry of these expressions with the definitions of ${\bf
R}_{W}$ and ${\bf R}_{T}$, we have again, ${\bf C}_{B}={\bf C}_{T}-{\bf
C}_{W}$, and can therefore write the inter-subject correlation in the non-
linear space as a function of the new parameters $\bm{\alpha}$ as follows:

$\rho=\frac{{\bf v}^{\top}{\bf R}_{B}{\bf v}}{{\bf v}^{\top}{\bf R}_{W}{\bf
v}}=\frac{\bm{\alpha}^{\top}{\bf C}_{B}\bm{\alpha}}{\bm{\alpha}^{\top}{\bf
C}_{W}\bm{\alpha}}.$ (76)

This can be solved again as an eigenvalue problem but now with the within- and
between-subject covariance matrices of the kernel vectors $\bar{\bf
k}_{i}^{*l}$ and the corresponding projection vector $\bm{\alpha}$. Note that
a similar result is obtained for the kernel version of CCA, namely, the
conventional CCA algorithm is applied to kernel matrices instead of the
original data [67].

To demonstrate the method, we generated a simple 2-dimensional example where
the shared dimension is non-linear, and this non-linear relationship has to be
discovered from the data alone. As with the example for linear CorrCA (Figure
1), we have simulated $N=2$ “subjects”, each providing two signals ($D=2$),
now with $T=40$ samples in time (Figure 11 A & B). In order to test whether
the algorithm finds the correct non-linear mapping we generate data with a
known non-linear relationship. In this case the two subjects share the
amplitude (distance from the origin in the 2-dimensional input space of ${\bf
x}$; see Figure 11 C). However, the phase angle in this 2-dimensional plane is
selected randomly at each time point and for each subject. We apply kernel-
CorrCA using Gaussian kernels [63] keeping the time axis unchanged
($T=T^{\prime}$; see also our companying code). The first component dimension
found by kernel-CorrCA, $y_{1}$, is approximately linear with the amplitude
and is independent of phase in the original 2D plane (Figure 11 G & H). The
algorithm therefore discovered in the first component amplitude as the shared
non-linear dimension. The second component dimension, $y_{2}$ is essentially a
nonlinear function of component $y_{1}$ (Figure 11 F), hence it also has high
ISC. The specific non-linearity between the two components is arbitrary, but
does enforce that the two variables are uncorrelated. We obtain nearly
identical results regardless of whether we use Gaussian or Tanh kernels [see
code and 63], although the specific relationship between $y_{1}$ and $y_{2}$
differs.

Figure 11: Kernel-CorrCA. In this demonstrative 2D example, two subjects share
amplitude but have a random phase. The goal is for the algorithm to discover
the common amplitude signal — a non-linear transformation of the data — from
60 samples. (A-B) Two dimensions of the original data with one trace for each
subject. (C) Samples are connected between the two subjects with gray lines,
as in Figure 1. Circles of varying radii are drawn for visualization purposes.
Note that connected points have the same distance from the centers, i.e.
subjects share an amplitude. (D-E) After kernel-CorrCA using Gaussian kernels,
signals between subjects are strongly correlated (ISC=1). (F) Points are in a
new non-linear space, which is necessarily curved to decorrelate the two
dimensions. (G-H) Colored lines correspond to the circles in panel C. The
first dimension $y_{1}$ is independent of phase and grows linearly with
amplitude. This means that the algorithm identified the non-linear dimension
that is shared between the two subjects.

##  Appendix G Further remarks on Inter-Subject Correlation

###  G.1 Relationship of ISC to Intra-Class Correlation

The definition of inter-subject correlation in the present work is similar to
the definition of intra-class correlation (ICC) first introduced by [68].
Fisher’s ICC is referred to as pairwise inter-class correlation and differs
from Pearson’s correlation coefficient [69]. A more modern definition of
intra-class correlation, which in the present context would allow for a
different number of subjects per sample (i.e. $N_{i}$), was introduced by
[70]:

$r_{p}=\frac{\sum_{i=1}^{T}\frac{1}{N_{i}-1}\sum_{k=1}^{N_{i}}\sum_{l=1,k\neq
l}^{N_{i}}(y_{i}^{l}-{\bar{y}}_{*}^{*})(y_{i}^{k}-{\bar{y}}_{*}^{*})}{\sum_{i=1}^{T}\sum_{l=1}^{N_{i}}(y_{i}^{l}-{\bar{y}}_{*}^{*})^{2}}.$
(77)

The definition of pairwise ICC, $r_{p}$, differs from our definition of ISC,
$\rho$, in that it subtracts the total mean $\bar{y}_{*}^{*}$, whereas we
subtract the sample mean of each subject, $\bar{y}_{*}^{l}$. The two are
identical in the case of unbiased ratings, i.e. ${\bf S}_{M}=0$ defined in Eq.
(21). A definition comparable to ISC can be found in [6] in the context of
MCCA under the term “modified inter-class correlation”, but is not commonly
seen in the literature. Instead, the generalized eigenvalue solution to MCCA
presented here has been often derived using the sum of correlations, subject
to a normalization constraint [e.g. 28, 71, 29].

In Table 2 we summarize various uses of these correlation measures. For each
case the indices play different roles. In the context of ICC, index $i$
represents different classes whereas $l,k$ are different exemplars of a class.
Each class can have a different number of exemplars $N_{i}$. In the context of
ISC, $l,k$ are subjects and $i$ represent measures taken, perhaps in time. In
the case of ratings, $l,k$ are the raters and $i$ are the items being rated.
In Karlin et al. [70], the measure $r_{p}$ is used to assess intra-familial
correlations (correlations between siblings). In this case, $k,l$ indexes the
siblings, whereas $i$ enumerates the families. The algorithms listed in Table
2 aim to extract components in multi-dimensional data that maximize the
corresponding statistic. Where no algorithm is listed, either CorrCA or LDA
could be used.

Table 2: Meaning of different indices used in the definitions of various measures of correlation.  Algorithm |  Statistic |  i=1 …T |  l,k=1 …N |  Reference
---|---|---|---|---
CorrCA | inter-subject correlation | samples | subjects | present paper
JD | Signal-to noise ratio | samples | repeats | [23]
LDA | between class separation | classes | exemplars | [18]
| intra-class correlation | classes | exemplars | [69]
| intra-rater correlation | items | raters | [2]
| intra-familial correlation | families | siblings | [70]
MCCA | modified intra-class correlation | samples | datasets | [6]
MCCA | inter-subject correlation | samples | subjects | present paper

The definition of ISC/IRC in Eq. (2) can be thought of as a Pearson
correlation coefficient, where we concatenate the signals (after subtracting
the individual subject/rater means) corresponding to all possible pairs of
subjects. In contrast to conventional Pearson correlation, with this symmetric
definition (where we include pairs $lk$ as well as $kl$), the signals of
different subjects have to have identical scale in order to achieve perfect
ISC/IRC (see section G.2). In the case of ratings, this means that different
raters can only differ in their mean ratings but have to otherwise provide
identical deviations from the mean in order to achieve perfect inter-rater
correlation. This means that this metric is sensitive to multiplicative bias
(i.e. each rater/subject has a different scale). The issue of bias in the
context of inter-rater agreement is discussed in more detail by [2], who argue
that Pearson correlation should be used to assess test-retest reliability to
account for learning effects of subjects from one test to the next. If there
is an undesirable proportional bias, then one can easily correct for this bias
by standardizing the data, i.e. dividing by the standard deviation for each
subject/rater in each dimension prior to applying CorrCA.

As we have shown here, maximizing IRC is equivalent to maximizing inter-rater-
reliability, where we have defined reliability in Eq. (25) as the variance of
the mean over the mean of the variance. For instance, if the rating is the
speed of performing a task (as in the application 3.2), this would mean that a
100 millisecond variation across repeats is fairly reliable for a task that
takes 5 seconds to complete, but really quite unreliable if the task only
takes 0.5 second in average.

It should be noted that, when ratings are on a discrete numerical scale, IRC
may not be an appropriate measure of reliability. For instance, a rating on an
integer scale of 1, 2, 3, 4, 5, may have an inter-rater variability of $\pm
1.0$. With the present definition of IRC — $\rho$ in Eq. (2) — this same
variability would be considered half as reliable if the mean rating is 4 as
compared to a mean of 2. Whether this is desirable or not may depend on the
specific ratings statistic. Finally, IRC can only be used meaningfully for
ratings that are on a proportional scale. For categorical ratings one should
use instead the conventional Cohen’s $\kappa$ or similar inter-reliability
measures.

###  G.2 ISC is normalized

To our knowledge, the present definition of $\rho$ in Eq. (2) does not appear
in this exact form in the literature, except for our previous work [14].
Therefore, we take some time here to show that this definition of correlation
has indeed a maximal value of 1, and that this occurs only when the all
signals are identical across subjects (not just proportional as in Pearson’s
correlations), except for a subject-dependent mean, which is permitted.

Without loss of generality, let us assume that all signals are zero mean,
${\bar{y}}^{l}_{*}=0$. To show that $\rho\leq 1$, we can also demonstrate,
equivalently, that the following expression is non-negative:

$\displaystyle 0$ $\displaystyle\leq$ $\displaystyle(N-1)r_{W}-r_{B}$ (78)
$\displaystyle=$
$\displaystyle(N-1)\sum_{l}\sum_{i}(y^{l}_{i})^{2}-\sum_{l}\sum_{k\neq
l}\sum_{i}y^{l}_{i}y^{k}_{i}$ $\displaystyle=$ $\displaystyle
N\sum_{l}\sum_{i}(y^{l}_{i})^{2}-\sum_{i}\left(\sum_{l}y^{l}_{i}\right)^{2}=E.$

The last expression is abbreviated as $E$. To find the $y^{l}_{i}$ that
minimizes $E$ we find the location where the gradient is zero and check
whether the curvature at that location is non-negative.

$\displaystyle\frac{\partial E}{\partial y_{j}^{k}}$ $\displaystyle=$
$\displaystyle 2Ny^{k}_{j}-2\sum_{i}\sum_{l}y^{l}_{i}\sum_{l}\delta^{lk}_{ij}$
(79) $\displaystyle=$ $\displaystyle 2Ny^{k}_{j}-2\sum_{l}y^{l}_{i},$
$\displaystyle\frac{\partial^{2}E}{\partial y_{j}^{k}\partial y_{j}^{l}}$
$\displaystyle=$ $\displaystyle 2N\delta^{kl}_{ji}-2\sum_{n}\delta^{kn}_{ij}$
(80) $\displaystyle=$ $\displaystyle 2\delta_{ij}(N\delta^{kl}-1).$

Solving for $\partial E/\partial y=0$ yields an unique solution:

$y_{j}^{k}=\frac{1}{N}\sum_{l}y_{j}^{l}=\bar{y}^{*}_{j}.$ (81)

This solution is a linear manifold, as one can pick arbitrary scale for
$\bar{y}^{*}_{j}$. Inserting this into (78) yields $E=0$ at that point. To
determine whether this is indeed a minimum value for $E$, and hence $\rho\leq
1$, we need to show that the Jacobian (80) has only non-negative eigenvalues.
Because of the $\delta_{ij}$ in (80), we can find the eigenvalues for each
$ij$ separately, i.e. the eigenvalues of $J=N{\bf I}-{\bf 1}$, where ${\bf I}$
is the $N$-dimensional identity matrix and ${\bf 1}$ is a $N$-dimensional
vector with all values set to 1. This ${\bf J}$ is the same for all $ij$, and
the eigenvalue equation is for ${\bf y}\in\mathbb{R}^{N}$:

${\bf J}{\bf y}=\lambda{\bf y}.$ (82)

The characteristic polynomial for these eigenvalues reads: