#  The MCC approaches the geometric mean of precision and recall as true
negatives approach infinity.

Jon Crall
Kitware Inc.
<EMAIL_ADDRESS>https://github.com/Erotemic/<EMAIL_ADDRESS>
###### Abstract

The performance of a binary classifier is described by a confusion matrix with
four entries: the number of true positives (TP), true negatives (TN), false
positives (FP), and false negatives (FN).

The Matthew’s Correlation Coefficient (MCC), F1, and Fowlkes–Mallows (FM)
scores are scalars that summarize a confusion matrix. Both the F1 and FM
scores are based on only three of the four entries in the confusion matrix
(they ignore TN). In contrast, the MCC takes into account all four entries of
the confusion matrix and thus can be seen as providing a more representative
picture.

However, in object detection problems, measuring the number of true negatives
is so large it is often intractable. Thus we ask, what happens to the MCC as
the number of true negatives approaches infinity? This paper provides insight
into the relationship between the MCC and FM score by proving that the FM-
measure is equal to the limit of the MCC as the number of true negatives
approaches infinity.

_Keywords_ Confusion Matrix $\cdot$ Binary Classification $\cdot$
Fowlkes–Mallows Index $\cdot$ Matthew’s Correlation Coefficient $\cdot$ F1

##  1 Introduction

Evaluation of binary classifiers is central to the quantitative analysis of
machine learning models [1]. Given a finite set of examples with known real
labels, the quality of a set of corresponding predicted labels can quantified
using a $2\times 2$ confusion matrix. A confusion matrix counts the number of
true positives (TP), true negatives (TN), false positives (FP), and false
negatives (FN) a model predicts with respect to the real labels. The confusion
matrix is often written as:

$\begin{bmatrix}{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}&{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\\\
{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}&{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\end{bmatrix}$
(1)

This matrix provides a holistic view of classifier quality, however, it is
often desirable to summarize performance using fewer numbers. Two popular
metrics defined on a classification matrix are precision and recall.

Precision — also known as the positive-predictive-value (PPV) — is the
fraction of positive predictions that are correct.

$\texttt{{PPV}}{}=\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}}$
(2)

Recall — also known as the true positive rate (TPR), sensitivity, or
probability of detection (PD) — is the fraction of real positive cases that
are correct.

$\texttt{{TPR}}{}=\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}$
(3)

One of the most popular confusion metrics is the F1 score. It can be defined
as the harmonic mean of precision and recall.

$\texttt{{F1}}{}=\frac{2\texttt{{PPV}}{}\cdot\texttt{{TPR}}{}}{\texttt{{PPV}}{}+\texttt{{TPR}}{}}=\frac{2{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{2{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}$
(4)

A similar, but lesser used metric is the Fowlkes–Mallows index [2], which was
originally developed for measuring the similarity between two clusterings of a
set of points. It can be defined as the geometric mean of precision and recall
[3].

$\texttt{{FM}}{}=\sqrt{\texttt{{PPV}}{}\cdot\texttt{{TPR}}{}}=\sqrt{\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}}\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}}$
(5)

In [1], Powers notes that the F1 score (and consequentially any metric that
only includes precision and recall) only takes into account three of the four
measures in a confusion matrix. Powers, introduces modifications of precision
and recall he refers to as informedness and markedness. Additionally he
advocates for the use of the MCC over the F1 measure.

The Matthews Correlation Coefficient (MCC) [4] accounts for all four terms in
the confusion matrix and is defined as:

$\texttt{{MCC}}{}=\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}\cdot{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}-{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\cdot{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})({\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}}$
(6)

While the MCC is a desirable measure due to its balanced inclusion of all
terms in the confusion matrix, it requires that the number of true negatives
is measurable. In the case of object detection problems [5], this is often
intractable as the number of the number of predicted boxes and missed true
boxes is dwarfed by the total number of boxes that the system correctly did
not predict. One can see this by considering the set of all $N\times M$ boxes
centered at each pixel, most of which will be considered true negatives. If
the width and height of the boxes are allowed to extend outside the image,
then the number of predictable boxes actually is unbounded (and even if they
must be contained in the image, there will still be a very large number of
them in real world cases).

Because calculating the number of true negatives is difficult for open-world
problems like object detection, it is conceptually simpler to ignore true
negatives and simply focus on the much smaller set of true positives, false
positives, and false negatives, which can be used to compute PPV, TPR, F1, and
FM. While these measures have proven themselves to be effective, simply
ignoring true negatives is somewhat unsatisfying. We seek to remedy this
noting that in these open-world problems the number of true negatives is so
large it is effectively infinite and thus we ask the question: what happens to
the MCC as the number of true negatives approaches infinity?

The main contribution of this paper is to highlight a relationship between the
MCC and the FM score. The MCC reduces to FM as the number of true negatives
approaches infinity. While this is not a difficult result to show, to the best
of the author’s knowledge, this was first shown in a blog post [6], but has
not yet been published. This paper is a more formal description of this
result.

##  2 The Relationship Between MCC and FM

### Taking the limit of the MCC

Consider the limit of the MCC as the number of true negatives approaches
infinity.

$\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\texttt{{MCC}}{}=\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}\cdot{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}-{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\cdot{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})({\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}}$
(7)

We can take this limit by applying some algebra to the body of the limit. We
multiply the numerator and denominator by
$\frac{1}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}}$:

$=\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\frac{\frac{1}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}}({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}\cdot{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}-{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\cdot{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}{\frac{1}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}}\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})({\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}}$
(8)

We distribute the $\frac{1}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}}$
term in the numerator and denominator:

$=\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\frac{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}-{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\cdot\frac{{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}})}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})(\frac{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}})(\frac{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}})}}$
(9)

The
$\frac{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}}$
terms in the denominator cancel:

$=\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\frac{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}-{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\cdot\frac{{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}})}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})(1+\frac{{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}})(1+\frac{{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}})}}$
(10)

The terms involving TN are fractions of simple rational polynomials (w.r.t.
TN) and in each case the degree of the denominator is greater than that of the
numerator, so in the limit each of these terms simplifies to $0$. Thus, the
entire equation simplifies to:

$=\frac{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}-{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}\cdot
0)}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})(1+0)(1+0)}}$
(11)

Thus we find that the limit of the MCC as true true negatives approach
infinity is:

$=\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}}$
(12)

### Rearranging the FM

Now rearranging the equation for FM, we find it is equivalent to the limit of
the MCC as the number of true negatives approaches infinity.

FM
$\displaystyle=\sqrt{\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}}}\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}}}}$
(13)
$\displaystyle=\sqrt{\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}^{2}}{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}}$
(14)
$\displaystyle=\frac{{\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}}{\sqrt{({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.95703125,0.15625,0.2109375}\texttt{{FP}}})({\color[rgb]{0.2421875,0.68359375,0.16796875}\texttt{{TP}}}+{\color[rgb]{0.9375,0.46484375,0.140625}\texttt{{FN}}})}}$
(15)
$\displaystyle=\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\texttt{{MCC}}{}$
(16)

### Verifying the proof

The correctness of these claims can be verified using SymPy [7]. We define a
symbolic expression for the definition of the MCC and FM score. We then use
SymPy to determine the limit of the MCC as
${\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty$. Finally we subtract
expressions that we claim are equal, which will result in zero only if they
are equal.

⬇

1

2from sympy import sqrt, symbols, simplify

3from sympy.series import limit

4

5tp, tn, fp, fn = symbols("tp tn fp fn",

6 integer=True, negative=False)

7

8# The definition of the MCC

9numer = (tp * tn - fp * fn)

10denom = sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn))

11mcc = numer / denom

12

13# The definition of FM

14FM = sqrt((tp / (tp + fn)) * (tp / (tp + fp)))

15

16# Compute the limit of the MCC definition

17mcc_lim = limit(mcc, tn, float(’inf’))

18

19# We claim the limit of the MCC and the FM are equivalant to:

20mcc_lim_claim = tp / sqrt((tp + fn) * ((tp + fp)))

21

22# Check the claim is equal to FM

23assert simplify(FM - mcc_lim_claim) == 0

24# Check the claim is equal to the MCC limit

25assert simplify(mcc_lim - mcc_lim_claim) == 0

The above program does not raise an AssertionError, thus we have proven
$\lim_{{\color[rgb]{0,0.40625,0.78125}\texttt{{TN}}}\to\infty}\texttt{{MCC}}{}=\texttt{{FM}}{},\square{}$.

##  3 Conclusion

This paper proves that the limit of the MCC as the number of true negatives
goes to infinity is equivalent to the Fowlkes–Mallows index (i.e. the
geometric mean of precision and recall).

This is a useful insight in open world problems where the number of true
negative cases grows faster than the number of other confusion categories. It
validates the use of precision and recall as a way of describing the quality
of object detection results and hints that the FM score may be a preferable
alternative to the more standard F1 score.

## References

  * Powers [2011] David Martin Powers.  Evaluation: from precision, recall and F-measure to ROC, informedness, markedness and correlation.  volume 2, pages 37–63, February 2011.  URL http://dspace2.flinders.edu.au/xmlui/handle/2328/27165.  00000\.
  * Fowlkes and Mallows [1983] E. B. Fowlkes and C. L. Mallows.  A Method for Comparing Two Hierarchical Clusterings.  volume 78, pages 553–569, 1983.  doi:10.2307/2288117.  URL https://www.jstor.org/stable/2288117.  Publisher: [American Statistical Association, Taylor & Francis, Ltd.].
  * Tharwat [2020] Alaa Tharwat.  Classification assessment methods.  _Applied Computing and Informatics_ , 17(1):168–192, January 2020.  ISSN 2210-8327.  doi:10.1016/j.aci.2018.08.003.  URL https://doi.org/10.1016/j.aci.2018.08.003.  Publisher: Emerald Publishing Limited.
  * Matthews [1975] B. W. Matthews.  Comparison of the predicted and observed secondary structure of T4 phage lysozyme.  volume 405, pages 442–451, October 1975.  doi:10.1016/0005-2795(75)90109-9.  URL https://www.sciencedirect.com/science/article/pii/0005279575901099.
  * Zou et al. [2023] Zhengxia Zou, Keyan Chen, Zhenwei Shi, Yuhong Guo, and Jieping Ye.  Object detection in 20 years: A survey.  _Proceedings of the IEEE_ , 2023.
  * Crall [2019] Jon Crall.  Closed form of the mcc when tn->inf, 2019.  URL https://erotemic.wordpress.com/2019/10/23/closed-form-of-the-mcc-when-tn-inf/.  [Online; accessed 30-April-2023].
  * Meurer et al. [2017] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz.  Sympy: symbolic computing in python.  volume 3, page e103, January 2017.  doi:10.7717/peerj-cs.103.  URL https://doi.org/10.7717/peerj-cs.103.


# Message-Relevant Dimension Reduction of Neural Populations

Amanda Merkley and Pulkit Grover  Department of Electrical & Computer
Engineering
Carnegie Mellon University, Pittsburgh, USA
Email: {amerkley<EMAIL_ADDRESS>Alice Y. Nam and Y. Kate Hong
Department of Biological Sciences
Carnegie Mellon University, Pittsburgh, USA
Email: {aynam<EMAIL_ADDRESS>
###### Abstract

Quantifying relevant interactions between neural populations is a prominent
question in the analysis of high-dimensional neural recordings. However,
existing dimension reduction methods often discuss communication in the
absence of a formal framework, while frameworks proposed to address this gap
are impractical in data analysis. This work bridges the formal framework of
$M$-Information Flow with practical analysis of real neural data. To this end,
we propose Iterative Regression, a message-dependent linear dimension
reduction technique that iteratively finds an orthonormal basis such that each
basis vector maximizes correlation between the projected data and the message.
We then define ‘$M$-forwarding’ to formally capture the notion of a message
being forwarded from one neural population to another. We apply our
methodology to recordings we collected from two neural populations in a
simplified model of whisker-based sensory detection in mice, and show that the
low-dimensional $M$-forwarding structure we infer supports biological evidence
of a similar structure between the two original, high-dimensional populations.

##  I Introduction

Noisy, high-dimensional data are collected in abundance in neuroscience
experiments. It is now common to record tens to hundreds of neurons from
multiple neural populations [1, 2], where each population is a functionally
distinct group of neurons in the brain. These advances are enabling
neuroscientists to ask questions about communication and information flow in
real neural data to gain insight into population-level interactions. There is
also a parallel interest to address these questions using information-
theoretic formalisms [3, 4, 5]. A specific direction we took in earlier works
[6, 7, 8] was “$M$-Information Flow”, a formal framework in which we defined
information flow of $M$, a message of interest, and provided a technique to
track the flow of $M$ through a network in time.

However, the rigorous definitions we laid out in [8] are impractical in the
analysis of noisy, high-dimensional neural datasets. To address this issue, we
consider the problem of message-relevant dimension reduction from first
principles. We seek an estimable representation of neural activity that
summarizes information about the message, $M$, in neural data based on a
measure of $M$-relevance. We propose Iterative Regression (IR) (Sec. II-B) as
a linear dimension reduction technique that uses sample correlation between
$M$ and the projected data as the measure of $M$-relevance. IR iteratively
finds an orthonormal basis to represent the neural data by computing the
linear projection that maximizes correlation with the message.

Our second contribution is inference of a low-dimensional communication
structure by conditional dependence testing using data reduced with IR.
Inspired by $M$-Information Flow, we define ‘$M$-forwarding’ in Sec. III as a
conditional independence relationship to formally capture the notion of how a
message $M$ is forwarded from one neural population to another. While
$M$-forwarding generally describes a relationship between neural populations
(i.e., vectors of activity), we estimate this quantity with one-dimensional
(1D) representations of neural populations (i.e., scalars). Inferring an
$M$-forwarding structure using IR-reduced data not only overcomes the
challenge of directly applying the definition to high-dimensional datasets,
but can also reveal the existence (or absence) of a low-dimensional
representation of $M$-forwarding.

To assess the performance of our approach on real data, we analyze a
simplified network of the whisker-based sensory detection system in mice. We
first apply IR to neural population data we acquired and show that IR captures
strong correlation with $M$ in each population (Sec. IV-B). We compare the
correlations obtained from IR to those from demixed PCA (dPCA) [9] and model-
based Targeted Dimension Reduction (mTDR) [10], two techniques developed for
message-dependent dimension reduction in neuroscience. We also compare to a
denoising sparse autoencoder (DSAE). The DSAE of [11] is a shallow
autoencoder, suitable for our dataset size, used in neural signal
classification. We adapt this model to see if nonlinear dimension reduction
provides an advantage in spite of not explicitly accounting for $M$. We find
that IR maintains clear $M$-relevance, while it is diminished in other methods
(Sec. IV-B).

Next, we use conditional dependence testing to determine if 1D representations
of both neural populations do _not_ harbor an $M$-forwarding structure (Sec.
IV-C). These tests allow us to infer a direct pathway between $M$ and a
dimension-reduced neural population, and that the message is not $M$-forwarded
from the other reduced population. Using IR to obtain the 1D representation of
each area allows us to conclude dependence of $M$ with both of the low-
dimensional neural populations as well as a direct pathway between $M$ and one
of the IR-reduced populations, suggesting a low-dimensional structure
supported by previous anatomical evidence [12]. Furthermore, most other
methods do not even detect unconditional dependence with one population. Only
dPCA detects the same communication structure as IR, and to less confidence
than IR.

##  II Linear Dimension Reduction

###  II-A A measure of $M$-relevance

We first address the problem of obtaining dimensions relevant to the message,
$M$, from an $n$-dimensional dataset111In the rest of this work, we assume $M$
and $X$ are mean-centered without loss of generality., $X$. Naturally, the
choice of dimension reduction depends on how relevance is measured. In this
work, we measure $M$-relevance as sample correlation of the reduced-
dimensional data with $M$ since correlation is both estimable and
interpretable. As such, we seek projections of the data that retain the most
correlation with $M$. The first projection is

$\displaystyle v_{1}=\operatorname*{argmax}_{u\in\mathbb{R}^{n}}r(u^{T}X,M)$
(1)

where $r(u^{T}X,M)$ is the sample correlation between $u^{T}X$ and $M$. In
fact, we show in Appendix A that (1) is equivalent to finding the linear
minimum mean square error (LMMSE) estimate of $M$, to which the solution is
ordinary least squares regression, i.e.,
$v_{1}=(XX^{T})^{-1}XM=\Sigma_{X}^{-1}\Sigma_{XM}$. When $X$ and $M$ are
jointly Gaussian, the regression vector is also a sufficient transformation of
$X$, which we show in Theorem 2 in Appendix A.

However, the regression vector may not be sufficient for representing neural
data in general. As such, $M$ may not be recoverable by a single linear
projection of $X$. A key observation is the existence of nonzero correlations
in other dimensions. In other words, data projected onto $u\neq v_{1}$ may
still be significantly correlated to $M$. Indeed, for $\mathbb{E}[u^{T}XM]=0$
to hold, the projection vector $u$ must be orthogonal to $\Sigma_{XM}$. This
motivates IR as a method for recovering an orthonormal basis such that each
basis vector iteratively maximizes $M$-relevance in decreasing dimensions
after removing contributions from previous regression vectors. We consider
linear projections to align with the dominant form of dimension reduction on
neural data [13]. The focus on linear dimension reduction in neuroscience is
partly due to the statistical versatility of linear techniques. Each dimension
is also interpretable as a ‘conceptual neuron’ that linearly combines real
neurons to summarize some aspect of the data.

###  II-B Iterative Regression

We now describe IR in detail as a technique for extracting $M$-relevant
projections one dimension at a time. Suppose $X\in\mathbb{R}^{n\times Q\times
T_{0}}$ is the dataset of $n$ neurons recorded for $T_{0}$ time points over
$Q$ trials, where a trial is a time series of neural activity in response to a
randomized stimulus. Let $M\in\mathbb{R}^{Q}$ be the message sampled over $Q$
trials. Define $X_{n}\in\mathbb{R}^{n\times Q}$ to be a slice of the dataset
at one time point. The first IR dimension, $v_{1}^{*}$, is the normalized
regression vector from (1), i.e., $v_{1}^{*}=v_{1}/\|v_{1}\|_{2}$. We assume
$n<Q$, an assumption usually satisfied in experiments and is otherwise
resolved by removing unresponsive neurons.

IR builds an orthonormal basis from the frame of reference of the first
regression vector by finding correlations in a reduced $q$-dimensional space
after removing contributions from the previous $(n-q)$ IR dimensions. To
illustrate this process for the second IR dimension, consider the dataset
after removing the contributions from $v_{1}^{*}$:

$\displaystyle X_{n-1}=X_{n}-v_{1}^{*}(v_{1}^{*})^{T}X_{n}.$

Whereas $X_{n}$ was the dataset of real neural activity, $X_{n-1}$ is now
interpreted as the activity of conceptual neurons. By removing a column from
$X_{n-1}$ (the redundant conceptual neuron), the regression vector for the
lower-dimensional space is again computed for the reduced set of conceptual
neurons to find the next IR dimension $v_{2}\in\mathbb{R}^{n-1}$. Then,
$v_{2}$ is normalized, projected back into the $n$-dimensional space, and
orthogonalized with respect to $v_{1}^{*}$. This process, described in
Algorithm 1, is repeated for all $n$ IR dimensions at each time.

Algorithm 1 Iterative Regression

1:$X\in\mathbb{R}^{n\times Q\times T_{0}}$, $M\in\mathbb{R}^{Q}$

2:$M\leftarrow M-\bar{M}$ $\triangleright$ Mean center

3:for $t=1,...,T_{0}$ do

4: for $i=1,...,n$ do

5: $X_{0}\leftarrow X(t)-\bar{X}(t)$ $\triangleright$ Mean center wrt $n$-axis

6: for $1\leq j<i$ do

7: $s_{j}\leftarrow v_{j}(t)^{T}X_{0}$

8: $s_{j}\leftarrow s_{j}/\|s_{j}\|_{2}$

9: $X_{0}\leftarrow X_{0}-X_{0}s_{j}^{T}s_{j}$

10: end for

11: $v_{q}\leftarrow(X_{0}X_{0}^{T})^{-1}X_{0}M$

12: $v_{GS}\leftarrow$ Gram-Schmidt($v_{q},\\{\varnothing,...,v_{i-1}(t)\\}$)

13: $v_{i}(t)\leftarrow v_{GS}\;/\;\|v_{GS}\|_{2}$

14: end for

15:end for

Since neural spiking data is often collected as a $n\times Q\times T_{0}$
tensor, the dataset must be converted into a matrix so it is amenable to
dimension reduction with IR. A common technique for reshaping the tensor is to
average over trials so the new dataset $X_{average}\in\mathbb{R}^{n\times
T_{0}}$ indexes neurons and time [14]. Dimension reduction is then performed
on $X_{average}$. This operation effectively treats each time point as an
independent sample of neural activity, when, by design of the experiment,
activity within each trial is not even stationary.

Instead, we apply dimension reduction per time point. This is to better
respect the nonstationary nature of the trial, acknowledging that there may be
dependence between time points. Furthermore, IR is a correlation-maximizing
method applied to neurons that may be anti-correlated with $M$. Applying IR on
the averaged dataset $X_{average}$ returns dimensions that maximize
correlation over the whole trial, which may actually result in dimensions that
are the most uncorrelated representation of $M$. By performing dimension
reduction per time point, the effect of anti-correlated neurons is accounted
for simply by flipping signs at each time.

###  II-C Related methods: Linear dimension reduction

IR is related to several classes of linear dimension reduction techniques,
including principal component analysis (PCA), canonical correlation analysis
(CCA), and partial least squares (PLS) regression. To see this, consider the
IR, PCA, CCA, and PLS objectives for finding the first dimension222We assume
that $M\in\mathbb{R}$ and $X\in\mathbb{R}^{n}$.
$\Sigma_{X}\in\mathbb{R}^{n\times n}$ is the covariance of $X$ and
$\Sigma_{MX}\in\mathbb{R}^{1\times n}$ is the cross covariance between $M$ and
$X$.:

  * •

IR objective

$\displaystyle\max_{u}\frac{\Sigma_{MX}u}{\sqrt{\sigma_{m}^{2}}\sqrt{u^{T}\Sigma_{X}u}}$
(2)

  * •

PCA objective [15]

$\displaystyle\max_{u}\frac{u^{T}\Sigma_{X}u}{u^{T}u}$ (3)

  * •

CCA objective [16]

$\displaystyle\max_{v,u}\frac{v^{T}\Sigma_{MX}u}{\sqrt{v^{T}\sigma_{m}^{2}v}\sqrt{u^{T}\Sigma_{X}u}}$
(4)

  * •

PLS objective [17]

$\displaystyle\max_{v,u}\frac{v^{T}\Sigma_{MX}u}{\sqrt{v^{T}v}\sqrt{u^{T}u}}$
(5)

We group the techniques by how they treat the message. IR and PLS are both
message-dependent dimension reduction techniques, while PCA and CCA are not
message-dependent. Although (4) suggests that CCA is also message-dependent
through $\Sigma_{MX}$, CCA is ill-suited for dimension reduction of $X$ when
$M$ is 1D. Standard CCA was developed for dimension reduction of two
homogeneous objects (e.g., both variables represent neural activity) unlike
PLS, which was developed for robust regression [18] (e.g., one variable is a
stimulus and another is neural activity). The algorithmic description of CCA
and PLS manifests this difference: the max number of dimensions standard CCA
can find is limited to the minimum of dimensions of the two variables [19]
(i.e. one dimension for a 1D message) while PLS is a sequential algorithm that
squeezes out covariance after successive ‘deflations’ of the data matrix $X$
(i.e. up to $n$ dimensions).

Although PCA is not message-dependent, $M$-relevance can still be assigned to
each Principal Component (PC). In fact, PCA is related to IR by the projection
of each PC onto the first IR dimension. For a given PC, the magnitude of this
projection, scaled by the PC’s explained variance, is proportional to its
$M$-relevance. Theorem 1 summarizes this relationship.

######  Theorem 1.

Let $X\in\mathbb{R}^{n\times Q}$ be the mean-centered dataset and
$M\in\mathbb{R}^{Q}$, the mean-centered message vector. Suppose
$v=(XX^{T})^{-1}XM$ is the regression vector and $\\{p_{i}\\}_{i=1}^{n}$ is
the set of principal components of $\Sigma_{X}\propto XX^{T}$. Then,

$\displaystyle r(p_{i}^{T}X,M)=\alpha\sqrt{\lambda_{i}}\ \langle
v^{*},p_{i}\rangle,$ (6)

where $r(p_{i}^{T}X,M)$ is the sample correlation coefficient between
$p_{i}^{T}X$ and $M$, $\lambda_{i}$ is the eigenvalue associated to $p_{i}$,
$v^{*}$ is the first IR dimension, and $\alpha=\|v\|/s_{M}$ is a constant with
$s_{M}$ being the standard deviation of $M$.

###### Proof.

See Appendix A. ∎

PLS optimizes symmetrically over $M$ and $X$, where the additional
optimization parameter $v$ of PLS yields a low-dimensional representation of
$M$. However, since $M$ is a design variable of the experiment, we are only
interested in the projection $u$ of the data $X$. In summary, IR can be viewed
as: (1) a message-dependent formulation of PCA, (2) an asymmetric, sequential
formulation of CCA, or (3) an asymmetric, correlation-based formulation of
PLS.

###  II-D Related methods: Dimension reduction in neuroscience

dPCA [9] and mTDR [20] are respective adaptations of PCA and regression for
message-dependent dimension reduction in the context of neuroscience. dPCA
‘demixes’ the contributions of different messages (e.g., experimental stimuli,
animal behavioral features, and time indices) in the low-dimensional
representation of neural activity [9]. dPCA achieves this separation by first
decomposing the original neural data into components marginalized with respect
to a specific message (or grouping of messages), then performs a low-rank
matrix factorization on each component with the goal of minimizing overall
reconstruction error [9], akin to regular PCA on the marginalized covariance.
This differs from IR, which performs dimension reduction with $M$ directly,
rather than a representation of $M$ derived from marginalized activity.

mTDR is a regression-based method [20]. mTDR assumes a high-dimensional linear
regression model of neural activity in response to multiple messages for each
neuron and each time point, then performs a low-rank matrix factorization of
the model. Unlike IR, mTDR finds all dimensions simultaneously. Although this
may be advantageous for evaluating the entire low-dimensional model, the
$M$-relevance of specific dimensions may be spread across multiple dimensions.
Because IR is necessarily sequential, each dimension can be optimized for
$M$-relevance. As shown in Sec. IV-C, this may be crucial for inferring an
accurate summarization of communication.

dPCA and mTDR are the primary message-dependent techniques in a large
literature on dimension reduction in neuroscience [13]. While other techniques
do not explicitly account for a message, they still discuss communication as
low-dimensional interactions between neural populations by using various
correlational measures [21, 22, 23]. Some of these works [23] can be viewed as
dimension-reduced extensions of Granger causality [24] and directed
information, which are also correlational measures of interaction, and have
been of interest within information-theoretic literature on neuroscience [25].
Nonlinear approaches in the form of autoencoders have also been explored,
ranging from shallow DSAEs [11] to deep, sequential networks [26]. Although
not message-dependent, these approaches may offer advantages in more flexible,
nonlinear modeling.

##  III $M$-Forwarding

We now propose $M$-forwarding as a formal definition to embody the concept of
how a message $M$ is forwarded between two disjoint populations
$X_{1}\in\mathbb{R}^{n_{1}}$ and $X_{2}\in\mathbb{R}^{n_{2}}$, for $n_{1}$
neurons in $X_{1}$ and $n_{2}$ neurons in $X_{2}$. Three scenarios are
possible: $X_{1}$ forwards $M$ to $X_{2}$, $X_{2}$ forwards $M$ to $X_{1}$, or
neither (i.e., no forwarding occurs or both forward to each other). Inspired
by the measure for $M$-Information Flow on an edge as positive conditional
mutual information [8], we distinguish between the first two scenarios by
formalizing what it means for one area to forward a message to another:

######  Definition 1 ($M$-forwarding333Definition 1 is the same as that of a
physically degraded broadcast channel.).

For $M\in\mathbb{R}$, $C_{1}\in\mathbb{R}^{n_{1}}$, and
$C_{2}\in\mathbb{R}^{n_{2}}$, where $n_{1},n_{2}\in\mathbb{N}$, we say $C_{1}$
$M$-forwards to $C_{2}$ if the Markov chain $M-C_{1}-C_{2}$ holds, i.e.,
$I(M;C_{2}|C_{1})=0$.

We emphasize that Definition 1 requires _all_ the information $C_{2}$ has
about $M$ must be available in $C_{1}$ in order to convey the intuitive idea
of message forwarding444We note that this interpretation is not accurate if
$C_{2}$ is an invertible and deterministic function of $C_{1}$, but is an
unlikely scenario in neural data.. We note that $M$-forward is a notion of
communication that is built on a dependence measure, unlike previous works[23,
24], whose notions of communication rest on correlational measures.

##  IV Experimental results

We apply our proposed methodology to data we collected to investigate activity
of two neural populations in the whisker-based sensory detection system in
mice. We show how IR, dPCA, mTDR, and DSAE applied to each population reveal
correlations across the trial. Then, we infer the absence of $M$-forwarding on
the dimension-reduced populations. We provide additional comparisons to PCA
and PLS in Appendix C.

###  IV-A Experiment design and data collection

We consider a simplified network of whisker-based sensory detection consisting
of the primary somatosensory cortex (whisker S1) and the superior colliculus
(SC). Whisker S1 is a specialized neural population in the mouse cortex that
processes tactile stimuli from whisker deflections [27]. SC is a sub-cortical
deep-brain region that is also involved in sensory processing [12]. Anatomical
evidence has shown that SC receives direct, bottom-up input from the brainstem
as well as top-down modulation from S1 [28]. The former is a direct
transmission of $M$ to SC and happens within 15ms of whisker stimulation,
while whisker tactile information likely reaches SC from S1 in the 15-40ms
window after stimulation[29]. However, it is unclear whether SC $M$-forwards
to S1 in this interval [30]. The network is shown in Fig. 1 on the left.

Figure 1: (Left) One whisker is deflected while population activity in S1 and
SC is recorded. We infer the existence of the red and blue arrows from S1 to
SC and SC to S1, respectively. (Right) Trial structure. The whisker is
deflected at 0s and the response window opens at 0.3s. The blue shaded region
depicts the window between 15-40ms after deflection in which communication is
posited to occur, and where we test for $M$-forwarding.

We conducted a pilot experiment on one mouse to understand how S1 and SC
interact with respect to activation from whisker deflections. The message $M$
is a stimulus taking discrete values that represent varying angles of whisker
deflection designed to invoke near-certain detection (large deflection angle),
ambiguous detection (medium to small deflection angle), and near-certain
absence of detection (zero deflection). The mouse is trained to respond as
accurately as possible by licking left or right, and receives a reward for
correct detection. Each trial consists of two key events: whisker deflection,
followed by licking response (Fig. 1 right). We investigate the possibility of
low-dimensional $M$-forwarding in the 15-40ms interval after deflection, as
evidence indicates that $M$ may be forwarded from S1 to SC within this period
[29]. Appendix B provides further details about the experiment and data
analysis.

###  IV-B Marginal $M$-relevance in S1 and SC

Fig. 2 illustrates how each method captures $M$-relevance in S1. IR, dPCA, and
mTDR capture similar $M$-relevant activity in their first dimension. DSAE is
the only method that is not message-dependent, and this is visibly apparent in
smaller and noisier correlations in each dimension. As in S1, IR captures
$M$-relevant activity in SC in its top dimension, highlighted in three clear
peaks (Fig. 3). Peak 1 is recognized as neural activity responsive to
stimulation [29]. Peaks 2 and 3 likely represent arousal and motor activity
related to licking [31]. The first mTDR dimension identifies a similarly
strong peak 2, but diminished peak 1 relative to IR. Peak 3 is evident only in
later dimensions (it does not appear as a distinguishable feature in the first
dimension). Peaks 1 and 2 are also present in dPCA’s first dimension, but are
noisier than in IR. Besides peak 1, DSAE shows almost no $M$-relevance in SC.

Due to the heterogeneous activity of neurons in SC [31], low-dimensional
correlations in SC exhibit higher variation relative to S1. While whisker S1
neurons are primarily activated by whisker deflection only, SC is more
functionally diverse. $M$-relevant activity in SC is not isolated to whisker
stimulation [31, 32], as indicated by the three peaks. IR demonstrates
robustness to this noise by maintaining clear $M$-relevance in its first
dimension. In Appendix C, we further show that IR reveals more nuanced
activity than the peristimulus time histogram (PSTH), which is another method
commonly used to visualize relevant population activity.

Figure 2: $M$-correlations of low-dimensional representations of S1 for the
top three dimensions, d1, d2, and d3. The solid vertical line at 0s marks the
time of whisker deflection and the dashed line at 0.3s marks response start.

###  IV-C Inferring low-dimensional dependencies

We present an application of IR for inferring the existence of a direct
pathway between $M$ and low-dimensional representations of S1 and SC. As
discussed in Sec. IV-A, we restrict our attention to the 15-40ms window after
deflection. First, we identify the top IR dimension for the S1 dataset,
$X_{S1}$, and the top IR dimension for the SC dataset, $X_{SC}$. We project
onto these linear dimensions, which yields the single ‘conceptual neuron’ $A$
of S1 and $B$ of SC that represent population-level activity most correlated
to $M$ in the relevant window. Then, we test for dependence between $M$ and
$A$, as well as $M$ and $B$ to verify that the dimension-reduced data still
retain relevance with the message [33]. Afterwards, we perform two tests of
conditional dependence [34] for both possible $M$-forward structures: $M-A-B$
and $M-B-A$. This analysis is repeated for dPCA, mTDR, and DSAE.

Because the null hypothesis of both tests is (conditional) independence,
rejecting the null of both independence and conditional independence implies
the existence of a _direct_ pathway from $M$. For example, for the data to
suggest a direct pathway from $M$ to $A$, we would need to reject the null
hypotheses $A\perp\\!\\!\\!\perp M$ and $A\perp\\!\\!\\!\perp M\;|\;B$. In
summary, rejection indicates that $M$ is _not_ forwarded from the other
conceptual neuron. The data are binned into 15ms intervals, requiring two
hypothesis tests to cover the window between 15-40ms after stimulation. We
test at a level of $\alpha=0.05$ and apply Bonferroni correction on the
estimated p-values. Table I summarizes the results from each hypothesis test,
where the top row lists the null hypothesis. $A$ is the dimension-reduced
representation of S1 and $B$ is the dimension-reduced representation of SC.
‘S’ means the p-value is significant (i.e., we reject the null hypothesis) and
(conditional) dependence can be concluded, while ‘NS’ means not significant
and the test is inconclusive. The p-values of each test are reported in Table
II in Appendix C.

TABLE I: Hypothesis tests in 100ms window after deflection Method | $A\perp\\!\\!\\!\perp M$ | $B\perp\\!\\!\\!\perp M$ | $A\perp\\!\\!\\!\perp M|B$ | $B\perp\\!\\!\\!\perp M|A$
---|---|---|---|---
IR | _S_ | _S_ | _S_ | NS
dPCA | _S_ | _S_ | _S_ | NS
mTDR | _S_ | NS | _S_ | NS
SAE | _S_ | NS | _S_ | NS
Figure 3: $M$-correlations of low-dimensional representations of SC for the
top three dimensions. Arrows 1, 2, and 3 in IR mark three distinct peaks
identified by IR.

The hypothesis tests using data pre-processed with IR allow us to confirm that
its reduced representations of S1 and SC indeed remain unconditionally
dependent with $M$ (consistent with visual observation from Fig. 3).
Furthermore, there exists a direct pathway from $M$ to $A$. Although this
conclusion is made on a low-dimensional representation of S1, it is consistent
with existing biological evidence [12]. There is not enough evidence for the
test $B\perp\\!\\!\\!\perp M|A$ to be conclusive, suggesting the possibility
of weak or nonexistent $M$-forwarding from SC to S1 in this window relative to
the evidence collected for S1. Hypothesis tests with dPCA reveal the same
structure as IR, but estimate a larger p-value than for IR (Appendix C). mTDR
and DSAE also detect (un)conditional dependence of S1 with $M$, but are unable
to conclude that SC depends on $M$, even unconditionally.

##  V Discussion

We propose a new framework for inferring low-dimensional communication between
two neural populations that identifies $M$-forwarding structures in IR-reduced
dimensions. IR may be viewed as a pre-processing step which effectively
defines noise relative to the message. Neural activity found to be irrelevant
to $M$ is removed, even if that activity is part of another biological
process. IR is one approach to $M$-relevant dimension reduction that is useful
when the message is known and measurable, a condition often satisfied in many
experimental designs. However, IR is not suitable for settings in which the
message is ambiguous, even when the question is one of interaction. For
example, following traumatic brain injuries, pathological areas of silence can
be ‘communicated’ to other brain areas, where the definition of the message
becomes unclear[35].

The alignment of IR conclusions with previous biological results suggests its
treatment of noise is effective for inferring an accurate summary of
communication. This highlights the importance of choosing an appropriate
measure of relevance to address the scientific question. In principle, our
framework could employ other measures of $M$-relevance, such as mutual
information or distance covariance [22], and could extend to nonlinear models
with kernel forms. Here, we compute $M$-forwarding using only one dimension.
Future work will investigate $M$-forwarding computed over a larger set of
dimensions.

## Acknowledgment

We thank Aswin Sankaranarayanan for helpful discussions.

## References

  * [1] N. A. Steinmetz, C. Koch, K. D. Harris, and M. Carandini, “Challenges and opportunities for large-scale electrophysiology with neuropixels probes,” _Current opinion in neurobiology_ , vol. 50, pp. 92–100, 2018.
  * [2] I. H. Stevenson and K. P. Kording, “How advances in neural recording affect data analysis,” _Nature neuroscience_ , vol. 14, no. 2, pp. 139–142, 2011.
  * [3] S. Kim, C. J. Quinn, N. Kiyavash, and T. P. Coleman, “Dynamic and succinct statistical analysis of neuroscience data,” _Proceedings of the IEEE_ , vol. 102, no. 5, pp. 683–698, 2014.
  * [4] G. Schamberg and T. P. Coleman, “Measuring sample path causal influences with relative entropy,” _IEEE Transactions on Information Theory_ , vol. 66, no. 5, pp. 2777–2798, 2019.
  * [5] C. J. Quinn, A. Pinar, and N. Kiyavash, “Bounded-degree connected approximations of stochastic networks,” _IEEE Transactions on Molecular, Biological and Multi-Scale Communications_ , vol. 3, no. 2, pp. 79–88, 2017.
  * [6] P. Venkatesh, S. Dutta, and P. Grover, “How should we define information flow in neural circuits?” in _2019 IEEE international symposium on information theory (ISIT)_. IEEE, 2019, pp. 176–180.
  * [7] ——, “How else can we define information flow in neural circuits?” in _2020 IEEE International Symposium on Information Theory (ISIT)_. IEEE, 2020, pp. 2879–2884.
  * [8] ——, “Information flow in computational systems,” _IEEE Transactions on Information Theory_ , vol. 66, no. 9, pp. 5456–5491, 2020.
  * [9] D. Kobak, W. Brendel, C. Constantinidis, C. E. Feierstein, A. Kepecs, Z. F. Mainen, X.-L. Qi, R. Romo, N. Uchida, and C. K. Machens, “Demixed principal component analysis of neural population data,” _eLife_ , vol. 5, p. e10989, 2016.
  * [10] M. C. Aoi, V. Mante, and J. W. Pillow, “Prefrontal cortex exhibits multidimensional dynamic encoding during decision-making,” _Nature neuroscience_ , vol. 23, no. 11, pp. 1410–1420, 2020.
  * [11] Y. Qiu, W. Zhou, N. Yu, and P. Du, “Denoising sparse autoencoder-based ictal EEG classification,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 26, no. 9, pp. 1717–1726, 2018.
  * [12] Y. K. Hong, C. O. Lacefield, C. C. Rodgers, and R. M. Bruno, “Sensation, movement and learning in the absence of barrel cortex,” _Nature_ , vol. 561, no. 7724, pp. 542–546, 2018.
  * [13] J. P. Cunningham and B. M. Yu, “Dimensionality reduction for large-scale neural recordings,” _Nature neuroscience_ , vol. 17, no. 11, pp. 1500–1509, 2014.
  * [14] J. K. Chapin and M. A. Nicolelis, “Principal component analysis of neuronal ensemble activity reveals multidimensional somatosensory representations,” _Journal of neuroscience methods_ , vol. 94, no. 1, pp. 121–140, 1999.
  * [15] C. M. Bishop and N. M. Nasrabadi, _Pattern recognition and machine learning_. Springer, 2006, vol. 4, no. 4.
  * [16] W. K. Härdle, L. Simar, W. K. Härdle, and L. Simar, “Canonical correlation analysis,” _Applied multivariate statistical analysis_ , pp. 443–454, 2015.
  * [17] A. Höskuldsson, “PLS regression methods,” _Journal of chemometrics_ , vol. 2, no. 3, pp. 211–228, 1988.
  * [18] S. Wold, A. Ruhe, H. Wold, and W. Dunn, Iii, “The collinearity problem in linear regression. The partial least squares (PLS) approach to generalized inverses,” _SIAM Journal on Scientific and Statistical Computing_ , vol. 5, no. 3, pp. 735–743, 1984.
  * [19] F. R. Bach and M. I. Jordan, “A probabilistic interpretation of canonical correlation analysis,” Department of Statistics, University of California, Berkeley, Technical Report 688, 2005.
  * [20] M. Aoi and J. W. Pillow, “Model-based targeted dimensionality reduction for neuronal population data,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [21] J. Semedo, A. Zandvakili, A. Kohn, C. K. Machens, and B. M. Yu, “Extracting latent structure from multiple interacting neural populations,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [22] B. Cowley, J. Semedo, A. Zandvakili, M. Smith, A. Kohn, and B. Yu, “Distance covariance analysis,” in _Artificial Intelligence and Statistics_. PMLR, 2017, pp. 242–251.
  * [23] E. Gokcen, A. I. Jasper, J. D. Semedo, A. Zandvakili, A. Kohn, C. K. Machens, and B. M. Yu, “Disentangling the flow of signals between populations of neurons,” _Nature Computational Science_ , vol. 2, no. 8, pp. 512–525, 2022.
  * [24] M. Ding, Y. Chen, and S. L. Bressler, “Granger causality: basic theory and application to neuroscience,” _Handbook of time series analysis: recent theoretical developments and applications_ , pp. 437–460, 2006.
  * [25] C. J. Quinn, T. P. Coleman, N. Kiyavash, and N. G. Hatsopoulos, “Estimating the directed information to infer causal relationships in ensemble neural spike train recordings,” _Journal of computational neuroscience_ , vol. 30, pp. 17–44, 2011.
  * [26] C. Pandarinath, D. J. O’Shea, J. Collins, R. Jozefowicz, S. D. Stavisky, J. C. Kao, E. M. Trautmann, M. T. Kaufman, S. I. Ryu, L. R. Hochberg _et al._ , “Inferring single-trial neural population dynamics using sequential auto-encoders,” _Nature methods_ , vol. 15, no. 10, pp. 805–815, 2018.
  * [27] J. F. Staiger and C. C. Petersen, “Neuronal circuits in barrel cortex for whisker sensory perception,” _Physiological reviews_ , vol. 101, no. 1, pp. 353–415, 2021.
  * [28] R. Aronoff, F. Matyas, C. Mateo, C. Ciron, B. Schneider, and C. C. Petersen, “Long-range connectivity of mouse primary somatosensory barrel cortex,” _European Journal of Neuroscience_ , vol. 31, no. 12, pp. 2221–2233, 2010.
  * [29] M. A. Castro-Alamancos and M. Favero, “Whisker-related afferents in superior colliculus,” _Journal of Neurophysiology_ , vol. 115, no. 5, pp. 2265–2279, 2016.
  * [30] S. Gharaei, S. Honnuraiah, E. Arabzadeh, and G. J. Stuart, “Superior colliculus modulates cortical coding of somatosensory information,” _Nature communications_ , vol. 11, no. 1, p. 1693, 2020.
  * [31] U. C. Drager and D. H. Hubel, “Physiology of visual cells in mouse superior colliculus and correlation with somatosensory and auditory input,” _Nature_ , vol. 253, no. 5488, pp. 203–204, 1975.
  * [32] J. Lee and B. L. Sabatini, “Striatal indirect pathway mediates exploration via collicular competition,” _Nature_ , vol. 599, no. 7886, pp. 645–649, 2021.
  * [33] A. Kraskov, H. Stögbauer, and P. Grassberger, “Estimating mutual information,” _Physical review E_ , vol. 69, no. 6, p. 066138, 2004.
  * [34] J. Runge, “Conditional independence testing based on a nearest-neighbor estimator of conditional mutual information,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2018, pp. 938–947.
  * [35] A. Chamanzar, M. Behrmann, and P. Grover, “Neural silences can be localized rapidly using noninvasive scalp EEG,” _Communications Biology_ , vol. 4, no. 1, p. 429, 2021.
  * [36] S. Seabold and J. Perktold, “Statsmodels: econometric and statistical modeling with Python,” in _9th Python in Science Conference_ , 2010.
  * [37] A. Yang, A. Ghassami, M. Raginsky, N. Kiyavash, and E. Rosenbaum, “Model-augmented estimation of conditional mutual information for feature selection,” _Proceedings of Machine Learning Research_ , 2019.
  * [38] H. Shimazaki and S. Shinomoto, “A method for selecting the bin size of a time histogram,” _Neural computation_ , vol. 19, no. 6, pp. 1503–1527, 2007.

## Appendix A

Regression is LMMSE and correlation-maximizing vector.
We show that the vector that maximizes sample correlation also achieves the
LMMSE estimate.

  * •

Let $\hat{M}=u^{T}X$ be a linear approximation of $M$ from $X$ for
$u\in\mathbb{R}^{n}$, where we assume $M$ and $X$ are zero mean. Then, the
LMMSE vector is given by

$\displaystyle v$
$\displaystyle=\operatorname*{argmin}_{u}\mathbb{E}[(\hat{M}-M)^{2}]$
$\displaystyle=\operatorname*{argmin}_{u}\mathbb{E}[(u^{T}X-M)^{2}].$ (7)

The solution to (7) is $v=\Sigma_{X}^{-1}\Sigma_{XM}$.

  * •

Let $r(u^{T}X,M)$ be the sample correlation coefficient between a linear
projection of $X$ and $M$. The optimization from (1) is re-expressed as

$\displaystyle v$ $\displaystyle=\operatorname*{argmax}_{u}r(u^{T}X,M)$
$\displaystyle=\operatorname*{argmax}_{u}\frac{u^{T}XM}{\sqrt{\sigma_{m}^{2}}\|u^{T}X\|}.$
(8)

Perform a change of variable by letting $c=X^{T}u$. Then, (8) becomes

$\displaystyle v$
$\displaystyle=\operatorname*{argmax}_{c}\frac{c^{T}M}{\|c\|}.$ (9)

The problem in (9) is optimized when $c=M$. Using this in our change of
variable expression, $v=(XX^{T})^{-1}XM=\Sigma_{X}^{-1}\Sigma_{XM}$.

Thus, IR can also be interpreted as finding the best average linear
approximation to $M$ at each iteration. Theorem 2 provides an additional
interpretation of IR vectors as the sufficient representation of $M$ when $X$
and $M$ are jointly Gaussian.

######  Theorem 2.

Suppose $M\in\mathbb{R}$ and $X\in\mathbb{R}^{n}$, and
$P(M,X)=\mathcal{N}(0,\Sigma)$, where $X$ has marginal covariance
$\Sigma_{X}$. If $\Sigma_{X}$ is full rank and
$v=\gamma\Sigma_{X}^{-1}\Sigma_{XM}$ for some constant $\gamma\neq 0$, then
$X$ projected onto $v$ is a sufficient representation of $M$. In other words,
the Markov chain $M-v^{T}X-X$ holds.

###### Proof.

We first define the relevant distributions: $P(M,X)$, $P(M,X,v^{T}X)$, and
$P(M,X|v^{T}X)$. Then, we derive the condition for sufficiency and show that
the condition holds for $v=\gamma\Sigma_{X}^{-1}\Sigma_{XM}$, i.e.
$I(M;X|v^{T}X)=0$.
For $P(M,X)$, define its covariance matrix as

$\displaystyle\Sigma=\begin{bmatrix}\sigma_{M}^{2}&\Sigma_{MX}\\\
\Sigma_{XM}&\Sigma_{X}\end{bmatrix}$

where $\Sigma_{X}\in\mathbb{R}^{n\times n}$ is the marginal covariance matrix
of $X$ and $\Sigma_{XM}\in\mathbb{R}^{n}$ is the cross covariance vector of
$X$ and $M$. Now, let $P(M,X,v^{T}X)=\mathcal{N}(\mu_{0},\Sigma_{0})$. The
covariance $\Sigma_{0}$ is

$\displaystyle\Sigma_{0}=\begin{bmatrix}\sigma_{M}^{2}&\Sigma_{MX}&\Sigma_{MX}v\\\
\Sigma_{XM}&\Sigma_{X}&\Sigma_{X}v\\\
v^{T}\Sigma_{XM}&v^{T}\Sigma_{X}&v^{T}\Sigma_{X}v\end{bmatrix}.$

Finally, we derive the conditional distribution $P(M,X|v^{T}X)$. Using
properties of multivariate Gaussian distributions, the covariance $\Sigma_{C}$
of the conditional distribution is

$\displaystyle\Sigma_{C}$
$\displaystyle=\begin{bmatrix}\sigma_{m}^{2}&\Sigma_{MX}\\\
\Sigma_{XM}&\Sigma_{X}\end{bmatrix}-\begin{bmatrix}\frac{\Sigma_{MX}v}{v^{T}\Sigma_{X}v}\\\
\frac{\Sigma_{X}v}{v^{T}\Sigma_{X}v}\end{bmatrix}\begin{bmatrix}v^{T}\Sigma_{XM}&v^{T}\Sigma_{X}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\sigma_{m}^{2}-\frac{\Sigma_{MX}vv^{T}\Sigma_{XM}}{v^{T}\Sigma_{X}v}&\Sigma_{MX}-\frac{\Sigma_{MX}vv^{T}\Sigma_{X}}{v^{T}\Sigma_{X}v}\\\
\Sigma_{XM}-\frac{\Sigma_{X}vv^{T}\Sigma_{XM}}{v^{T}\Sigma_{X}v}&\Sigma_{X}-\frac{\Sigma_{X}vv^{T}\Sigma_{X}}{v^{T}\Sigma_{X}v}\end{bmatrix}.$

Because zero correlation implies independence in jointly Gaussian
distributions, setting the off-diagonals of $\Sigma_{C}$ to zero implies
$I(M;X|v^{T}X)=0$. Without loss of generality, we derive the condition for
sufficiency using the top right matrix. This yields the condition:

$\displaystyle\Sigma_{MX}-\frac{\Sigma_{MX}vv^{T}\Sigma_{X}}{v^{T}\Sigma_{X}v}=0$
$\displaystyle\Leftrightarrow
v^{T}\Sigma_{X}v\Sigma_{MX}=\Sigma_{MX}vv^{T}\Sigma_{X}.$ (10)

When $v=\gamma\Sigma_{X}^{-1}\Sigma_{XM}$, the LHS of (10) equals the RHS:

$\displaystyle v^{T}\Sigma_{X}v\Sigma_{MX}$
$\displaystyle=\gamma^{2}\Sigma_{MX}\Sigma_{X}^{-1}\Sigma_{X}\Sigma_{X}^{-1}\Sigma_{XM}\Sigma_{MX}$
$\displaystyle=\gamma^{2}\Sigma_{MX}\Sigma_{X}^{-1}\Sigma_{XM}\Sigma_{MX}$
$\displaystyle=\gamma^{2}\Sigma_{MX}\Sigma_{X}^{-1}\Sigma_{XM}\Sigma_{MX}\Sigma_{X}^{-1}\Sigma_{X}$
$\displaystyle=\Sigma_{MX}vv^{T}\Sigma_{X}.$

∎

Proof of Theorem 1.

###### Proof.

The correlation coefficient is

$\displaystyle r(p_{i}^{T}X,M)=\frac{\sum_{j=1}^{T}\langle p_{i},x_{j}\rangle
M_{j}}{\bar{s}_{M}\sqrt{\sum_{j=1}^{T}\langle p_{i},x_{j}\rangle^{2}}},$

where $x_{j}$ is the $j$th column of $X$ and

$\displaystyle\bar{s}_{M}=\sqrt{\sum_{j=1}^{T}M_{j}^{2}}.$

The correlation coefficient is expressed in matrix form as

$\displaystyle r(p_{i}^{T}X,M)=\frac{p_{i}^{T}XM}{\bar{s}_{M}\|p_{i}^{T}X\|}.$

The regression vector is $v=(XX^{T})^{-1}XM$. Hence,

$\displaystyle XM=XX^{T}v=(N-1)\Sigma_{X}v,$

where $\Sigma_{X}$ is the sample covariance matrix, expressed in terms of $X$
as

$\displaystyle\Sigma_{X}=\frac{XX^{T}}{N-1}.$

Note that the denominator of the correlation coefficient is
$\|p_{i}^{T}X\|=\sqrt{(p_{i}^{T}X)(p_{i}^{T}X)^{T}}=\sqrt{(N-1)p_{i}^{T}\Sigma_{X}p_{i}}$.
Making the appropriate substitutions for $XM$ and $\|p_{i}^{T}X\|$, the
correlation coefficient can be simplified as

$\displaystyle r(p_{i}^{T}X,M)$
$\displaystyle=\frac{(N-1)p_{i}^{T}\Sigma_{X}v}{\bar{s}_{M}\sqrt{(N-1)p_{i}^{T}\Sigma_{X}p_{i}}}$
$\displaystyle=\frac{(N-1)\lambda_{i}\langle
v,p_{i}\rangle}{\bar{s}_{M}\sqrt{(N-1)\lambda_{i}\langle p_{i},p_{i}\rangle}}$
$\displaystyle=\frac{\sqrt{\lambda_{i}(N-1)}}{\bar{s}_{M}}\langle
v,p_{i}\rangle,$

where the second equality is due to $p_{i}$ being an eigenvector of
$\Sigma_{X}$ with singular value $\lambda_{i}$. The first IR dimension is
$v^{*}=v/\|v\|$. Projecting with respect to $v^{*}$ instead of $v$, the
expression simplifies as

$\displaystyle r(p_{i}^{T}X,M)$ $\displaystyle=\alpha\sqrt{\lambda_{i}}\
\langle v^{*},p_{i}\rangle,$

where $\alpha=\|v\|/s_{M}$ and $s_{M}$ is the standard deviation of $M$. ∎

##  Appendix B Experiment Methodology & Data Analysis

Experiment methodology. The study protocols were approved by the CMU
Institutional Animal Care and Use Committee. Neural spiking activity in one
mouse was recorded at 1kHz in a sub-sampled population of neurons in S1
($n=38$) and SC ($n=35$) over a total of 208 trials. Each trial is a 5s
recording in which the whisker is deflected (i.e., stimulation). The mice are
trained to respond after a delay period of 300ms following stimulation to
indicate whether it detected a deflection (right lick) or no deflection (left
lick) for a water reward. In this paper, we consider the 1s interval between
200ms pre-stimulation and 800ms post-stimulation.

The whisker deflection values are from the set $\\{0,3,4,6,10\\}$. The value 0
means no deflection, 10 means max deflection, and the remaining values are
designed to induce ambiguous detection ability. Deflection values $0,3,4,6$
and $10$ are randomly presented with probabilities 0.4, 0.1, 0.1, 0.2, and
0.2, respectively. This is to include enough rewarded (i.e. unambiguous)
trials to keep the mouse engaged for the duration of the experiment. Figure 1
(left) was created with BioRender.com.

Data analysis: pre-processing. Unresponsive neurons are first removed, where a
neuron is considered unresponsive if it produces less than 15 spikes in the
100ms period after whisker stimulation over all trials with deflection value
greater than 0. A total of 31 neurons in S1 and 25 neurons in SC remain after
removing unresponsive neurons. The trials are binned into 15ms time bins. We
perform 4-fold cross validation so each fold tests on 25% of the dataset.
Figs. 2, 3, and 5 show average correlations found on the test data over all
folds.

Data analysis: model parameters. We train mTDR over all time points as this is
its model specification. For the remaining methods, we train a new model at
each time point. Methods with additional special parameters are summarized
below:

  * •

_IR_. The regression vectors are found using the statsmodels OLS function in
Python [36].

  * •

_dPCA_. The model is trained in Python with stimulus and time parameters, and
without regularization. Three components are identified and their weights from
the stimulus-only set are extracted.

  * •

_mTDR_. The optimal rank (i.e. number of dimensions) is found by mTDR’s greedy
estimation procedure. This finds 6 dimensions in S1 and 4 dimensions in SC.
The top 3 representations are extracted. Two covariates are used to model the
one-dimensional stimulus $M$ and a bias term.

  * •

_DSAE_. The model closely follows the network given in [11]. The number of
units in the input layer is the number of responsive neurons, which varies
between S1 and SC. There is only one hidden layer, which consists of three
units. All activation functions are sigmoid functions. Our loss function is
the mean-square error with a KL divergence loss which encourages sparsity. The
KL sparsity parameter is set to $\rho=0.1$. The penalty weight on the KL loss
is set to $\beta=0.01$. We apply a corruption operation to achieve a denoising
autoencoder. The noise distribution is Poisson with rate parameter set for
each input unit to the mean firing rate of the neuron in a given time bin. The
noise ratio is set to 10%, meaning 10% of the input nodes are corrupted by
Poisson noise. The model is trained over 300 epochs with a learning rate of
0.01. As there is no natural ordering of the units, we choose to order them by
total absolute correlation with $M$ over the whole trial.

Data analysis: hypothesis tests. Both independence and conditional
independence tests are based on projections obtained by training over the
entire dataset in the relevant time window between 15-40ms after stimulation.
For mTDR, we train the model over the whole trial before picking out
projections from the relevant window. We use the pycit implementation (based
on [37]) of k-nearest neighbors dependence and conditional dependence testing,
drawing 10,000 samples from each population and using 5 neighbors.

##  Appendix C Additional comparisons

Peristimulus time histogram (PSTH). The PSTH is a method commonly applied to
spiking data for visualizing key events over the trial [38]. A PSTH is
computed by summing the total number of times a neuron, or population of
neurons, spikes in a given time bin. Fig. 4 shows the PSTH computed over all
trials and all neurons in S1 and SC. Aside from peak C in the PSTH of SC, no
other peaks are immediately evident, in contrast to peaks 1, 2, and 3 in IR
(Fig. 3). This suggests that correlations in IR-projected data can reveal
$M$-dependent neural activity that may not be evident in the PSTH. Another
feature in Fig. 4 is peak B in S1. This is not evident in any of the dimension
reduction methods in Fig. 2. Rather, the correlations of the dimension-reduced
data suggest that the entire period between the first and second peaks is
highly correlated with the message. However, this may also highlight a
limitation in using correlation as the measure of $M$-relevance, which is
limited to capturing linear dependencies.

Figure 4: PSTH for S1 (left) and SC (right) over all trials. The S1 PSTH
exhibits two distinctive peaks, labeled A and B, while the SC PSTH only has
one evident peak, labeled C.

PCA and PLS comparison. We provide additional comparisons with PCA and PLS. In
S1, the first dimension of both PCA and PLS identify a strong, decaying
correlation with $M$, and minimal correlations in the other two dimensions.
PLS identifies a longer period of $M$-relevance than PCA. In SC, PLS
identifies the three peaks of $M$-relevance similarly to IR in Fig. 3, albeit
noisier. The only visibly identifiable feature in PCA is peak 1 at 0s. The
remaining dimensions of PCA only detect small correlations.

Figure 5: $M$-correlations of low-dimensional representations of S1 and SC
using PCA and PLS.

Table of p-values.

TABLE II: p-values for hypothesis tests Method | $t$ | A$\perp\\!\\!\\!\perp M$ | B$\perp\\!\\!\\!\perp M$ | A$\perp\\!\\!\\!\perp M\;|\;$B | B$\perp\\!\\!\\!\perp M\;|\;$A
---|---|---|---|---|---
IR | $p_{1}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $0.0027$
$p_{2}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $0.0364$
dPCA | $p_{1}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $0.0374$
$p_{2}$ | $<10^{-4}$ | $0.0023$ | $<10^{-4}$ | $<10^{-4}$
mTDR | $p_{1}$ | $<10^{-4}$ | $0.0004$ | $<10^{-4}$ | $0.0101$
$p_{2}$ | $<10^{-4}$ | $0.4376$ | $<10^{-4}$ | $0.9801$
DSAE | $p_{1}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $0.1551$
$p_{2}$ | $10^{-4}$ | $0.1723$ | $<10^{-4}$ | $0.0004$
PCA | $p_{1}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $0.0159$
$p_{2}$ | $<10^{-4}$ | $0.7223$ | $<10^{-4}$ | $0.5678$
PLS | $p_{1}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $0.0023$
$p_{2}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$ | $<10^{-4}$

Table II summarizes the p-values for the two time bins covering the 15-40ms
interval after stimulation. Column 2 of the top row indicates the time bin and
columns 3-6 list the null hypothesis. $A$ and $B$ are the dimension-reduced
representations of S1 and SC, respectively. Since we test at a level of
$\alpha=0.05$ with Bonferroni correction, we reject the null hypothesis if
each p-value is less than 0.025.

PCA identifies a direct pathway from $M$ to S1, but is inconclusive regarding
dependence of SC and $M$. PLS is the only method that attains significance in
all four hypothesis tests. This result means the low-dimensional PLS
representation identifies a direct pathway from $M$ to S1 and $M$ to SC. While
this functional structure is biologically plausible, direct connectivity from
$M$ to SC has been found to occur within 15ms of the stimulus [29], which is
outside of the tested window.