Gaussian Partial Information Decomposition: Bias Correction and Application to High-dimensional Data
Authors: Praveen Venkatesh
Corbett Bennett
Sam Gale
Tamina K. Ramirez
Greggory Heller
Severine Durand
Shawn R Olsen
Stefan Mihalas
Venue: NeurIPS 2023 spotlight (2023)
Keywords: partial information decomposition
estimation
bias
inter-area interaction
neuroscience
TL;DR: We present a new bias-corrected estimator of the BROJA-PID for multivariate Gaussian distributions, show that it agrees with ground truth even at high dimensionalities, and use it to characterize inter-areal interactions in neural data

Abstract: Recent advances in neuroscientific experimental techniques have enabled us to simultaneously record the activity of thousands of neurons across multiple brain regions. This has led to a growing need for computational tools capable of analyzing how task-relevant information is represented and communicated between several brain regions. Partial information decompositions (PIDs) have emerged as one such tool, quantifying how much unique, redundant and synergistic information two or more brain regions carry about a task-relevant message. However, computing PIDs is computationally challenging in practice, and statistical issues such as the bias and variance of estimates remain largely unexplored. In this paper, we propose a new method for efficiently computing and estimating a PID definition on multivariate Gaussian distributions. We show empirically that our method satisfies an intuitive additivity property, and recovers the ground truth in a battery of canonical examples, even at high dimensionality. We also propose and evaluate, for the first time, a method to correct the bias in PID estimates at finite sample sizes. Finally, we demonstrate that our Gaussian PID effectively characterizes inter-areal interactions in the mouse brain, revealing higher redundancy between visual areas when a stimulus is behaviorally relevant.

Gaussian Partial Information Decomposition: Bias Correction and Application to High-dimensional Data
Praveen Venkatesh Corbett Bennett', Sam Gale', Tamina Ramirez',
Greggory Heller Séverine Durand', Shawn Olsen', Stefan Mihalas
Allen Institute; Uniivessity of Washington; 'Columbia University;
Masssacuusetts Institute Technology
praveen. "praveen. www.weeeeeaaalaaeeiniiccccocaaee org; Stefanm@@ll1enneee org
Abstract
Recent advances in neuroscientific experimental techniques have enabled us simultaneously record the activity of thousands of neurons across multiple brain regions. This has led to growing need for computational tools capable of analyzing how task-relevant information is represented and communicated between several brain regions. Partial information decompositions (PIDs) have emerged one such tool, quantifying how much unique, redundant and synergistic information two or more brain regions carry about task-relevant message. However, computing PIDs is computationally challenging in practice, and statistical issues such as the bias and variance of estimates remain largely unexplored. In this paper, propose a new method for efficiently computing and estimating PID definition on multivariate Gaussian distributions. We show empirically that our method satisfies an intuitive additivity property, and recovers the ground truth in a battery of canonical examples, even at high dimensionality. We also propose and evaluate, for the first time, a method to correct the bias in PID estimates at finite sample sizes. Finally, we demonstrate that our Gaussian PID effectively characterizes inter-areal interactions in the mouse brain, revealing higher redundancy between visual areas when a stimulus is behaviorally relevant.
1 Introduction
Neuroscientific experiments are increasingly collecting large-scale datasets with simultaneous record- ings of multiple brain regions with single-unit resolution [1-3]. These experimental advances call for new computational tools that can allow us to probe how multiple brain regions jointly process relevant information in behaving animal.
Partial Information Decompositions (PIDs) offer a new method for studying how different brain regions carry task-relevant information: they provide measures to quantify the amount of unique, redundant and synergistic information that one region has with respect to another. The information itself could pertain to task-relevant variables such as stimuli, behavioral responses, or information contained in third region. For example, we may be interested in how much information about stimulus is communicated or shared (i.e., redundantly present) between two brain regions over time. Or, we might be interested in the extent to which one region's is activity uniquely explains that of another, while excluding information corresponding spontaneous behaviors.
Ideas such as redundancy and synergy have long history in neuroscience, having been proposed for understanding noise correlations [4] and to understand differences encoding complexity between different brain regions [5]. PIDs have also been suggested for quantifying how much sensory information used execute behaviors [6] and for tracking stimulus-dependent information flows between brain regions [7, 8]. Outside of neuroscience, PID has been used to understand interactions
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

between different variables in financial markets [9], to quantify the relevance of different features for the purpose of feature selection machine learning [10], and to define and quantify bias in the field of fair Machine Learning [11].
An important constraint that has limited the broader adoption of PIDs in neuroscience is the com- putational difficulty of estimating PIDs for sigh-dimensional data. Many PID definitions that are operationally well rell-motivated involve solving an optimization problem over a space of probability distributions: the number of optimization variables can thus be exponential in the number of neu- rons [12]. This has led to the use of poorly motivated PID definitions that are easy to compute (such as the "MMI-PID" of [13], in works such as [9, 14-16)), or limited analyses to very few dimensions [17]. Furthermore, due to the limited exploration of estimators for PIDs, issues such as the bias and variance of estimates have received no attention far, to our knowledge.
In this paper, we make the following contributions:
1. We provide a new and efficient method for computing and estimating a well-known PID definition called the PID or the BROJA-PID [18] on Gaussian distributions (Section 3). By restricting our attention to Gaussian distributions, we are able to significantly reduce the number of optimization variables, so that this is just quadratic in the number of neurons, rather than exponential.
2. We present set of canonical examples for Gaussian distributions where ground truth is known, and show that our method outperforms others (Section 4).
3. We also raise (for what we believe is the first time) the issue of bias in PID estimates, propose method for correcting the bias, and empirically evaluate its performance (Section 5).
4. Finally, we show that our Gaussian PID estimator closely agrees with ground truth, even on non-Gaussian distributions, and show an example of its on real neural data (Section 6).
Related work
Our method is based on carlier work [12], where we also examined PIDs for Gaussian distributions. Our current work differs in a few key aspects: (i) we estimate the PID of a different PID definition, the -PD rather than the d-PID, because the 8-PID does not satisfy a basic property called additivity [19] (defined in Sec. 2); (ii) our current method provides an exact upper bound to the PID definition being computed, rather than approximate upper bound: (iii) we now consider the problem of estimation, not just computation, and explore the issue of the bias of PID estimates; and (iv) our current method much faster, and we demonstrate agreement with ground truth at much higher dimensionality.
Several other studies have also considered methods for efficiently estimating PIDs: Banerjee et al. [20] and Makkeh al. [21] address computing discrete PIDs, but their method does not scale to higher dimensions; Pakman et al. [17] estimate PIDs using copulas, but their method would also potentially be computationally prohibitive high dimensionalities: Liang et al. [22] use convex optimization to directly estimate the -PDI for general high-dimensional distributions, but they do not compare with ground truth at high dimensionality or examine bias in their estimates.
Background: An Introduction to PIDs and the -PD
In this section, we provide an introduction to the concept of partial information decomposition along with an illustrative example. Let M, X and Y be three random variables with joint distribution PMXY. PID decomposes the total mutual information between the message M and two constituent random variables X and Y into sum of four non-negative components that satisfy [18, 23]:
TX; Y)) = UT(M: X \ Y) + UT(M : Y \X) + RT(M : X;Y) + ST(M X:Y)
(1)
I(M; UI(M: X (\ Y) + RIMM : X;
(2)
I(M;Y)- = UIM Y X) + RI : X;Y)
(3)
Here, [AA; BB is the Shannon mutual information [24] between the random variables A and B, and the four terms in the RHS of (1) are respectively the information about M that is (i) uniquely present in X and not in Y: (ii) uniquely present in Y and not in X; (ii) redundantly present in both X and Y and can be extracted from either; and (iv) synergistically present in X and in Y. i.e., information which cannot be extracted from either of them individually, but can be extracted from their interaction. For the sake of brevity, we may also refer to these partial information components as UTx Uly RI and ST respectively. Notwithstanding notation, they should properly be understood to functions of the joint distribution PAXY.
2

---

Now, Ulx, Uly, RI and SI consist of four undefined quantities, subject to the three equations in (1)- In addition, they are typically assumed to be non-negative, RI and SI are each constrained to be symmetric in X and Y, and the functional forms of UTx and UTy should be identical when exchanging X for Y. Despite the number of constraints, many definitions satisfy all of them, each differing in its motivation and interpretation [18, 23, 25-27] (see [27, 28] for a review), and we need to formally define one of these partial information components to determine the other three.
Example 1. Before we jump into a specific definition, we provide intuition into what these terms mean using simple example, Suppose M = [A, B, C), X = [A, B, C and Y - [B, Z], where A, B, C, Z ~ i.i.d. Ber(O.5)) Then, X has 1 bit of unique information about M, i.e., A; Y has no unique information; X and Y both have 1 bit of redundant information, i.c., B, since it can be obtained from either X or Y; and X and Y have bit of synergistic information, i..., C., which cannot be obtained from either X or Y individually (since CΦZ Z 1 C), but can only be recovered when both X and Y are known. For more examples on binary variables, refer the reader to [18]. O
In this manuscript, we consider a definition that we refer to as the PDD [18, 25J, which is defined below. We chose to build an estimator for this definition for two reasons: (i) it is Blackwellian PID definition, i.e,, it has well-defined operational interpretations based on concepts from statistical decision theory (e.g., see [18, 27, 29] for details); and (ii) it satisfies many desirable properties (e.g.. see [18, 30D, and in particular, a property that we call additivity of independent components.
Definition 1 (~~PID [18)). The unique information about M present in X and not in Y is given by UI(M: : X x)) min QEap
Xx(M;X] Y), |Y),
where Ap = {QMXY: QMX = PMX QMY = }My and lo( 1-) is the conditional mutual information over the joint distribution QMXY. The remaining PID components, UI(M: : Y Y X).
RI(M: X; and SI(M Y.; follow from equations (1)(3).
Property 1 (Additivity of independent components). Suppose M = [M M2) X = x= X2], and
Y= [Y,, Y2), such that (M1, X1, Y1) 1 (M2, X2, Y2). Then, additivity implies that
UI(M: : X (Y)= UIM X2 \8) + UI(M2 X2 \ Y2),
(5)
and similarly for the other three partial information components, Uly, RI and S1.
Property 1 stipulates that we should be able to compute the PIDs of two independent systems separately, and then add the components across both systems. In effect, additivity implies that the PID of an isolated system should not depend on the PID of another isolated system, making it an intuitive and highly desirable property (see App. A.1 for concrete examples). Of the many PID definitions examined by Rauh et al. [19], only the P-II satisfied additivity (as proved in [18)..
3 Computing the PID for Gaussian Distributions
The first contribution of this paper a method to efficiently compute bounds on the PPDI for jointly Gaussian random vectors M, X and Y. To be precise, our method computes an upper bound for Ulx and Uly, and lower bounds for Rl and SI. Similar to our earlier work [12], we present a new PID definition that we call the ~PI which characterizes an upper bound on the unique information of the PPD by restricting the optimization space to jointly Gaussian QMXY:
Definition 2 (~G-PID). Let PMXY be jointly Gaussian. Then, the unique information about M present in X and not in Y is given by
Ulc(M:X(M)
min 12(M;X QEar
(6)
where Ap = {QMXY: QMXY jointly Gaussian, QMx = PMX: QMY = ParY) and lo is the conditional mutual information over the joint distribution QMXY.
If the optimal QmXY in the unrestricted optimization of Definition I happens to be Gaussian for some PMXY, then the PPD would be identical to the PID for that PMXY. We conjecture
ii..d. stands for "independent and identically distributed": X 1 Y means X and Y are independent. "This PID is also sometimes referred to as the BROJA PID (after the authors of [18]) or the ciinimum-synergy PID in the literature. We prefer to use an agnostic nomenclature as introduced in our earlier work [12]. because this PID was also introduced contemporaneously by [25].
3

---

that this happens whenever PMXY is Gaussian: for example, in a similar optimization problem for computing the information bottleneck [31], the optimal distribution is Gaussian whenever P is Gaussian [32, 33]. We provide empirical evidence in favor of this conjecture through sequence of examples with Gaussian in Sec. 4, where we recover the ground truth even when Q restricted to be Gaussian. However, we leave theoretical examination of this conjecture for future work.
In practical terms, restricting the search space Gaussian QMXY reduces the number of optimization variables from being exponential in the dimensionality to quadratic (see Appendix A.2), allowing us to compute the ~GDPID for much higher dimensionalities of M, X and Y. what follows, we show how the optimization problem for the ~PPID can be written out in closed-form and then solved using projected gradient descent.
3.1 Notation and Preliminaries
Suppose M, X and Y are jointly Gaussian random vectors of dimensions das, dx and dy respectively, with joint covariance matrix given by ΣMXY. We will make extensive use of the submatrices of ExxY, so we explain their notation here:
Σxy will denote the (dx +dy) x (dx++dy) joint (auto-)covariance matrix of the vector [XT, EXTYYT YTIT.
Ex,y (note the comma) will denote the dx x dy cross-covariance matrix between X and Y.
Exy.M will denote the (dx dy) dM cross-covariance matrix between the concatenated vector [X",Y]] and the vector M.
In general, groupings of vectors without commas represent joint covariance, while a comma represents cross-covariance between the groups on either side of the comma. The same notation will also be used for conditional covariance matrices: for example, M |M is the conditional joint covariance of (X, Y) given M, while 1 ur is the conditional cross-covariance between X and Y given M.
We will also use equivalent notation for the joint distribution [12], where PMXX is parameterized as a "broadcast channel" [24, Ch. 15.6] from to X and Y: X H M + Nx and Y M++++
Here, Hx == Σx.M and Hy == Σy.M represent channel gain matrices, while Nx and Ny represent additive noise and are not necessarily independent of each other: ANJ] (0, Exyjar)
Remark 1. Without loss of generality, we can assume that M, X and Y are all zero-mean, and that Σy = I. Further, we explicitly assume that the X and Y channels are individually whitened, i.c., that ExIa I and MyJas = I. This assumption precludes deterministic relationships between M and or Y, and is required to ensure that information quantities remain finite [12].
3.2 Optimizing the Union Information
Bertschinger et [18] showed that the minimizer for the unique information also the minimizer for the "union information", T X;Y)= == UTx + UTy + RT. In other words, we can also solve the following optimization problem, which yields simpler expressions for the objective and gradient:
TM : X;Y) = min Io(M;X,Y) s.t. QMx = PMx. QMY = PMY
(8)
QuXY
Now, suppose PM is Gaussian with covariance Dxrxy and the solution QM is also assumed to be Gaussian with covariance XQ MXY- Then, the constraint in (8) implies that 2M E and 1 Ety In other words, ΣM. Σx, Σr, and ΣM,xy are all constant across P and Q. Therefore, the only part of وº 1 that is variable is Mx,y, or equivalently, Nº " follows, we will drop the superscripts denoting the distribution, as this will be clear from context. Generally speaking, we will discuss the optimization problem and thus the distribution will be Q.
what
Proposition 1. The union information for the ~GPPID of Definition 2 is given by
18 = min log det (I + EM TT E-1 XXyIM EXYM
s.t. EXYIM T
(9)
FXX.YIM
2
'e can Σx.y as in place of Σx,y because they differ by constant: ΣX.YIM- Ex,y off-diagonal block in Exy EXYIM, which is equal to M'y, 2dd T MXY.M. which is constant across P and Q.

---

where the optimization variable Σx is an off diagonal block embedded within EXvIM; all other matrices in the objective are constants that are derived from > XY-
We solve the above optimization problem using projected gradient descent: we analytically derive the gradient and the projection operator for the constraint set as shown below. Then, we use the RProp algorithm [34, 35] for gradient descent, which independently adjusts the learning rates for each optimization parameter (derivations, and details of implementation and complexity are in App. A). Code for our implementation is available on GitHub [36], and details of the compute configuration given in App. E.
Proposition 2. The objective in Proposition I can be simplified to
f(Ex.rIM) = log det [ + H + B), 2 where B = (Hx Exy Hy) and S = (I- (I-ExyIIM* XX.II
(10)
The gradient of the objective with respect to yx,Y|I is given by
(E = S-tBBI++ + H Hy + A projection operator on to the constraint set ΣxyIM > can be obtained as follows: let E Σx = VAVT be the eigenvalue decomposition of ΣXYMM, with A =: diag(X.). Let X == ma X,) represent the rectified eigenvalues, and A = diag(A, (A,). Then, define
-H+)-
(11)
ΣxyIM VÄVT,
(12)
y Mrof X.YIM = 1--1/25 IX 5/1
E X YIM
(13)
where E EX|M. EvJM and Ex are submatrices of E Ex
4 Canonical Gaussian Examples
In this section, we show how well our PPD estimator performs on a series of Gaussian examples of increasing complexity, which have known ground truth. Barrett [13] showed that, for Gaussian distributions, the ~-PID reduces to the MMI-PID (defined below), whenever M is scalar. These also happen to be cases when the optimal distribution Q is Gaussian [12], and thus the ~G-PID should recover the ground truth. We then leverage additivity (Property 1) to combine two or more simple examples into complex ones, where ground truth continues to be known.
Definition 3 (Minimum Mutual Information (MMI) PID). Let the redundant information be defined as the minimum of the two mutual informations:
RIuM M X; = (I)
(14)
The remaining MMI-PID components, M (M : X \ Y), U MM Y X) and S (M
X; Y). follow from equations (1)-(3).
We first provide a Gaussian analog of Example in Examples 2-4 (for dı = dx = dy 1). We will use the channel notation described in Equation (7). Complete derivations for these examples (and some nuances that are omitted here) are presented in App. B.
Example 2 (Pure uniqueness: variable A from Example 1). Suppose M ~ (N, 1), Hx = 1 and Hy = O. with Nx, Ny i.i.d. (0, 1). Here, only X receives information about M, while Y is pure noise. Thus, X has unique information about M (UIx = X(M > 0), with no unique information in Y, and no redundancy or synergy (Uly = RI SI = 0).
0
Example 3 (Pure redundancy: variable B from Example 1). Ideally, we would set M ~ N 1). X M and Y = M. However, for continuous random variables, T(M; X) = de when M X.
So instead, we set M ~ NN(0, 1), Hx = 1 and Hy = 1, with Nx ~ ~N(0, 1) while Ny- = Nx (i.e..
X x=Y- Y, so they are both the same noisy version of M). In this case, X and Y are fully redundant since they both contain exactly the same information about M. Thus, RI 1(M; (X, > 0.
while UIx = U = SI 0.
Example 4 (Pure synergy: variable C' from Example 1). We cannot replicate pure synergy for Gaussian variables, but we can approach it a limit, Let M ~ (, 1), Hx = and Hy 0, with Nx ~ NOO o²) and Ny Nx (i.e., X = M+Y). Further, let o² - co. In this case, I(M; Y) 0
5

---

PIDs for Canonical Gaussian Examples
Unique and Redundant
Unique and Synergistic
Redundant and Synergistic
PID component
Partial information (bits)
K[M; Y3)
2.5
TTx
2.0
R
1.5
1.0
PID definition
0.-
MMI-PID Ground truth
0.00
2.00
1.75
3.16
5.62
10.00
17.78
31.62
58.23
100.00
0.00
0.10
0.20
0:70
0.00
0.90
0.99
0.00
0.30
0.50
0.99
Noise correlation,
Noise correlation,
Noise given ay
Figure 1: PID values for Examples 5, 6 and 7. The - and the Pc-PID agree exactly with the MMI-PID, which is known be the ground truth, since M scalar [13].
Increasing Gain in X,
Error in gain example
Rotation of X w.r.t.
PID component H[M, (X,
2.5
Partial information (bils
10
Partial informationn (bils)
2.0
10
1.50
27ty
1.5
1.00
1.0
10-41
definition
10
fp-PDD
0.25
MMI-PID
0.0
0.00
0:33
0.57
2.00
2.30
0.00
0:00
0.30
020
157
Gain
Angle,
Figure 2: Left and right panels respectively show PID values for Examples 8 and which combine two scalar examples with known ground truth, using Property 1. The middle panel shows the absolute error between each PID definition and the ground truth. The ~G-PID diverges from the dØ- and MMI-PIDs, and is the only one that agrees with the ground truth. The Pg-PID maintains an error less than 10-7 bits, whereas other definitions have errors greater than 0.1 bits.
and [[M; X) + 0 as d² - oc, so X and Y individually convey little to no information about M. However,
recover information about M from X and Y together by taking their difference, since X - Y = M. Thus, SI > O, while Uly RI = 0 and Ulx +
Examples 2, and 4 have been provided solely for intuition. Their PIDs can be inferred directly from Equations (1)-3). We next describe three one-dimensional examples that each have two non-zero only provide brief description and defer details to App. B.
PID components. For lack of space, We estimate the PPPI (as well as the PPD [12] and the ground MMI-PID [13]) for these examples and show that all three are equal (see Fig. 1).
Example 5 (Unique and redundant information). Let X be a noisy representation of M, and let Y be noisy representation of X with standard deviation øy When Y X (zero noise), this example reduces to Example 3. As Tyix + RI reduces while Ulx approaches X).
Example 6 (Unique and synergistic information). Let M ~ (N,, 1). Hx = 1, Hy = 0 and -O( o²) such that their correlation is p. When o² is finite and p 0, this example reduces to Example 2, since there can be no synergy between X and Y. As p + 1. Y + M: so the total mutual information T(M; (X, Y)) - oo, driven by synergy growing unbounded, while the unique component remains constant at X):
Ny Ny
L
= Example 7 (Redundant and synergistic information). Let M N 1), Hx = = and (0. 1) such that their correlation is p. When p < 1, X)M;X) and I(M;Y) are both equal by symmetry, and thus equal to RI (see Def. 3 for the MMI-PID, which ground truth here). As o reduces, the two channels X and Y have noisy representations of M with increasingly independent noise terms. Averaging the two, (X + Y)/2, will provide more information about M than either one of them individually (i.c., synergy), and thus SI increases as p reduces.
Ny Ny ~
D

---

Stability of the PeDPID over increasing dimensionality
I[M;; X))
Ulx
ST
Timing analysis
Partial Information (bits)
1024
taken (seconds)
PID
10
PID
128
10
0.25
10-2
Dimension ef and Y: dr
Figure 3: The first five plots on the left show PID values for Example 10. Different shadings represent different values of gain in X1 (a) from Example 8. The ~G-PID doubles every time d doubles as seen by the constant 45 slope on the base-2 log-log plot, even when das dx = dy 1024. The right-most plot shows a timing comparison between the ~G- and the 8c-PIDs on this example. The SG-PID failed for d 64; the ~G-PID's timing performance up to d 1024 is shown App. B.4.
The next set of examples will use the examples presented above in different combinations. This ensures that, where possible, the ground truth remains known in accordance with Property 1. These examples are also designed to reveal the differences between the ~.PID, the MMI-PID and the 5-PID: in particular, they show how the MMI-PID and the 5-PID fail where the ~-PID does not. These examples use -wo-dimensional M, X and Y, i.e., (dM, dx, dy) (, A diagrammatic representation of Examples 8 and 9 is given in App. B.2.
Example 8. Let X1 a&1 + Nx,1, Y1 = M1 + Ny.1, X2 M2 + Nx,2 and Y2 3M2₂ Ny.2, where M1, M2, Nx,,, Ny.; ~ i..dd (O0, i 1,2. Here, (M1, V1) is independent of (M2, X2, Y2), therefore using Property 1, we can add the PID values from their individual decompositions (which each have known ground truth via the MMI-PID since Mi and M2 are scalar). Fig. 2(1) compares the ~P-I the PGDD and the MMI-PID for the joint decomposition of [(;; (X, Y)). at different values of the gain in XI. Only the -PD matches the ground truth, as it the only definition here that is additive.
Example 9. Let M and Y be as in Example 8. Suppose X = Hx R(0) M, where Hx is a diagonal matrix with diagonal entries and 1, and R(0) is a 2 2 rotation matrix that rotates M by an angle O. When 0 = 0, X has higher gain for M1 while Y has higher gain for M2. When 0 increases to x x/2, X and Y have equal gains for both Mj and M2 (barring a difference in sign). Since X1, Yi) ) is not
independent of (M2, X2, Y2) for all 0, we know the ground truth only at the end-points. Nonetheless, the example shows a difference between the three definitions, as shown in Fig. 2(r).
O
Example 10. In this example, test the stability of the ~G- -PID as the dimensionality, d dy = dx = dy increases. By Property 1, if we take two i.idd systems of variables (M, X,XY, Y) at dimensionality d and concatenate their respective variables, every PID component of the composite system of dimensionality 2d should be double that of the original. This process can be repeated, taking two independent 2d-dimensional systems and concatenating them create 4d-dimensional system. Fig. 3 shows precisely this process starting with the system in Example with d = 2, and continually doubling its size until d 1024. The -PID accurately matches ground truth by doubling in value, and remains stable with small relative errors (shown in App. B.5). The PPD also runs much faster than the PID this example, providing speed-up of more than 1000 at d = 64 (see right-most plot in Fig. 3, and extended results in App. B.4).
Remark 2. In the above examples, the ground truth referred the PPIIID of Def. 1, since it was inferred using Def. 3 and additivity (Property 1). Since our method for computing the PPD recovers the ground truth PID, we can infer that the distribution of the optimal Q for the --PPDD is in fact Gaussian in these examples. This provides empirical evidence in support of the conjecture stated in Sec. 3.
Estimation and Bias-correction for the -PD
Having discussed how to compute the ~GPPID and shown that agrees well with ground truth in several canonical examples, we discuss how the PPD may be estimated from data. Given a sample
7

---

Blas correction the "Bit-of-alli" upt
Bias correction the "Fully-redundant" set up
2.00
Partial information value (bits)
1.75
Blased
Blas-corrected
Z7x
(X,;
information (bits)
B clarecoreected
Blased
2/7x
(X,)),
Ground truth
Miy
3.5
Ground truth
3.0
1.25
.57
2.5
1.00
1.5
0.50
1.0
0.25
10
0.00
0.0
100
200
300
600
1000
100 200
300
600
1000
Sample size
Sample size
Figure 4: Empirical evaluation of bias-corrected PID estimates with increasing sample size for two configurations with dır = dx = dy = 10, as described Section 5. Solid horizontal lines represent ground truth (as computed from the covariance matrix); dark colored bars represent biased PID components and light colored bars represent bias-corrected PID components, as estimated from the sample covariance matrix. Overlaid box plots indicate results from 100 random draws. Empirically, we find our bias-corrected estimates are both unbiased and consistent.
of n realizations of M, X' and Y drawn from PArXY, we may estimate the sample joint covariance matrix Easxy. We therefore use the straightforward "plug-in" estimator for the ~g-IID, by using ÊxxY in place of M in the optimization problem in equation (9).
However, it is well-known that estimators of cntommttion-theoreee quantities suffer from large biases for moderate sample sizes [37]. Cai et al. [38] characterized the bias in the entropy of d-dimensional Gaussian random vector, for fixed sample size n.
Proposition 3 (Bias in Gaussian entropy [38]). Suppose M R has auto-covariance matrix ΣM. The entropy of M is (M) 2 log det det(2neE when ΣM is known [24]. For the sample covariance matrix Ess, the bias is given by:
dar
Bias [MM) = Σ log(1 - k/n).
(15)
k=1
For a proof, we refer the reader to [38, Corollary 2]. This result may be naturally extended to compute the bias of each of the mutual information quantities in the LHS of equations (1)-(3):
Corollary 4 (Bias in Gaussian mutual information). For the joint mutual information I(M; (X, Y)),
dxs
dx+dy
dis + +dy
Bias [X, Y))
Σ log(1 - k/n) +
Σ log(1 - k/n) - Σ log(1 k/n) (16)
This follows directly from the fact that T(M: (X, Y))
Similarly, we can compute the bias of i(M; X) and Y(, But this does not uniquely determine the bias in the individual PID components, and as with defining PIDs, need to decompose the bias in Corollary 4 across the four PID components such that they agree with these constraints. We solve this problem by defining bias-corrected version of the union information from Proposition 1.
Definition 4 (Bras-corrected Union Information). We assign the bias in the union information to be the same fraction as the bias in the joint mutual information T(M; X,Y)). This gives rise correcorde estimate of the union information:
(M X;Y) T8(M (.MM::Y(())) Bias [M: (X,Y))]
(17)
bias-ComT
i(M; (X,
We do not analyze theoretically whether the bias-corrected union information is consistent and unbiased, thus, it may still have some residual bias relative to the true union information. However, we find empirically that this bias correction process works reasonably well and appears both consistent and unbiased, in number of examples. Figure 4 shows biased and bias-corrected PID values for

---

values (left) and values normalized by ((M:(X) (right) in multivariate Poisson spike-tount simulation
Normalized partial information
companent
Partial information (bits)
(X.)X))))
9.9
0.6
0.4
Myy
S7
9.2-
PID definition
3.2
>.PPPD
81
A--PID
MII
Ground truth
01 022 03 DA 05 d.6 0.7 0.5 0.
0.1 0.2
03 5-1
555
0.5 0.7 0.8 0.0
Polsson PPPDD
Weight from Mj to X
Weight from M4 to X
Figure 5: A comparison of the ~g-PID, the Pg-PID and the MMI-PID for a multivariate Poisson system. The ground truth is discrete PPDI computed using the package of Banerjee al. [30]. The ~G-PID comes closest to the ground truth (possibly because they compute the same PID definition). despite the fact that the ~GG only uses the covariance matrix of the Poisson distribution, whereas the ground truth uses knowledge of the distribution itself.
Redundancy betw 3 visual cortical areas (top 20 principal companentss)
Redundancy fraction betw visual areas (top principal carn ponents)
Redundancy, RE/VSpp VISI, VISaa) (blls)
Change
Ri(VISp: Viia Visally tyy (VIS) Visai)
Change Non-change
Non-change
14-
ns
f15
115
10-
4-
0.4
2
50
100
150
50 100
150
Time after sti muluus us onset (ms)
Time after sti mtimulus onset (ms)
Figure 6: Bias-corrected redundancy estimates for information about VISp activity that is shared between VISI and VISal: redundancy in bits (left) and redundancy as fraction of total mutual information (right). Data spread across 40 mice. Statistical comparisons use a two-sided Mann- Whitney-Wilcoxon test. Observe that there is greater and more sustained redundancy flashes corresponding to image changes, which are behaviorally relevant and linked to rewards in this task.
100 runs of two configurations called "Bit-of-all" (with a little bit of each PID component) and Feulda-redundatt" (which has predominantly redundancy), with dM dx = dy = 10 (details and additional setups in App. C). We find that bias correction brings the PID values closer to their true values even at small sample sizes. In App. C.3, also include preliminary analysis of the variance of PID estimates using bootstrap [39, Ch. 8].
6 Application to Simulated and Real Neural Data
So far, we have only considered the PPD when applied to Gaussian PMXY- Although Def. 2 strictly applies only to Gaussian PMXY, the estimation process in Sec. 5 relies only on sample covariance matrix, which is well-defined for a wide variety of non-Gaussian distributions. Many applications where PIDs could be useful have non-Gaussian data. For instance, there is great interest in applying PIDs in neuroscience (e.g., to understand how multiple brain regions jointly encode and communicate information [6, 40]), but spike-count distributions are non-Gaussi...
Simulated neural data. To show that our PID estimates provide reasonable results on non- Gaussian spiking neural data, we first simulate spike-count data using Poisson random variables (following [12]; described in App. D.1). We evaluate the ground truth ~PPI for this distribution using the discrete PID estimator of Banerjee et al. [20]. The ~ is estimated from sample covariance matrix using 10th realizations of M, X and Y. We find that the PPD closely matches

---

the ground truth for a range of parameter values, despite the fact that the PGD is effectively computed on a Gaussian approximation -a-PID applied to non- Non-Gaussian data are provided in App. D.2 (also see App. A.2). We conclude that it is reasonable to use and interpret the ~GDP on non-Gaussian spike count data.
Poisson distribution (Fig. 5). More examples of the
Real neural data. We then applied our bias-corrected ~gDPID estimator to the Visual Behavior Neuropixels dataset collected by us the Allen Institute [41]. We recorded over 80 mice using six neuropixels probes targeting various regions of visual cortex, while the mice were engaged in visual change-detection task. In the task, images from a set of 8 natural scenes were presented in 250 ms flashes, at intervals of 750 ms; the image would stay the same for a variable number of flashes after which it would change to a new image. The mouse had to lick to receive a water reward when the image changed. Thus, given image flash could be behaviorally relevant target if the previous image was different, or not, if the previous image was the same.
We used our bias-corrected PID estimator to understand how information is processed along the visual hierarchy during this task. We estimated the PPD to understand how information contained in the spiking activity of primary visual cortex (VISp) was represented in two higher-order visual cortical areas, VISI and VISal. We aligned trials to the onset of a stimulus flash, binned spikes in 50 ms intervals and considered the top 20 principal components (to achieve reasonable estimates these sample sizes; explained further in App. D.5) from each region in each time bin. We computed the PPDI the sample covariance matrix of these principal components (shown in Fig. 6). We found that there was a significantly larger amount of redundant information about VISp activity between VISI and VISal for stimulus flashes corresponding to image change, compared to flashes that were not changes (Fig. 6()). The larger redundancy was also sustained slightly longer for flashes corresponding changes, than non-change flashes. Both of these effects were maintained even when the redundancy was normalized by the joint mutual information, suggesting that the effect was not purely due to an increase in the total amount of information (Fig. 6(r)). Our results suggest that the visual cortex propagates information throughout the hierarchy more robustly when such information relevant for behavior.
Discussion
In this paper, we proposed a new and efficient method for estimating the PGD for Gaussian distributions. We showed that our method recovers the ground truth and suitably corrects for bias through a series of examples. In particular, Fig. 4 showed how large the biases can be at small sample sizes, which makes bias correction particularly important in neuroscientifie settings where sample sizes are often small.
We focused Gaussian distributions, as they have historically been starting point for many estima- tors (c.g., correlation is used as a measure of dependence, but zero correlation implics independence only in the Gaussian case). We were able to show ground-truth validation for our ~PPPID estimator at high dimensionalities only thanks to the existence of closed-form results on scalar Gaussians.
While central claims and results applied Gaussian distributions, our method for computing PIDs did not immediately break down for distributions close to Gaussian in some limit (e.g.. Poisson). The effective spiking rate used in our multivariate Poisson distribution was also not particularly large (see App. D.1); we would expect our estimates to improve for higher firing rates, since the Poisson distribution will then be more Gaussian.
Our work has several limitations that require further theory and simulations to resolve, the most important of which are: (1) Our estimator is technically bound the PID values because we assume Gaussian optimality in Definition 2; (2) Our bias-correction method is heuristic: we do not provide rigorous theoretical characterization of the bias of PID values.
Linitations.
Broader impacts. Our work is mainly methodological, so the scope for negative impacts depends on how the methods might be used. For example, incorrect interpretations drawn from the use of our PID estimators may affect scientific conclusions. In particular, the PID is inherently a correlational quantity and carries the same caveats: it may not be appropriate to make causal interpretations on the basis of observed PID values. Also, despite our best efforts to explore a variety of systems, we cannot tell how accurate our bias-correction method will be in novel configurations.
10

---

Acknowledgments and Disclosure of Funding
We thank Lukasz Kušmierz for providing a valuable reference on the bias of Shannon entropy estimates. We also thank Gabe Schamberg and Christof Koch for helpful discussions.
We wish to thank the Allen Institute founder, Paul G. Allen, for his vision, encouragement, and support. Praveen Venkatesh Interface of Data and Neuroscience, supported in part by the Allen Institute. Stefan Mihalas was in part supported by NSF 2223725, NIH R01EB029813, and RF1DA055669 grants.
supported by Shanahan Family Foundation Fellowship at the
References
[1] Saskia EJ de Vries, Jerome A Lecoq, Michael A Buice, Peter A Groblewski, Gabriel K Ocker, Michael Oliver, David Feng, Nicholas Cain, Peter Ledochowitsch, Daniel Millman, et al, large-scale standardized physiological survey reveals functional organization of the mouse visual cortex. Nature neuroscience, 23 (1):138-155 151, 2020.
[2] Joshua H Siegle, Xiaoxuan Jia, Séverine Durand, Sam Gale, Corbett Bennett, Nile Graddis, Greggory Heller, Tamina K Ramirez, Hannah Choi, Jennifer A Luviano, et al. Survey of spiking in the mouse visual system reveals functional hierarchy. Nature, 592(7852)886-92, 2021.
[3] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Ken- neth D Harris. Spontancous behaviors drive multidimensional, brainwide activity. Science, 364(6437); caav7893, 2019.
[4] Elad Schneidman, William Bialek, and Michael Berry. Synergy, redundancy, and independence population codes. Journal of Neuroscience, 23(37):11539-11$5 2003.
[5] Itay Gat and Naftali Tishby, Synergy and redundancy among brain cells of behaving monkeys. Advances in neural information processing systems, 11, 1998.
[6] Giuseppe Pica, Eugenio Piasini, Houman Safaai, Caroline Runyan, Christopher Harvey, Mathew Diamond, Christoph Kayser, Tommaso Fellin, and Stefano Panzeri. Quantifying how much sensory information in a neural code is relevant for behavior. Advances in Neural Information Processing Systems, 30, 2017.
[7] Giuseppe Pica, Mohammadreza Soltanipour, and Stefano Panzeri. Using intersection information map stimulus information transfer within neural networks. RioSystems, 185:104028, 2019.
[8] Jan Bim, Vito De Feo, Daniel Chicharro. Malte Bieler, Ileana Hlagaanu-Opate Andrea Brovelli, and Stefano Panzeri. A non-negative measure feature-related information transfer between neural signals. BioRxiv, page 758128, 28, 2019.
[9] Tomas Seagliarini. Luca Faes, Daniele Marinazzo, Sebastiano Stramaglia, and Rosario Mantegna. Synergistic information transfer in the global system of financial markets. Entropy, 22(9):1000, 2020.
[10] Patricia Wollstadt, Sebastian Schmitt, and Michael Wibral. A rigorous information -etororticcccoiicicccc definition of redundancy and relevancy in feature selection based on (partial) information decomposition. arXiv preprint 1rXiv:2105.04185 2021.
[11] Sanghamitra Dutta, Praveen Venkatesh, Piotr Mardziel, Anupam Datta, and Pulkit Grover. An information- theoretic quantification of discrinination with exempt features. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 3825-388 3833, 2020.
[12] Praveen Venkatesh and Gabriel Schamberg. Partial information decomposition via deficiency for multivari- ate Gaussians. In 2022 IEFE International Symposium on Information Theory (ISTT), pages 2892-2897. IEFF, 2022.
[13] Adam B Barrett. Exploration of synergistic and redundant information sharing in static and dynamical Gaussian systems. Physical Review E, 91(5):05280, 2015.
[14] Nigel Colenbier, Frederik Van de Steen, Lucina Q Uddin, Russell A Poldrack, Vince D Calhoun, and Daniele Marinazzo. Disambiguating the role of blood flow and global signal with partial information decomposition. Neuroimage, 213:116699, 2020.
[15] Tjeerd W Boonstra, Luca Facs, Jennifer N Kerkman, and Daniele Marinazzo. Information decomposition of multichannel EMG map functional interactions the distributed motor system. Neurolmage, 202: 11609, 2019.
[16] Jana Krohova, Luca Faes, Barbora Czippelova, Zuzana Turianikova, Nikoleta Mazgutova, Riccardo Pernice, Alessandro Busacca, Daniele Marinazzo, Sebastiano Stramaglia, and Michal Javorka. Multiscale informa- tion decomposition dissects control mechanisms of heart rate variability at rest and during physiological stress. Entropy, 21(5):526, 2019.

---

[17] Ari Pakman, Amin Nejatbakhsh, Dar Gilboa, Abdullah Makkeh, Luca Mazzucato, Michael Wibral, and Elad Schneidman. Estimating the unique information of continuous variables. Advances neural information processing systems, 34:20295-20307, 2021.
[18] Nils Bertschinger, Johannes Rauh, Eckehard Olbrich, Jürgen Jost, and Nihat Ay. Quantifying unique information. Entropy, 16(4):2161 2183, 2014.
[19] Johannes Rauh, Pradeep Kr Banerjee, Eckehard Olbrich, Guido Montáfar, and Jürgen Jost. Continuity and additivity properties of information decompositions. arXiv preprint arXiv:2204. 10982, 2022.
[20] Pradeep Kr Banerjee, Johannes Rauh, and Guido Montúfar. Computing the unique information. In 2018 IEEE International Symposium on Information Theory (ISIT), pages 141-145, IEEE, 2018
[21] Abdullah Makkeh, Dirk Oliver Theis, and Raul Vicente. BROIA-2PID. A robust estimator for bivariate partial information decomposition. Entropy, 20(4):271, 2018.
[22] Paul Pu Liang. Yun Cheng. Xiang Fan, Chun Kai Ling, Suzanne Nic, Richard Chen, Zihao Deng. Faisal Mahmood, Ruslan Salakhutdinov, and Louis-Philippe Morency. Quantifying modcling feature interactions: An information decomposition framework. arXiv preprint arXiv:2302. 12247, 2023.
[23] Paul L Williams and Randall Beer. Nonnegative decomposition multivariate information. arXiv preprint arxi 1004. 2010.
[24] Thomas M Cover and Joy A Thomas. Elements of Information Theory. John Wiley & Sons, 2012.
[25] Virgil Griffith and Christof Koch. Quantifying synergistic mutual information. In Guided self-organization: inception, pages 159-190. Springer, 2014.
[26] Malte Harder, Christoph Salge, and Daniel Polani. Bivariate measure of redundant information. Physical Review E, 87(1):01210,, 2013.
[27] Artemy Kolchinsky. A novel approach to the partial information decomposition. Entropy, 24(3):403, 2022.
[28] Joseph T Lizier, Nils Bertschinger, Jürgen Jost, and Michael Wibral. Information decomposition of target effects from multi-source interactions: perspectives previous, current and future work. Entropy, 20(4): 307, 2018.
[29] Praveen Venkatesh, Keerthana Gurushankar, and Gabriel Schamberg. Capturing and interpreting unique information. arXiv preprint arXiv:2302. 1873, 2023.
[30] Pradeep Kr Banerjee, Eckehard Olbrich, Jürgen Jost, and Johannes Rauh. Unique informations and deficiencies. In 2018 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton). pages 32-38. JEEE, 2018.
[31] Naftali Tishby, Fernando C Pereira, and William Bialek. The information bottleneck method. arXiv preprint phyxics/0004057, 2000.
[32] Gal Chechik, Amir Globerson, Naftali Tishby, and Yair Weiss. Information bottleneck for Gaussian variables. Advances in Neural Information Processing Systems, 16. 2003.
[33] Amir Globerson and Naftali Tishby. On the optimality of the Gaussian information bottleneck curve. The Hebrew University of Jerusalem, Tech. Rep, page 22, 2004.
[34] Martin Riedmiller and Heinrich Braun. A direct adaptive method for faster backpropagation learning: The RPROP algorithm. In IEEE International Conference on Neural Networks, pages 586-591. IEEE, 1993.
[35] Geoffrey Hinton. Coursera neural networks for machine learning, lecture 6. 2018. URL https:// www. cs. toronto. edw/eeduuaaaaa "tij csc3cc http///ddddseeeerreeeeee satrre__ssliie slides_lect .ddf. Also see https:// optimizas on. cbe. cornell edu/ index. wwtp:e/ppppptttteeespoooo.
[36] Code this paper. URL https: http://gtthhab..ooaaaeee//////// cow.ppravenn253//g2 gpid.
[37] Liam Paninski. Estimation of entropy and mutual information. Neural computation, 15(6): 1990:119911100 2003.
[38] T Tony Cai. Tengyuan Liang, and Harrison H Zhou. Law of log determinant of sample covariance matrix and optimal estimation of differential entropy for high-dimensional Gaussian distributions. Journal of Multivariate Analysis, 137:161-172, 2015.
[39] Larry Wasserman. All of statistics: a concise course statistical inference, volume 26. Springer. 2004.
[40] Nicholas M Timme and Christopher Lapish. A tutorial for information theory in neuroscience. eneuro, 5
(3), 2018.
[41] Allen Institute. Visual behavior neuropixels dataset overview, 2022. URL htttttss://portae http:///ootaallreee brain map. org/expii http///poooee/criit////////////e www.aaao--eaaar--eeeeeeeeeeeeeee neuropixxx
[42] Volker Strassen. Gaussian elimination is not optimal. Numeerisce mathematik, 13(4):354-355 1969.
12

---

Gaussian Partial Information Decomposition: Bias Correction and Application to High-dimensional Data
Appendices
Praveen Venkatesh, Corbett Bennett, Sam Gale, Tamina K. Ramirez, Greggory Heller, Séverine Durand, Shawn Olsen, Stefan Mihalas
A Supplementary Material for Sections 2 and 3
A.1 Implications of the Additivity Property
As mentioned in the main text, additivity states that the PID values of an isolated system should not depend on the PID values of another isolated system. Without additivity, it not possible to examine independent systems in isolation, since broadening your view to include a different isolated system could change the PID values of the first system.
Hypothetically, suppose we have two separate individuals (labeled 1 and 2) receiving completely independent stimuli, and we examine the PID between the activity in brain regions M, X and Y in each of their brains. Then the total unique information that X1 and X2 have about M1 and M2 with respect to Y and Y2 should be equal to the sum of the unique information in each of the individuals taken separately.
As another example, suppose we are trying to examine visual information flow and auditory informa- tion flow in multi-sensory integration task. We may want to understand the degree to which the activity in the sub-regions of the visual and auditory systems depend on each other. A reasonable null model of independence between the two systems would be that the PID values of the joint system will be equal to the sum of the PID values in the two individual systems. Then, measuring the actual degree to which the joint PID) value is not equal to the will be meaningful measure of dependence between the systems. This would only be possible with a PID definition that guaranteed additivity of independent sub-systems. However, since the 5-PID does not satisfy the additivity (Property 1), we cannot guarantee that the aforementioned null would be the correct null model, and we would not be able to perform such an analysis.
A.2 Explaining the Exponential Reduction in the Number of Optimization Variables
In Def. 2, we restrict QMXY to be jointly Gaussian. This reduces the number of optimization variables from being exponential in the dimensionality to quadratic, as we show here.
Previous discrete PPDI estimators [20, 21J have found efficient methods to solve the optimization problem of Def. without reducing the number of optimization variables. If PMXX is discrete in Def. 1, and each dimension of M, X and Y has support of size K (i.e., M,, X, and Y, can each take one of K discrete values), then the total number of degrees of freedom in Ap is (KK(dart (xx+))) This corresponds to the total support of a discrete QMXY, which can be an arbitrarily complex discrete distribution.
contrast, the PPD in Def. reduces the dimensionality of the optimization space by assuming that QMXY is Gaussian. QMXY is then completely parameterized by its covariance matrix, NQ M In fact, given the additional constraints in Def. 2, the only part of 2& M xy that is variable is Xx.v. as explained in Sec. 3.2. This matrix has only dx x dy variables. Thus, by restricting the optimization space to jointly Gaussian distributions, we have reduced the number of variables from exponential in the dimensionality to quadratic.
As concrete example, consider a neuroscientific setting (as in Sec. where M, X and represent the spiking activity of neurons in three different brain regions, so that d dx and dy represent the number of neurons, and K represents the maximum number of spikes. Makkeh et al. [21] demonstrated their method for a maximum support size of K = 18 with du dx = dy = 1, for total support of approximately 183 = 5831 in QMXY- we assume a single neuron can have at most K = 10 spikes, then get effective dimensionality of (5831) < 4. In other words.
13

---

these discrete r--PID estimators can handle around 4 neurons across all three regions M, X and Y (which is also what we use the simulation of Fig. 5 in Sec. 6), whereas our PG- estimator can handle several hundreds if not thousands of neurons (as seen in Example 10).
A.3 Proofs of Propositions and
Proof of Proposition 1. Firstly, the differential entropy of a Gaussian random variable M with covariance matrix Σn is given by [43, Thm. 8.4.1]:
((M) log (2Te
(18)
Secondly, for a joint Gaussian distribution PMY parameterized by covariance matrix MMXY, the conditional covariance matrix can be written as [44, Sec. 8.1.3]:
ΣXrIM Exy E EN E XXYM
(19)
Exy = E "xyas + E
XXYM D I "XY
(20)
Using these equations, derive the mutual information between M and (X,Y) follows:
I(M;(X,Y)) = h(X,Y) - h(X, M)
(21)
1
log det (et((x)) log det(22e2ixx M)
(22)
1
log ((2Te) (Exx)) log (2T((2 ane del (El((i)(MM)
(23)
2
1
det(Exx))
log
(24)
2
det((Ex M)
det(2ixVIM + Ixy, 2
(
log
(25)
2
det(ExyIM)
dde((Ex M) det(! Exy EXY,M-XX,
En' IT 2xy,))
1
log
(26)
2
det(Exy]M)
1
log det(/ + E-1 X XXYM My MM M )
(27)
2
( log det(1 + M YY レー X X 2xx,,r))
(28)
where in (a) we used equation (18), (b) we used equation (20), and in (c) we used the fact that det(/ + AB) det(/ BA).
The remainder of the proof follows from the arguments presented below equation (8). The constraint in Proposition arises because, when optimizing over Σx.y|a we require ΣMpx to be a valid positive semidefinite covariance matrix, i.e., ΣMXY > 0. This happens if and only if Eas and its Schur complement in Σarx are both positive semidefinite, i.e., Eas > and En - Exy TT EXYJM > O.
M,
Proof of Proposition The proof divided into three parts consisting of derivations for the objective, the gradient and the projection operator.
Objective. After whitening the Px m and the PYJM channels, and assuming that EM Remark 1), without loss of generality we have that
(see
Σx.M = E (H M + Mx) = Hx HMM = H
(29)
I
ΣX.yIM 1
ΣXYIM =
(30)
TT
ΣX I
Hx
+
Σn -T XXYMM NYMM HJ H..
ET
xx
Hy
(31)
14

---

For the sake of brevity, let Σ represent the optimization variable Ex,y,. and let S be its Schur complement in Exy M. 1 E.T Then, the inverse of 2xylas can be written as [44, Sec. 9.1.5]:
8-1
-5-1E
E E XXIM
(32)
-
1+
Therefore, we get:
EN It XXYM = Hi HP 1
-T - I+ T+E
8
-
][
Hx Hy
(33)
HH + -ily) (Hx-) -EHy)
(34)
Thus, setting B = Hx EHy, the optimization problem in Proposition reduces to
min
log det (1 +H + BBB
(35)
s.t. EXY M 0
Gradient. Let the objective derived in the previous section be called f(Z), where Σ == Σx,y Ex.yIM as before. We can compute the gradient of f with respect to ) using standard identities from matrix calculus. First, note that the gradient of scalar function with respect matrix is itself matrix with entries as follows:
0f
V((2)
(E). EEij
(36)
Considering each element of this matrix:
0f
a
(Σ)
log det(I + Hy + BBS-B))
(37)
2 20
Σ
(a) Tr (T H Hy + S-
(1 H Hy + BJ s-i
(38)
1 2 In (I+ Hy B
(B B)
(39)
Σ
where in (a), we have used the identity d = aa()) [44, Sec. 2], while in (b), we use the fact that only and S depend on Σ implicitly, with the other terms being constants.
Expanding the partial derivative alone, we get:
a S B + BT. a - B + BTS O
(B& S S-BBBB
(B
(B)
(40)
FEg
5
EEij
FEg
FEi
x
wherein
0
(B)
= وΣ
(Hx ΣHy)
(41)
5
_ Hy,
(42)
as
(5-1)
S-1 HEi 5-1
(43)
S-1
(1 ΣΣ 1
(44)
- --
(45)
where JM is the single-entry matrix, containing I at location (i.j) and 0's everywhere else; in (b) and (d), use the fact that aX/ax.j Ju [44, Sec. 9.7.6]; and in (c) we use the identity
15

---

0(x-)= = =x- - 0 (BTS-BB) 12ny
Therefore, (40) becomes
(46)
-- + (J + EJ + (-J
(47)
- + B'S-TJEE'S + B - J
(48)
Putting all together, and letting A = (1 H Ily + BB B). (39) becomes
ff
(x)
TT H B +
B
+ =
(49)
- ",,"
2
+ "=AA""55" TrAA-1' ||
(50)
(c) TKS--BA-- + THA"S-BBA-115--11
2 + TASS Tr (51)
where in (e), we have used the fact that the trace of a matrix product is invariant under cyclic permutations of the matrices within the product.
Finally, using the fact that [ زツح Tr(WJ Wi for any matrix W [44, Sec. 9.7.5],
05 1
[ 2 - + BA-
(52)
ØEg
BA BA (B' 11}}
(53)
= Vf(E) = 5 BA (B -BB - H
(54)
= =S- +
H4 H + BB ( -H
(55)
Projection operator. Recall that the optimization variable, ا = Ex,yJas is an off-diagonal block of E'y(M Mr which the matrix upon which the constraint is defined:
Σ
ΣXYIM = ΣT T -
(56)
wherein the diagonal blocks are identity due Remark 1. For the purposes of this section, let us suppose EXyIM is function of Σ, ΣxyIM =: g(E), so that the constraint may be written as g(2) > 0. A suitable projection operator, therefore, will accept a value Eg (that may violate (())) 0) and find a point projj close to it that satisfies the constraint, i.c., g(2 b 0.
We do not find the "orthogonal" projection operator, which has the minimum distance ||_ Eoll in some norm. Instead, we propose a simple heuristic to find Droj which satisfies the constraint.
Eg satisfies the constraint, then we are done, so let us assume that gEE)) * O. Then, we can find a matrix ΣXYJM which close to (Eo)) and satisñiss ΣXyJM >= 0 as follows: let the eigenvalue decomposition of (Eo) be given by VAVT, with A =: (Aag being the diagonal matrix consisting of its eigenvalues X; Then, since g (Eo) is not positive semidefinite, s.t. X; < 0. We set X, == 0 for all such i; effectively, X, = max 0, } V We then reconstruct the matrix using these "rectified" eigenvalues and set it to be ΣXYJM == VÄVT, where = diag(Ax).
Now, we need to find Σ such that g(E) = ExymM- However, ΣxrΙM may not have identity matrices on its diagonal blocks, i.e., might not correspond to a whitened channel. We therefore whiten
16

---

Exy M as follows:
1-1 X|M 0
Σ
Σ-
0
Iihhitned XYYM
ExYIM
XIM
(57)
E-1/2 SYMM
Σ 1/2 -YM
where Σx M and SyJm are the diagonal blocks of EXYJM. Crucially, since the matrix multiplying Exy M on either side is itself (the inverse square-root of) a covariance matrix (and hence positive semidefinite), Enhhionndd is also positive semidefinite.
Now, the off-diagonal block of Whittened will satisfy g(-) = "X IM
Iiwhieened
"XYMM b 0. This off-diagonal block
forms the output of our projection operation and can be written as
Eproj X.Y|M
- > Y
- XIM
(58)
which comes directly from equation (57).
A.4 Details of PPD Optimization and RProp Implementation
The optimization problem for the ~c-PID, using projected gradient descent with RProp (mentioned in Section 3), implemented as follows:
1. Let Σx,yJa be shorthand for the optimization variable, and let Proj(-) represent the projection operator defined in Prop. 2. Let ( represent the value of Σ at iteration i of the optimization. Initialize ( = Proj(Hy, H where Hý is the pseudoinverse Hy.
2. Evaluate the objective and the gradient as defined in Prop. 2, the current value of Σ() Compute the sign of (each element of) the gradient,
(( = Sgn Sgn(V,(E(()),
(59)
When computing the objective and the gradient, add small regularization term to the computa-
tion of s-1 (as defined in Prop. 2): S = ((1 + t) -E where we take = 10-7.
3. Update:
( = Proj - a²1 o ()))))
(60)
where Proj is the projection operator defined in Prop. 2; () is time-varying learning rate vector of the same dimension as Σ, describing the learning rate for each element of Σ; o represents an clement-wise (or Hadamard) product between vectors; and = 0.999 is constant, which when raised to the power of imposes a slow overall decay of the learning rate to promote convergence.
The matrix inverses embedded in the projection operator are also regularized by modifying Equation (58) as follows:
sproj X M
(7I+ E1/2 - -E "ExyI E
(61)
where y is slowly increased from 10-12 12 by factor of 10 in each step, until g(_x,) > 0.
4. (0)) is initialized to 10-3 and (7()) is updated as follows:
(++1) r ( @ B W(E ()OW)))))
(62)
where 3 = 0.9 is a constant that determines how fast the learning rate increases or decreases; and all operations are carried out element-wise Note that when some element of the gradient changes in sign, that element of (+++)) @ (E()) will be positive, resulting in decrease in that element of ()) On the other hand, if the sign of some element of the gradient remains the same, then the learning rate for that component will increase by a factor of 1/0.9.
5. Stop when the absolute differences between the current objective and the previous objectives from the last 20 consecutive iterations are all less than 10-6 ("patience") )) or when the maximum number of iterations is exceeded (set to 10 iterations).
17

---

M 10, A.5 Computational Complexity of the -GP
Suppose for simplicity that M, X and Y all have the same dimensionality d= dM dx = dy. Then, Mx and all its submatrices have side of dimensionality O(d). Thus, the computational complexity of each gradient descent iteration determined by the complexity of the matrix operations in equations (10)-113) These operations include matrix multiplication, matrix inverses, log dct(-) and eigenvalue decomposition, which have a worst case complexity of (dd))
Matrix multiplication can potentially be performed at lower complexity for large matrix sizes (e.g.. see Strassen [42]), which has downstream implications for each of the other operations as well. However, for our purposes, we take the computational complexity of the objective, the gradient and the projection operator to all (dd))
B Supplementary Material for Section 4
First, observe that by subtracting equation (2) from equation (1), we have
I(M: (X, Y)) I(M;X) = Uly + SI
(63)
I(M;Y = Uly SI.
Similarly, subtracting equations (1) and (3), we get that T(M; X Y) UIx + SI. These two equations hold in general, and will be used in what follows.
B.1 Details and Derivations for Examples in Section 4
Example 2 (Pure uniqueness).
(64)
= M Nx
Nx, Ny ~ i..dd (N, NN) 1 M
(65)
Y = NY
(66)
Derivation of PID values in Example 2.
Y M
7
I(M;Y) =0
(67)
Uly + RT= =0
(68)
Uly RI 20
Uly RI =
(69)
UIx I(M;X)
(70)
ST
X(
(71)
I(M:X) XM); =0.
(72)
Example 3 (Pure redundancy).
M ~ N(O,1) 1) X M Nx
(73)
Nx N((01)
(74)
Y M + Nx
Nx 1 M
(75)
Derivation of PID values in Example 3.
[[M; X IY)
Ulx + =1 =
(76)
I(M;Y X) =
Uly SI SI =
(77)
RI = I(M: X(M;X).
(78)
Example 4 (Pure synergy).
M N(0,111) 1)
(79)
= = Nx
Nx N()
(80)
18

---

M2
12
1/2
1/2
X2
X2
Y1
M1
0
m1
0
X1
V1
X1
Figure 7: Diagrams explaining Examples and 9. See App. B.2 for details.
N=
Nx T M
(81)
Derivation of PID values in Example 4.
YLM
I(M;Y) =0
(82)
Uly RI
(83)
Uly, RI 20
Uly = RI = =0
(84)
Ulx ((M;))
(85)
SI I(M; (X, (MM;))
(86)
00 ((M;)) = 0c.
(87)
should be noted that certain nuances have been omitted in discussing Examples 2-4 above. For instance, in Example 3, ΣxyJM rank deficient and hence invertible,, which would be an issue when computing the objective in Equation (9). Also, in Example T(M; co, however this could be corrected by adding some noise to either X or so that their difference is noisy representation of M.
Example 5 (Unique and redundant information).
M N(O, 1)
(88)
X M Nx
Nx ~ N(O.1) N2 ~N(,00
Nx M NÝ L (Nx,M)
(89)
Y = M + N NN
(90)
Derivation of PID values in Example 5. Essentially, X is a noisy representation of M, while Y is noisy representation of X. Since M-X-Y forms Markov chain, I(M;Y Y X) = and hence Uly = SI O. When 2² = 0, this example reduces to Example 3 with only redundancy being present. For any finite non-zero value of o², both RI and U Ix are present and are non-zero. Since M is scalar, the redundancy for the P-D is identical to the MMI-PID's redundancy [13]:
RI min (MM; (M; I(M;Y)
(9))
since I(M;Y) X), by the data processing inequality. At the limit when o² I(M;Y) + 0, and therefore RI + 0. while U UIx will become equal to X).
oc.
Example 6 (Unique and synergistic information).
M A 1)
(92)
X M Nx
Nx, Ny Ny ~ ~((,00² (on(Nx,Ny)
(Nx, M
(93)
Y Ny
(94)
Derivation of PID values in Example 6. When p = and 0² + oc, this example reduces to Exam-
ple 4, with only synergy being present. In general, Y 1. M, therefore, I(M;Y) UTy + RT= = 0. meaning UTy RT = 0. For any finite value of o², X will have some unique information about M given by Ulx X(M;X) 0. Correspondingly, SI I(M;X Y) Ulx
19

---

~ vs. 6-PID at dyr
dx dy
Partial information (bits)
PID
60
20
0.00
0.33
0.67
0.99
1.67
2.00
2.33
2.67
3,00
in X,
Figure 8:
A comparison of the ~G- and PG-PIDs at dimensionality of d 64. The plot is virtually identical to that in Fig. 2 (which uses d = 2), except that the y-axis is scaled up by a factor of 32.
[[M;X - X). When = 0, [M;XX = X) and therefore SI = 0. As 1, Y M: so the total mutual information I [MM; (X, Y)) - co, driven by synergy growing unbounded, while the unique component remains finite at X).
Example 7 (Redundant and synergistic information).
M N(00 1)
(95)
X M + Nx
Nv, 10,,1), Ny)
(Nx, 1 M
(96)
Y M Ny
(97)
Derivation of PID values in Example 7. When p = 1, we once again reduce to Example 3 with only redundancy. When 1, we cannot infer the PID values using Equation (1) and non-negativity alone, since none of the individual mutual information values (or conditional mutual information values) go zero. Instead,
determine the redundancy using the MMI-PID since M is scalar. Note that I (M; and [(M; Y') are both equal by symmetry, and thus equal to RI. This also implies that both UIx and U ly must be equal zero. As p reduces, the two channels X and Y have noisy representations of M with increasingly independent noise terms. Therefore, their average, (X Y)/2 will be more informative about M than either of them individually, meaning that X and Y jointly contain more information than any one individually. This extra information about M is synergistic, given by SI [MM; X Y), and increases as p decreases, attaining its maximum possible value at 0.
B.2 Diagrams Explaining Examples 8 and 9
Examples and 9 can be understood diagrammatically as shown in Fig. 7()) and Fig. 7(r), respectively. In both diagrams, we represent the two-dimensional plane describing M, with axes mi and m2- The colored vectors shown on this plane represent Hx and Hy, i.e., the gain with which X and Y represent each value of M. For example, Y2 captures only M2. with gain corresponding to its length. The gains are directly representative of the signal-to-noise ratio (and hence the amount of information) in each variable, since the noise in each variable is i..dd, with unit variance. In Example 8, the gain in X₃ variable, while in Example 9, the angle at which X₃ and X₂ sample M1 and M2 is variable.
B.3 Comparison of the ~G- and PPDS in Example Higher Dimensionality
For completeness, we also present comparison between the ~a- and SG-PIDs for Example 8 at higher dimensionality of d = 64. The results can be found in Fig. 8, and show that both PIDs scale proportionally. We use gain of a 0.99 in place of 1, since gain of unity causes certain matrices to become ill-conditioned.
Although the PPD does not obey the additivity property in general (which is why it does not agree with the ~G-PID at gains of a > 1), it appears to double along with the ~p-PID (i.e., additivity across identical independent copies appears to hold, at least in this specific example).
20

---

Timing analysis the ~ --DD for the doubling example
101
10'
10-1
16
128 256
512
1024
Dimension of M, X
Figure 9: A plot showing the time taken compute the ~GGPID at different dimensionalities in Example 10. This plot extends what is shown in the right-most plot of Fig. 3 for the ~g-PID.
Absolute error the cDPID over increasing dimensionality
Abs. error in PID value (bits)
[[M; (X,Y))
Ulx
Uly
RJ
SI
10
10
10
10
10
16
64
420
258
512
1024
32
64
128
256
512
1024
1024
18
256
1022
32
68
120
256
512
Dimension of each of M, and y
Relative error in the ~ PPDIID over increasing dimensionality
(X.:XX)))
UIx
Uly
RI
S1
Rel. error in PID value
10
10:42
10442
32
64
128
256
512
1024
16
32
64
128
256
512
1024
16
32
512
32
64
128
256
512
1024
16
32
64
128
256
512
1024
Dimension of each of M, and Y
Figure 10: Absolute (top) and relative (bottom) errors in computed PID values from Example 10.
B.4 Run Time Analysis of the PPD in Example 10
In Fig. 9, we present complete picture of the time taken to compute the PG- function of dimensionality d, in Example 10. Note that we cannot compare with the Sc-PID beyond a dimensionality of d 64, because the computation for the PPD failed in setup for d 64. Neither the P nor the ~GDPID were carefully profiled and optimized, therefore, the runtime analysis in both Fig. and Fig. 3 should be considered preliminary. We leave a more thorough comparison of the run time both these methods to future work.
B.5 Absolute and Relative Errors in Example 10
Figure 10 shows how the absolute and relative errors in PID values scale with increasing dimensional- ity in Example 10. The absolute emrors increase in proportion to dimensionality, starting under 10-7
21

---

at = 2 and remaining under 10-4 at d = 1024. The relative errors are all roughly constant, and remain under 10
Supplementary Material for Section
C.1 Implementation Details for Bias-correction
We use a number of different setups based on sampling from random connectivity matrices for bias correction in Section 5, All of these setups assume that dx = dy.
The both-unique, Kully-redundan and high-synergy setups have the following in common:
MM = I
(98)
ΣxIM T
(99)
Σy|M
(100)
I
H Hx EXIM Hy + Ew
H Hy + Ew Hy Hy ΣyJM
ΣMXY
Hx Hx
(101)
Also, the elements of Hx are either zero or one, H Hi ~ ii.d. Ber(O.1) These three setups differ in their definitions of Hy (the channel gain from M Y) and Ew (which controls the extent of correlation between X and Y).
The both-unique setup draws Hy (i, ii.d Ber(O.1), with all elements of Hy independent the elements of Hx, and sets Σw = 0.
The fully redundant setup similar Example by setting Hy Hx and Ew 0.91 (note that Ew is square, since dx dy). By keeping Ww close to the identity matrix, we are effectively in the regime with high correlation p Example 7. This allows to come close emulating Example without suffering from the issue of non-invertibility of ExyJM, mentioned in App. B.
The high-synergy setup is similar to Example 6, by setting Hy 0 and Ew 0.8I. As with the relldated setup, by keeping Σw close to the identity matrix, we are in the high-p regime. This allows us come close emulating Example 4, while not making the synergy or the total mutual information infinite.
The zero-synergy setup is similar to Example and uses the following setup:
m = ΣxIM =
(102)
Σx Hx H + ExIM
(103)
Σy x =
(104)
I
H
HY Ex H + Eyix
ΣMXY
Hx
Σx
(105)
Hx
Here, (i, f) ~ i.i.d. Ber(0.1), while Hy = Hx, with (i,j) ~ i.i.d. Ber(O.1), 1), Hy 1 Hx- so that ST
Defined this way, M-X-y form Markov chain, ensuring that [M: (Refer equation (63)).
The bit-of-all setup is combination of equal parts of the high-synergy and zero-synergy setups The variables X and Y are swapped in the zero-synergy setup, so that both X and Y can have some unique information.
Remark (Rectification). In practice, we observed that the bias correction procedure prescribed in Definition 4 could lead to negative values for certain PID quantities. This occurred because the bias-corrected union information was not guaranteed to satisfy certain bounds, which we enforce below. To prevent the occurrence negative PID values after bias-correction, we require a form of rectification:
= max 1
i(M; X) bias-corr
i(M; Y)
(106)
bias-comt
bias-corr
22

---

(1) reat
i(M; X)] bias-corr
((; bias-corr
7(M;(X, Y)) bias-cort
(107)
where bias-cor tion equation above, the union information is bounded from below by the individual (bias-Corrected) mutual information values, and bounded from above by the sum of the individual mutual information values, and by the total mutual information.
biissooo represents a bias-correcter mutual information estimate. After the second rectifica-
C.2 Bias-correction Performance in Additional Setups and at Higher Dimensionality
Plots showing bias correction performance for all setups described in App. C.1 shown in Fig. for 10-dimensional, and in Fig. 12 for 20-dimensional M, X and Y.
Of the setups examine, only the case with both X and Y having purely unique information appears to have somewhat poor performance, where our bias-correction method appears to over-correct the bias in unique information, while insufficiently correcting the bias in redundancy and synergy.
C.3 Preliminary Analysis of the Variance of PID Estimates
In Figures 13 and 14, we present preliminary analysis of the variance of our PID estimates using bootstrap. The figures represent the true distribution of the PID estimates over multiple sample draws, over multiple bootstrap sample draws, in the form of box plots. In what follows, refer to these box plots as "confidence intervals". The true "confidence intervals" were estimated using 100 runs of bias-corrected PID estimates, i.e., by drawing 100 different samples, each of size The bootstrap "confidence intervals" were estimated using 100 bootstrap samples that were resampled from a single randomly drawn sample of size n.
colloquially
When correcting for bias in the PID estimates on bootstrap samples, we use the number of unique data points in each bootstrap sample in place of n (refer Corollary 4), rather than the total sample size. This leads to more stable bootstrap-PID estimates.
The quality of the bootstrap "confidence interval" is affected greatly by the quality of the individual sample used for bootstrap resampling. Nevertheless, we observe a reasonable degree of qualitative agreement between the true "confidence interval" and the bootstrap "confidence interval", particularly as the sample size increases. Future work will assess confidence intervals with greater care, using well-defined metries, and assess how well these confidence intervals are calibrated.
D
Supplementary Material for Section
D.1 Details Regarding the Multivariate Poisson Spike-count Simulation
We follow our previous paper [12], where this analysis first presented. In this simulation, M is two-dimensional, consisting of two independent and identically distributed Poisson random variables, Mj and M2. X and Y are each generated through a linear combination of Binomially thinning Mi and M2, along with some Poisson noise:
Mi, M2 Poiss(2)
(108)
X ~ Binom [M, a) + Binom
(nmom(Me (M2, 0.5) + Poiss(1)
(109)
Y Binom M1, 0.5) Binom(Ma 0.5) + Poiss(1)
(110)
While our results in Fig. show that our ~PDD estimator comes closest to obtaining the same values the estimator of Banerjee et al. [20], several remarks are warranted:
1. The method of Banerjee et al. [20] estimates the ~PPID and not the 8-PID. Since there is no discrete estimator for the 8-PID, it is entirely possible that the 8G-PID is equally (or more) accurate with respect to its own true value. That is, the difference between the P and the Banerjee al. [20] estimator is entirely due difference in definitions;
2. We stated in the introduction that the PID estimates approximate upper bound, whereas our PPD estimates an exact upper bound. The difference due to the approximate nature of the PPD upper bound is also not captured here, nor in Examples 8 and 9. In all these examples, we cannot tell how much of the difference is due the difference in definitions and the approximate nature of the estimators.
23

---

Bias correction in the "both-unique" set up (d= 10)
Bias correction the "lulll-rreuundan" set up (d 10)
Portial information value (bits)
Blased
(XM; X))
Information [bits]
Blased
(X, Y)
Bias-corrected
Uxy
Blas-corrected
2/y
Ground truth
tiyy
3.5
Ground truth
3.0
ST
2.0
1.5
1.0
0.0
100
200
300
600
1000
100
200
300
600
1000
Sample size
Sample size
Blas correction the "high-syneryy" set up (d 10)
Bias correction in the "zero-synergy" set up (d= 10)
Partial information value (bits)
Blased
I(M, (X,)))
4.0
Blased
(X,))
Blas correctdd
77x
8 Blas cas-ccorectedd
2/7
Ground truth
Uiy
Ground truth
E/y
3.0
.ST
information
2.5
1.0
0.5
44
100
200
300
600
1000
100
200
1000
Sample size
Sample size
Bias carrection the "bit-of-all" up (d=10)
3.0
Blased
(M,
Partial information value (bits)
Bias-correcte
DTy
2.5
SI
1.5
1.0-
100
200
300
600
1000
Sample size
Figure 11: Bias correction for various setups described in App. C. with d = dM dx = dy 10.
24

---

Bias correction in the "both-unique" set up (d 20)
Bias correction in the "fully-redundant set (d 20)
Portial information value (bits)
16
Blased
(X,, K))
Information (bits)
17.5
X[M (X,Y))
Bias-corrected
Uxy
Bias corrected
Ground truth
tiyy
Ground truth
U7y
12-
10
10.0
8-
7.5
6-
5.0
2.5
0.0
100
200
300
600
1000
100
200
300
600
1000
Sample size
Sample size
Blas correction the "nigh-synergy" set up (d 20)
Blas correction the "zero-synergy" set up (d=20)
Partial information value (bits)
Blased
I(M; (X,)))
Biased
(X),
Blas-correcter
77x
17.5
Bias- correctdd
2/x
20
Ground truth
Uiy
Ground truth
E/y
15.0
.ST
15
information
12.5
10.0
10-
5.0-
2.5
100
200
300
600
1000
100
200
300
600
1000
Sample size
Sample size
correction "bit-of-all" set up 20)
17.5
Blased
[MM; (X,)3)
Partial information value (bits)
15.0
Bias-corrected
DTy
12.5
SI
10.0
7.5
5.0
100
200
600
1000
Sample size
Figure 12: Bias correction for various setups described in App. C. with d = dM dx = dy = 20.
25

---

Bootstrap the "both-unique" setup (d= 10)
Bootstrap Cls "fully-reduuddan"" setup (d= 10)
True CI
(X.X))
True CI
(X.
3.5
PID value (bits) Cls
Bootstrap CI
UTy
Bootstrap CI
truth
Uiyy
3.0
Ground truth
R7
RI
SZ
2.5
(bits)
2-
1.5
1.0
19
100
200
300
600
1000
100
200
1000
Sample size
Sample size
Bootstrap Cis "high-synergy" setup
Bootstrap "zero-synergy" setup (d
74
True
(X,
True CI
(X,))
Bootstrap
Bootstrap
PID value (bits) with Cis
Ground truth
Uly
Ground truth
Uly
RZ
(bits) with
RI
4-
3-
21
200
1000
100
Sample size
Sample size
Bootstrap Cls the "bit-of-all" (d=
True CI
(X,
Bootstrap
D7x
PID value (bits) with Cls
Ground truth
CTy
6-
S1
4-
100
300
1000
Sample size
Figure 13: Bootstrap "confidence intervals" for various setups described in App. C.I with d == dxr = dx dy 10. Note that these are not true confidence intervals, but box-plot representations of the true variance of the estimator (over 100 runs) and the bootstrap estimate of the estimate's 's variance (from 100 bootstrap resamplings of single random sample).
26

---

Bootstrap Cls "both-unique" setup
Bootstrap "fully-redundanl" setup (d-
25
True CI
(X,)))
True CI
(X.
17.5
Bootstrap
UTy
Bootstrap
PID value (bits) with Cis
15.0
Ground
Miyy
Uly
R7
RI
12.5
S7
15-
value (bits)
10.0
5.0
5-
2.5
100
200
300
600
1000
100
200
300
600
1000
Sample size
Sample size
Bootstrap the "high-synergy" setup (d 20)
Bootstrap "zero-synergy" setup
25
True CI
(X,
20.0
True
(X,
PID value (bits) with Cis
Bootstrap
17.5
Bootstrap
Ground
Uly
Ground truth
Uly
RZ
15.0
RI
15-
ST
value [bits)
12.5
SI
10.0
10
7.5
5.0
2.5
0.0
1000
300
1000
Sample size
Sample size
Bootstrap the "bit-of-all" setup
True Ci
(X,.XX)))
PID value (bits) with Cls
14
Ground truth
Bootstrap
DTx
C7T,
12
.S1
10
6-
100
1000
Sample size
Figure 14: Bootstrap "confidence intervals" for various setups described in App. C.I with d == dx = dx dy 20. Note that these are not true confidence intervals, but box-plot representations of the true variance of the estimator (over 100 runs) and the bootstrap estimate of the estimate's 's variance (from 100 bootstrap resamplings of single random sample).
27

---

3. Finally, we do not know how accurate the "ground truth" itself is, since this is also assessed using the estimator of Banerjee et al. [20]. Differences between the various estimators could also be due to a poor estimate produced by the discrete PID estimator.
D.2 Additional Simulations on Non-Gaussian Data
In order to evaluate how the ~- performs on a greater variety of non-Gaussian distributions, we considered three more setups similar to the one described in App. D.I.
The first additional setup is a multivariate Binomial distribution:
M1, M2 Binom (4, 0.5)
(111)
X ~ Binom (,, ox) + Binom M2 0.5) + Binom 0.5)
(112)
Y Binom( 0.5) Binom M2, 0.5) + Binomm2, 0.5)
(113)
This setup has the same mean as the multivariate Poisson distribution, but has a smaller variance and is therefore more peaked, further from zero. The multivariate Binomial distribution is closer to Gaussian than the multivariate Poisson, and indeed, we see good agreement between the ~G-PlI and the "ground truth", assessed by the estimator of Banerjee et al. [20] (see Fig. 15, top panels).
The next two setups are "zero-inflated" versions of the multivariate Poisson and Binomial distributions, which are produced by passing each of the aforementioned variables through a "Z-channel", with zero-out probability 0.3. That is,
M M1 ZM,
X' = X 7x Y' Y x Zy
(114)
M2 M2 x ZM,
(115)
where Z, i.d. (0r(.. i E M,, The primed versions of these variables have an additional probability mass at zero, with a probability of 0.3. This is said to better reflect the statistics of Calcium imaging [45], although we use a different base distribution than the one in the cited paper. The pero-inflated versions of the Poisson and Binomial distributions one of the peaks at zero, and the Binomial version having more well separated peaks. For both zero-inflated versions, the PPD performs poorly at recovering the absolute P-DI values (Fig. 15, middle-left and bottom-left panels). However, the relative PID values (normalized by the total mutual information) are closer their true values (Fig. 15, middle-right and bottom-right panels), providing hope that a better mutual information estimate (e.g., [46, 47J) could correct the absolute PID) values to some extent.
bimodal, with
Further validation on non-Gaussian data is necessary to understand the impact of non-Gaussianin on GG-II estimates. However, a more extensive evaluation of non-Gaussian distributions (especially at higher dimensionalittes) is made difficult by the unavailability of ground truth. The "ground truth" in these examples obtained using the discrete PID estimator of Banerjee al. [20] (which we assume is more accurate for discrete variables). Our method provides a better result than the 8-PID and the MMI-PID in Fig. 5, least part of which is probably due the difference in definitions. Nevertheless, we believe these examples demonstrate that applying our method to non- non-Gaussian data need not be ruled out, although care should always be taken interpreting results.
D.3 Implementation Details of the Analysis Pipeline
The Visual Behavior Neuropixels data was analyzed follows:
1. We selected mice that had at least 20 units in each brain region of interest. Only mice with both familiar and novel sessions were selected.
2. From each region we selected units of 'good' quality, with SNR at least 1, and with fewer than 1 inter-spike interval violations.
3. Trials were aligned the start of each stimulus flash, and spikes were counted in bins of 50 ms, between 0 and 250 ms after stimulus onset (0-50 ms, 50-100 ms, etc.).
4. Trials corresponding to non-change flash were defined as those that occurred between 4 and 10 flashes after the start of a behavioral trial, such that the image remained the same as the original image in this behavioral trial. Flashes corresponding to an omission, flashes after an omission, and flashes during which the animal licked, were all removed. Only flashes that occurred while
28

---

PID values (left) and normalized IM; (J, (right) in multivariate Binomial spike-count simulation
0.7
partial information
0.6
PID component
(MM(X.
Partial Information [bits)
0.5
Uix
Rt
0.4
0.2
Normalized
definition
0.1
P:D
Ground truth
(Bancrjee al. ...PID]
0.0
0.0
0.2
02
0.3 0.4
0.5
0.6 0.7
0.8
0.9
0.1
0.2 0.3 04 0.5 0.6 0.7 0.8 0.9
Weight from
Weight from to
PID values (left) and PID normalized by IM; Y)) (right) zero-inflated Poisson spike-count simulation
0.25
Normallized lized partial information
PID companent
AMF;(X,
0.20
Partial information
UTy
RI
0.15
0.30
definition
0.06
0.1
Ground truth
(Banerjee al. -PDD
0.00
0.0
01
0.2
0.3
0.4
0,5
0.6
0.7
0.8 0.9
0.1
022 0.1 044 0.5 0.6 0.7 0.8 0.9
Weight from
X
Weight from M, to
PID values (left) PID normalized by IM; (X,Y)) (right)
zero-inflated Binomial spike-count simulation
0.25
Normalized partial information
0.6
PID component T.M(XX
0.20
0.5
U1x
Partial information
0.15
definition PID
0.00
Ground truth (Banerjee ot al. -.D)
0.00
0.0
01 0.3 0.4 0.5 0.6 0.7 0.8 0.9
01
0.3 0.5
0.7 0.9
Weight from M;
X
Weight from M; to
Figure 15: PID values and normalized PID values for three different non-Gaussian distributions, as described App. D.2. The top, middle and bottom panels are respectively for a multivariate Binomial, zero-inflated Poisson, and zero-inflated Binomial distribution. The left column represents the PID values bits, while the right column represents the PID values normalized by the total mutual information. All values plotted function of the gain from M1 to X (a) the z-axis. An explanation the results can found in the text.
29

---

Redundancy betw 3 visual cortical areas
Redundancy fraction betw visual cortical areas (top 10 principal components)
(top 10 principal components)
RI(VISp: ISp. VISI; VISal (bits)
12-
#+++
Change
RIVISS VISI, (1 (VIS), VISwJ)
Change
Non-change
Non-change
10
n5
2043005
20805
n5
0.8
115
Redundancy RIIV
ns5
0.6
4000
0.4
2-
0.2
100
150
200
50
100
200
Time after stimulus onset (ms)
Time after stimulus onset (ms)
Figure 16: Redundancy about VISp activity between VISI and VISal, function of time, for flashes corresponding to an image change (blue) and flashes corresponding to a non-change (orange). Data points are across 41 mice. The plot on the left shows the redundancy in bits, while the plot on the right shows the redundancy normalized by the total mutual information.
Redundancy betw 3 visual cortical areas (top 20 principal components)
Redundancy fraction betw visual cortical areas
(top 20 principal components)
RI(VISp: VISI; VISam) (BSaa)/I[(V) VISp (VISI, VISwm))
Redundancy, RM(VISp: VISI. VISam) [bits]
Change
Change
Non-change
Non-change
ns
12
ns
ns
10
0.8
ns
0.6
0.4
2-
0.2
100
150
200
100
150
200
Time after stimulus onset (ms)
Time after stimulus onset (ms)
Figure 17: Redundancy about VISp activity between VISI and VISam, function of time, for flashes corresponding to an image change (blue) and flashes corresponding to a non-change (orange).
Data points across 38 mice. The plot the left shows the redundancy in bits, while the plot on the right shows the redundancy normalized by the total mutual information.
30

---

PID values in bits; M VI8p, VISS, Y- ViSal (top PCS)
time
time 50
time 100
time 150
time - 200
PID value (bits)
Condition
Charge Nonichange
Auor 01x
Diy
61
Juxr 68x
Jasr Dig
Juce D0g
D3y 11
Juce Mix
PID components over time (in ms)
PID fraction of total mutual information: VISp, VISI, Y= Vi8al (top 20 PCs)
150
PID traction
0.4
avauurr
Charge Non-shange
0.2
biy
ole
dl
dr
d1
S1
PID components over time ms)
Figure 18: All PID components VISp activity, between VISI and VISal-in bits (top), and fraction of total mutual information (bottom), various times after stimulus onset, for change and non-change flashes. Here, the label L in the x-axis of the top plot refers to I [M; (X, Y)). These plots show that redundancy is the primary driver of mutual information.
PID values bits; M VISp, X VIS, V- VISam: (top 20 PCs)
time
time 50
time 100
time 150
time 200
PID value (bits)
10-
Condition
Charge change
Aucy Dix D/y
Tux
G3x
biy
Juse D1x
Juco D1x Dby RI
Fucr Dix
PID components over time (in ms)
PID fraction of total mutual information; M VISp, X VISI, VISam PCS)
time 150
0.8
PID fraction
0.6
eee Wet +
Condition
Charge Non-change
0.0
67
67
0/4
65,
21
PID components over time ms)
Figure 19: All PID components about VISp activity, between VISI and VISam in bits (top), and fraction of total mutual information (bottom), various times after stimulus onset, for change and non-change flashes. Here, the label 4Mx in the x-axis of the top plot refers to I (M; (X, Y)). These plots show that redundancy is the primary driver of mutual information.
31

---

the animal was engaged (as measured by an average reward rate of at least 2 rewards/min) were selected.
5. Trials corresponding to change flash were defined those during which an image change occurred, and the animal was engaged (as above).
6. The top 10 or 20 principal components of neural activity selected at each time bin, and for each brain region under consideration. Principal component analysis was carried out using the Scikit-learn [48] package in Python.
7. ~GPPID estimates were computed on the covariance matrix between principal components across regions, for each time bin.
8. Data were aggregated across 40 or 41 mice for the figures with VISal, and across 38 mice for the figures with VISam.
9. Statistical significance was assessed using two-sided unpaired Mann-Whitney-Willee test.
D.4
Additional Results
Fig. 16 shows the redundancy about VISp activity, between VISI and VISal. is identical to Fig. 6, except that it uses only the top-10 PCA components in each of the three visual cortical regions. This figure shows that the result presented in Fig. gracefully degrades with a reduction in the number of PCA components. When using fewer PCA components, the increase in redundancy on change trials is still statistically significant, but the effect size is not as large, and not as sustained, as it is in Fig. 6.
Fig. 17 shows the redundancy about VISp activity between VISI and VISam, the latter of which is different higher-oud order cortical region (see, e.g., [49]). This figure shows that the result in Fig. 6 can also be seen for another choice of higher-order cortical region (i.e., VISam rather than VISal).
Figures 18 and 19 show all PID components, not just the redundancy, for the same settings as in Figures 6 and 17 respectively. These show that redundancy the dominant partial information component, and appears to be the main driver of changes in the overall mutual information. This justifies why we include only plot of redundancy Figs. 6, 16 and 17.
D.5 The Necessity for PCA, and Maximizing Dimensionality
In analysis of real neural data in Sec. we used principal components analysis reduce dimensionality, despite the fact that our method could scale to much larger dimensionalities. As mentioned briefly Sec. 6, this was due the limited number of trials we had. For many mice, there were a larger number of neurons in each region than trials with which to compute the covariance matrix. In other words, directly computing the covariance of the neural activity would have resulted in covariance matrix that was aank-deficeent. We used PCA to reduce dimensionality and ensure that we obtained a reasonable estimate of the covariance matrix, and to minimize the error in our PID estimate, which would naturally be higher with greater number of PCA dimensions.
For completeness, we analyze here what would happen we used the maximum possible number of PCA components for each mouse, while ensuring that the covariance matrix is not rank deficient. We also ensure that had the same number of components across all regions and across both change and non-change conditions (within each mouse), so to perform fair comparison. This gives us on average 53+t PCA components, with a minimum of 22 and a maximum of 84 PCA components. We find that the results of higher and more sustained redundancy hold when using the maximum possible number of PCA components for each mouse (see Fig. 20). However, as expected, there is much greater variance across mice, possibly due to variability in the number of PCA components chosen across mice, and due to greater errors in our PID estimates. For the main results presented in Fig. 6, in order to be consistent across mice, we rounded down from the minimum and used the top 20 PCA components for all mice.
change flashes continues to
D.6 Differences between Change and Non-change Conditions are not an Artifact of Bias-correction
The number of trials corresponding to change flashes is much smaller than the number of trials corresponding to non-change flashes. Accordingly, the sample size used to estimate the covariance matrix is different in each of the two conditions. Bias correction performed using the appropriate
32

---

Redundancy betw 3 visual cortical areas
Redundancy fraction betw visual cortical areas
(max principal components)
(max principal components)
Redundancy, RR VISI, VISal) (bits)
Change
RI(VISp: VISI; VISI; VISal) /IWVISS;; (PIIVVSS VISal)
Change
20
Non-change
Non-change
ns
1.2
15-
1.0
0.8
10-
0.6
0.4
0.2
100 150
200
50 100
200
Time after stimulus onset (ms)
Time after stimulus onset (ms)
Figure 20: Redundancy about VISp activity between VISI and VISal, using the maximum possible number PCA components for each mouse. Data points are across 42 mice, however, not all mice are represented at all points in time and for both conditions, since individual data points may not have converged due ill-conditioned matrices. The plot on the left shows the raw redundancy in bits, while the plot on the right shows the redundancy normalized by the total mutual information. The increase in redundancy on change trials is still statistically significant, however, the errorbars have greater spread for reasons described in the text.
Redundancy betw 3 visual cortical areas (top 20 principal components)
Redundancy fraction betw visual cortical areas (top 20 principal components)
Redundancy, RI(IISp VISI, VISal) (bits)
16
Change
RI(VISp: VISI, V Visal) / /XIVISp: (VISI, VISal))
Change Non-change
Non-change
ns
ns
12
1.0
*K
n5
0.8
ns
200040
0.6
0.4
2-
0.2
50
100
150
200
50 100
150 200
Time after stimulus onset (ms)
Time after stimulus onset (ms)
Figure 21: Redundancy about VISp activity between VISI and VISal, as a function time, for flashes corresponding to image change (blue) and flashes corresponding to non-change (orange) with an equal number samples. Data points are across 41 mice. The plot on the left shows the redundancy in bits, while the plot on the right shows the redundancy normalized by the total mutual information. The observations made in the other figures continue to hold; the differences seen between the two conditions are. therefore, not result of differences in sample size.
33

---

sample size; however, as noted in Section 5, our bias correction process is not perfect, and may leave some residual bias.
In order to show that the results we observed were not an artifact of differences in residual bias caused by different sample sizes, we randomly subsampled the non-change flashes to produce dataset with equal numbers of trials for change and non-change flashes. Repeating the analysis as before, we found that our conclusions continued to hold even in the setting where both conditions have equal sample sizes, shown in Figure 21.
D.7 Note on Inhomogeneous Poisson Processes
Inhomogeneous Poisson processes are those whose mean firing- (or "emission-") rates change with time. Our method for estimating the PID does not consider the temporal characteristics of the original signal. Rather, the data analyst can choose the random variables M, X and Y to span some time range (or different time ranges) as they please.
In our analysis of real neural data in Sec. 6, we counted the number of spikes in a 50-125 ms window after stimulus onset. Even if the underlying spiking process was an inhomogeneous Poisson process, this spike count would be Poisson distributed, with a mean given by the integral of the emission rate over the fixed window. In general, while analyzing data and computing PID values, one will have to be aware of severe inhomogeneities (e.g.. if the distribution changes not just in the same way in each trial, but changes across trials), and account for them separately. For example, in our analysis, we excluded time periods when the mice were not "engaged"" in the task, as defined by not actively consuming rewards at a rate of at least rewards per minute.
E Compute Configuration Used and Code Availability
All analyses were performed on a workstation equipped with an Intel Core 17-10700KF CPU with 8 cores (16 threads), 48 GiB of RAM and data stored on TB PCIe NVMe solid state drive.
Analysis of the Visual Behavior Neuropixels data (for 84 sessions) with d = 20 took approximately 9 minutes to run. This included loading data for each session and computing 840 PID values on 60 x 60 covariance matrices, implying an average run-time of about 0.64s for each PID estimate (including amortized data-load time).
All code used to compute and estimate the PPD and correct for bias, including all examples in this paper and code for neural data analysis, is available GitHub [36].
Additional References
[43] Thomas M Cover and Joy A Thomas. Elements of Information Theory. John Wiley & Sons, 2012.
[44] K. B. Petersen and M. S. Pedersen. The Matrix Cookbook. Technical University of Denmark, 2012. URL http://www2.. compute. dtu. ck/ /akk///p///pppaaapeeaaeapbbaap// htcuhb/uuss//aa /2774 full .html. Version: November 15, 2012.
[45] Xue-Xin Wei, Ding Zhou, Andres Grosmark, Zaki Ajabi, Fraser Sparks, Pengcheng Zhou, Mark Brandon, Attila Losonczy, and Liam Paninski. A zero-inflated gamma model for post-deconvolved calcium imaging traces. Neurons, Behavior, Data Analysis, and Theory, 3(2):1 21, 2020.
[46] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger. Estimating mutual information. Physical review E, 696661661,3,, 2004.
[47] Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm. Mutual information neural estimation. In International conference on machine learning, pages 531 5400 PMLR, 2018.
[48] F. Pedregosa, G. Varoquaux, Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, Dubourg. J. Vanderplas, A. Passos, Cournapeau, M. Brucher, M. Perrot, and E, Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
[49] Joshua H Siegle, Xiaoxuan Jia, Séverine Durand, Sam Gale, Corbett Bennett, Nile Graddis, Greggory Heller, Tamina K Ramirez, Hannah Choi, Jennifer A Luviano, et al. Survey of spiking in the mouse visual system reveals functional hicrarchy. Nature, 592(7852) 2021.
34


Reviews:
[{"review_id": "8HaUqzWYHV", "reviewer": "Reviewer_jMyH", "summary": "First, I would like to thank the authors for the work they put into contributing to the field. \nIn their paper, the authors provide a new method to estimating a well-known measure, the Partial Information Decomposition (PID). Within the topic of PIDs, their method computes a specific version of PID, the ~-PID. \n\nWhile methods to compute the ~-PID already exist, the authors' method is advantageous in that \n1) it is more efficient to compute (quadratic in number of variables instead of exponential) and \n2) in that it corrects for the bias. \n\nThe author's then go on to show with simulations and examples that \n1) their choice of PID method, the ~-PID, outperforms other PID methods\n2) the bias correction works empirically\n3) their method is stable even for increasing dimensionality\n4) their method can be applied on data which is not Gaussian distributed such as simulated Poisson data and real neural recordings (neuropixel)", "questions": "- In the list of contributions, the authors state in line 50 that they are able to reduce the number of optimization variables from exponential to quadratic. Where can I see this in the equations and derivations? \n\n- The authors show that their method works on non-Gaussian data by applying it to simulated Poisson data and neuropixel spike count data. They explain this by stating that the Gaussian distributions is a good-enough approximation to the Poisson distribution. It is known, however, that neural firing does not follow a Poisson distribution. This is especially pronounced when recorded via 2photon imaging for which the distribution is strongly inflated at zero. Xue-Xin Wei et al. 2020 propose a zero-inflated gamma model to accurately capture calcium imaging traces. Is it feasible for the authors to test their method on such bimodal distributions which might be harder to approximate with a Gaussian distribution? If the rebuttal time for such experiments is too short, can the authors comment on the expected outcome?", "limitations": "The authors discuss and point out most limitations. The claim that their method works on non-Gaussian data is only backed by an experiment with Poisson data, however. In my opinion, there is still room for the method to fail for other relevant distributions (refer to point \"Questions\") which would be worth mentioning as a (possible) limitation.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- The paper addresses a topic which is becoming more and more relevant as the number of neurons that can be recorded from simultaneously is increasing: The question of how information is distributed among neurons and brain areas.\n- The paper is written in a very clear and structured way which successfully guides the reader through every logical step. \n- The mathematical formulation is sound\n- The quality of figures is excellent", "weaknesses": "- The paper is an improvement of an already existing method which makes it important and publishable yet disqualifies it for a very high grade.\n- The claim of generalizability to non-Gaussian data is not backed up strongly enough empirically or mathematically (refer to \"Questions\")\n- Small mistakes/typos in lines: 177, 206, 322", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n> The paper is an improvement of an already existing method which makes it important and publishable yet disqualifies it for a very high grade.\n\nWe thank the reviewer for raising this concern, and providing us with an opportunity to explain why we believe our advance is significant:\n\n1. In the overall rebuttal, we explain how our advance is significant with respect to other works in the field, and in particular the work of Venkatesh et al. (2022). In particular, our method works at higher dimensions, runs much faster and corrects for bias.\n\n2. Additivity is an extremely fundamental property: effectively, it states that the PID values of an isolated system should not depend on the PID values of another isolated system (please see our response to reviewer DDWz for an explanation). The ~-PID satisfies the additivity, whereas the delta-PID does not. This makes our method much more practically applicable than the method based on the delta-PID. Apart from additivity, several works have shown that the ~-PID also has other strong operational motivations (e.g., Kolchinsky, Entropy 2022), which has not been shown to the same degree for the delta-PID. These points (additivity in particular), make it extremely important to have a method capable of computing the ~-PID. We will add an explanation of the importance of additivity to the revised version of our paper.\n\n3. Bias correction is extremely important for obtaining correct PID values, as evidenced by the sheer amount of bias seen at small sample sizes (Fig. 4, and Figs. 9 and 10 in the appendix). For example, we find that synergy is often highly over-estimated. To our knowledge, ours is the first paper to demonstrate the severity of not correcting for bias when estimating PIDs from real data, and we are also the first to propose and evaluate a method to correct this bias. We will re-emphasize the importance of bias correction in the revised paper.\n\n> The claim of generalizability to non-Gaussian data is not backed up strongly enough empirically or mathematically (refer to \"Questions\")\n\nThank you for raising this point. We did not mean to claim generalizability to non-Gaussian distributions in general. Rather, since our Gaussian PID method uses only a covariance matrix, we tried to justify its application to non-Gaussian spiking neural data before demonstrating its utility in a practical neuroscientific setting. Accordingly, we presented a small multivariate Poisson simulation (which was the largest possible example in which we could compare with known estimators), which we have now extended to include a few more non-Gaussian cases. We hope our new simulations highlight the limitations of the applicability of our method to non-Gaussian data, and show that care needs to be exercised when drawing interpretations from the results of our method. A more detailed response is given below.\n\n> Small typos ...\n\nThank you, we will correct these.\n\n**Questions:**\n\n> In the list of contributions, the authors state in line 50 that they are able to reduce the number of optimization variables from exponential to quadratic. Where can I see this in the equations and derivations?\n\nWe thank the reviewer for raising this question, and we will add the following explanation to the paper:\nThe general form of the ~-PID involves optimizing over Q_MXY. If we consider M, X and Y to be discrete, with a support of size K in each dimension, then the number of degrees of freedom of Q will be O(K^(d_M + d_X + d_Y)). For general continuous distributions, if we discretize the distribution, the same analysis would hold. In the case of our Gaussian ~-PID method, the optimization variable is \\Sigma_{X,Y|M}, which has a dimensionality of just d_X * d_Y.\n\n> The authors ... applying it to simulated Poisson data ... It is known, however, that neural firing does not follow a Poisson distribution. This is especially pronounced when recorded via 2photon imaging for which the distribution is strongly inflated at zero. ... Is it feasible for the authors to test their method on such bimodal distributions ...?\n\nWe agree that further validation on non-Gaussian data to understand the impact of non-Gaussianity is important. As described in the overall rebuttal, we have now included an example on a multivariate Binomial distribution (which is also close to Gaussian), as well as zero-inflated versions of both the multivariate Poisson and multivariate Binomial distributions. Our observations are outlined in the overall rebuttal.\n\nWe note that a more extensive evaluation of non-Gaussian distributions (especially at higher dimensionalities) is made difficult by the unavailability of ground truth. The “ground truth” in our multivariate Poisson example is obtained using the discrete ~-PID estimator of Banerjee et al. (which we assume is more accurate for discrete variables). Our method provides a better result than the delta-PID and the MMI-PID, at least part of which could be attributed to the difference in definition. Nevertheless, we believe it demonstrates that applying our method to non-Gaussian data need not be ruled out, although care should always be taken in interpreting results.\n\nWe will add both the new simulations, as well as this discussion on careful interpretation to the revised version of our paper.\n\nWith the addition of new results showing the extent to which our method applies to non-Gaussian data, we request the reviewer to reconsider whether our paper meets the criteria for a higher score."}, {"review_id": "Rqd07FxrMV", "reviewer": "Reviewer_Tp5P", "summary": "This article proposes an upper bound on mutual information (more precisely: on the \"unique information\" / \"union information\" that appear in \"Partial Information Decomposition\") that is easier to compute. This is done by replacing an infimum over all distributions matching given marginals by an infimum over only the ones that are Gaussian.\nToy examples (with Gaussian distributions) are given, and an experiment on real data is performed.\nFurther theory also includes considerations about biases arising from covariance estimators.", "questions": "I would be willing to revise my score if the points in the section above were significantly addressed.\n\nEDIT: they were.", "limitations": "Cf the validation issue in the weakness section above.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is mostly clearly written, with simple examples to help the reader follow.\nThe paper is self-complete, with reminders of definitions.\nConsidering covariance estimator biases is a plus.", "weaknesses": "EDIT: after discussion with the authors, I see that I had misunderstood the scope and contributions of the paper. The weaknesses below do not longer hold (or not as strongly).\n\n-----------\n\nThe main issue is that the definition of the upper bound actually supposes that the distribution MXY that is studied is Gaussian, at least marginally. Indeed the infimum is performed over all possible distributions Q_MXY that are Gaussian and whose marginals satisfy Q_MX = P_MX and Q_MY = P_MY, which implies that P_MX and P_MY have to be Gaussian (otherwise the infimum is performed over an empty set).\nAs a consequence, this upper bound definition cannot be applied to datasets that are not Gaussian (at least marginally). This crucial point is not discussed in the paper.\n\nActually, from Sections 3.2 and 5, it seems that the upper bound is estimated based on covariance matrices only, without using distribution Gaussianity. If this is the case indeed, then the definition could be changed, to rely directly on these covariance matrices and removing the Gaussianity assumption. This however would probably require to significantly rewrite the paper.\n\n\nAnother issue is the lack of sufficient validation, either theoretically or experimentally, in particular regarding distributions that are not Gaussian (i.e., really not Gaussian). This comment could have been tempered if there had been a discussion about general approximative Gaussianity of distributions in the field of study, but there is none. \nAlso, the impact of the introduction of such a Gaussianity constraint in the upper bound estimation should be studied closely: how realistic is this assumption, how far are real optimal Q_MXY from their Gaussian versions, how tight is the upper bound thus obtained, etc.\n\nIf this was a theoretical paper, I would expect theoretical results (guarantees such as a bound on the error made by the upper bound, or proving the conjectures, etc.). If this was an experimental paper, I would expect an extensive validation on many distributions for which the actual mutual informations are estimated (brute force or known solutions) in settings as varied as possible (not just normal distributions + a single real dataset). But the paper, in its current state, does meet these expectations.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n> The main issue is that the definition of the upper bound actually supposes that the distribution MXY that is studied is Gaussian, at least marginally … This crucial point is not discussed in the paper.\n\nIt is true that our definition applies only to Gaussian P_MXY (in fact, it must be jointly Gaussian, not just marginally), and that the definition technically does not apply to non-Gaussian distributions. This is why we include the word “Gaussian” in the title of the paper.  We mention “we compute/estimate the PID on Gaussian distributions”, in the abstract (lines 9-11) and in the introduction (lines 48-49). We also make it clear that M, X and Y are jointly Gaussian just above Definition 2 (lines 112-113), and again in lines 132-133. Overall, we believe we have tried to convey the point that P_MXY is assumed to be jointly Gaussian.\n\nNevertheless, we understand that there is a potential for confusion between the Gaussianity of P_MXY and the assumption restricting the optimization variable, Q_MXY, to be Gaussian. To make the Gaussianity of P_MXY more explicit, we will state this assumption within Definition 2 itself, rephrase line 120, and we will also explicitly make the point that the definition does not apply to non-Gaussian distributions, strictly speaking. We will also make small edits as necessary throughout the whole paper to ensure that this distinction is clear. We welcome any suggestions that may help make this point more clearly.\n\nEven though our method is designed only for Gaussian distributions, we believe it is an important contribution. Gaussian distributions have historically been a starting point for many estimators (e.g., correlation is used as a measure of dependence, but is exact only in the Gaussian case). Our ~-PID estimator is essentially the best in the field in terms of providing ground-truth validation at high dimensionalities; this was made possible because of Gaussianity. While the central claims and results of the study focus on Gaussian distributions, we showed that when distributions are close to Gaussian (e.g., Poisson) our method for computing PIDs does not immediately break down.\n\n> Actually, from Sections 3.2 and 5, it seems that the upper bound is estimated based on covariance matrices only … definition could be changed … would probably require to significantly rewrite the paper.\n\nAs the reviewer correctly notes, our definition depends only on the covariance matrix, and can thus be applied to any data, even if the data itself is non-Gaussian. We do not believe this requires a rewrite of the paper; rather, we just apply an estimator that was designed for Gaussian distributions to data that is non-Gaussian. This is what we do in section 6: while Gaussian P_MXY forms the main scope of the paper in Sections 1-5 (please see the overall rebuttal), we justify the applicability of our method to spiking data using a Poisson simulation.\n\n> Another issue is the lack of sufficient validation, either theoretically or experimentally, in particular regarding distributions that are not Gaussian …\n\nGiven that the main scope of the paper in Sections 1-5 is Gaussian distributions, we believe we have demonstrated sufficient empirical validation on Gaussian P_MXY. This has been shown through Examples 5-10 and Figures 1-4. Compared to other papers, we believe our paper provides one of the most comprehensive ground-truth validations (please see the overall rebuttal for a detailed comparison with [17] and [21]).\n\nHowever, we agree that it is important to justify that our method can be applied to non-Gaussian neuroscientific data, which we do in Section 6 and in the additional analyses (please see overall rebuttal). However, we do not make claims in general non-Gaussian cases, which could be very far from Gaussian. We will make this explicit in the revised paper. While a broader and more rigorous analysis of non-Gaussian distributions is important, we believe it is beyond the scope of the current paper.\n\n> … If this was an experimental paper, I would expect an extensive validation on many distributions for which the actual mutual informations are estimated (brute force or known solutions) …\n\nThis is an empirically-focused paper, and we do have extensive validation for Gaussian distributions, which is the main scope of the paper.\n\nAs the reviewer states, comparisons on a wider variety of Gaussian and non-Gaussian examples are limited by the availability of ground truth. This is why we rely on the MMI-PID (Definition 3) and Barrett’s theorem [13] (lines 177-178), which provides a closed form expression for the ~-PID for Gaussian P_MXY, _for scalar M_. Using this, we show that we can recover the ground truth by restricting Q_MXY to be Gaussian in Examples 5-7 (Fig 1). We then construct more complex examples using additivity (Property 1), and show that even then, the restriction to Gaussian Q_MXY can recover the ground truth (Examples 8-10; Figs 2, 3). To go beyond these examples, we do not know of brute-force methods for computing the ~-PID for Gaussian P_MXY, except for [17], which applies to general distributions P_MXY with scalar M, X and Y (which is trivial for Gaussian P_MXY due to Barrett’s result), and [21], which is not sufficiently well-tested to consider as “ground truth”, and was also published too recently for us to perform comparisons. For the same reasons, comparisons with non-Gaussian distributions at higher dimensions (beyond small Poisson/Binomial) are also limited by the availability of ground truth.\n\n> I would be willing to revise my score...\n\nTo summarize:\n- We will make the distinction between the Gaussianity of P and Q clear\n- For the scope of Gaussians, we believe our evaluations are comprehensive\n- We have performed additional analyses on non-Gaussian distributions, and we will be explicit that our method does not extend in general\n\nPlease let us know if our arguments and updated results are convincing, and how we may improve our score."}, {"review_id": "tfHnhUAEJv", "reviewer": "Reviewer_wULj", "summary": "The authors propose a new, efficient method for computing Partial Information Decompositions (PIDs) on multivariate Gaussians. They build their approach around the $\\sim$-PID approach, as this allows them to preserve an additivity property (allowing PIDs to be computed on independent systems separately and then added later). They present a number of canonical examples for Gaussians, shaw that their Gaussian PID works even when distributions are non-Gaussian, and show an example of its use on real neural data from the Allen Institute. Finally, they address the issue of bias in PID estimates, propose a bias-correction method, and evaluate it empirically.", "questions": "Can the authors method run on thousands of neurons, consistent with the original motivation they provide? \nWhat are the limitations on dimensionality more explicitly for this method? \nWhat is the computational complexity?\nDoes the author's method require using PCA or having high firing rate neurons? \nHow does it deal with low firing rate neurons or highly inhomogenous Poisson variables?", "limitations": "The authors adequately address limitations, although it would be good to hear more related to the questions above.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is very clearly written and the topics well explained. The literature review seems comprehensive, and although of somewhat specific topical interest, the work seems original and useful.", "weaknesses": "Given one of the major motivating factors given in the introduction is the need for efficient estimators of PID that can accommodate higher dimensional neural data, I was disappointed that in the end the authors chose to apply their method to a PCA-reduced version of the Allen Institute neuropixel probe data... wasn't the point to be able to address higher dimensional problems? It's thus not completely clear to me that this method, as is, has delivered on the promise of providing a PID approach more applicable to data with thousands of neurons than other PID methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses**\n\n> Given one of the major motivating factors given in the introduction is the need for efficient estimators of PID that can accommodate higher dimensional neural data, I was disappointed that in the end the authors chose to apply their method to a PCA-reduced version of the Allen Institute neuropixel probe data... wasn't the point to be able to address higher dimensional problems? It's thus not completely clear to me that this method, as is, has delivered on the promise of providing a PID approach more applicable to data with thousands of neurons than other PID methods.\n\nWe thank the reviewer for raising this concern and providing us with an opportunity to address it. We have now performed the PID analysis on the Allen Institute data again, using a larger number of PCA components (as described below; Figs. 8,9 in attached PDF). We were forced to use PCA because there were often a larger number of neurons in each region, than the number of trials available for computing the covariance matrix (in other words, if we had directly computed the covariance on the neural activity, the covariance matrix would have been rank-deficient). We wanted to ensure that we obtained a reasonable estimate of the covariance, and minimize the error in our PID estimate, which would naturally be higher for higher PCA dimensions.\n\nIn the new analysis, for each mouse in the dataset, we used as many PCA components as possible, based on the number of neurons in each region and the number of trials available for computing the covariance matrix. We used the maximum number of PCA components subject to these constraints, while using the same number of components across all regions and across change and non-change conditions (within each mouse), so as to perform a fair comparison. This gave us on average $53 \\pm 16$ PCA components, with a minimum of 22 and a maximum of 84 PCA components. We find that the results of higher and more sustained redundancy on change flashes continue to hold in the new analysis (Figs. 8,9 in the attached PDF). However, as expected, there is much greater variance across mice, possibly due to variability in the number of PCA components chosen across mice, and due to greater errors in our PID estimates. We will add these new results to the revised supplementary material.\n\nFor consistency across mice, we used a common basis of 10 or 20 PCA components across all mice (and across all regions and conditions). Fig. 6 in the paper had 10 PCA components; Figs. 13 and 14 in the supplementary material of the paper had 20 PCA components. We will also include the rationale for selecting up to 20 PCA components in the revised version of the paper.\n\nApart from this, to show that we can go to higher dimensions, we have also extended our high-dimensional example (Example 10) to a dimensionality of 256 (equivalent to a total of 768 neurons across regions). It should also be noted that other PID methods do not come close (please see the overall rebuttal for comparisons): the discrete PID estimators are limited to low single-digit dimensions, while other methods do not show ground truth validation at high dimensions.\n\n**Questions**\n\n> Can the authors method run on thousands of neurons, consistent with the original motivation they provide? What are the limitations on dimensionality more explicitly for this method? What is the computational complexity? Does the author's method require using PCA or having high firing rate neurons? How does it deal with low firing rate neurons or highly inhomogenous Poisson variables?\n\n1. We have re-run our analysis on as many dimensions as possible in the Allen Institute dataset, as described above. We have also increased the dimensionality in our simulated data experiment and shown that our method works well up to 256 dimensions in M, X and Y each (equivalent to 768 total neurons).\n\n2. The method departs from the ground truth at 512 dimensions; we will identify the reason and provide an update.\n\n3. The computational complexity is O(ND^3) for the various matrix operations, where D=d_M+d_X+d_Y, and N is the number of iterations until convergence. We also provide a timing analysis in Fig 7 (attached PDF).\n\n4. Our method does not require the use of PCA, however, we used PCA to obtain more stable estimates of the covariance matrix, and to minimize bias and error, as described above.\n\n5. Our method also does not require particularly high firing rates: this is demonstrated in the multivariate Poisson simulation, where M1 and M2 had mean “spike-counts” of 2 (see Fig. 4a in the attached PDF). However, we expect that as firing rates increase, the data will appear more Gaussian, as a result of which our method will become more accurate. We will include a mention of this in the revised paper.\n\n6. Inhomogeneous Poisson processes are those whose mean firing- (or “emission-”) rates change with time. Our method for estimating the PID does not consider the temporal characteristics of the original signal. Rather, the data analyst can choose the random variables M, X and Y to span some time range (or different time ranges) as they please. In our example, we counted the number of spikes in a 50-125 ms window after stimulus onset. Even if the underlying spiking process was an inhomogeneous Poisson process, this spike count would be Poisson distributed, with a mean given by the integral of the emission rate over the fixed window. In general, while analyzing data and computing PID values, one will have to be aware of severe inhomogeneities (e.g., if the distribution changes not just in the same way in each trial, but changes across trials), and account for them separately. For example, in our analysis, we excluded time periods when the mice were not “engaged” in the task, as defined by not actively consuming rewards at a rate of at least 2 rewards per minute. We will add a discussion of inhomogeneity, as well as how to choose M, X and Y, to the paper or to the supplementary material."}, {"review_id": "10D0qBOM1X", "reviewer": "Reviewer_DDWz", "summary": "Partial Information Decompositions (PIDs) play an essential role in neuroscience research. One constraint of the broader usability of PIDs is the computational difficulty of computing PIDs for high-dimensional neural recordings. To address this concern, the authors propose a method to compute and estimate a PID efficiently. More specifically, they restrict the optimization space of PIDs to jointly Gaussian, which reduces the number of optimization variables, allowing them to compute PIDs for much higher dimensionalities of neural data. Then, the authors show their method could be written out in closed form and solved by projected gradient descent, and they use nine examples with increasing complexity to show their method has the ability to recover ground truth and stability over increasing dimensionality. The authors also claim it's the first time to correct the bias and variance of estimates in a heuristic way. Finally, the authors evaluate the performance of their method on both synthetic and real neural data.", "questions": "* Does $\\delta$-PID also agree with the ground truth even when distributions are not Gaussian?\n* Could you please discuss how the \"additivity\" property would affect the analysis of communication among brain regions?", "limitations": "I have highlighted technical limitations and weaknesses above. I have nothing further to add here.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "* Computational scalability of ~PID with a basic property called additivity.\n* The authors' exposition of their method and experiments is clear.\n* The first time to raise the issue of bias in PID estimates.", "weaknesses": "* No analysis is done to show the stability of $\\delta$-PID over increasing dimensionality. According to the paper, there are two differences between $\\delta$-PID and $\\sim$-PID. The first one is the \"additivity\" property, which is clearly shown by Examples 8-9 in Section 4. The second one is that $\\sim$-PID uses an exact upper bound, while $\\delta$-PID uses an approximate upper bound.  \nHowever, no analysis is performed to show how the second difference will affect the performance. In other words, people may be curious about the stability of $\\delta$-PID over increasing dimensionality if we use examples with no \"additivity\" property ($\\delta$-PID and $\\sim$-PID both agress with the ground truth).  \nCurrently, it's hard to see the differences between $\\delta$-PID and $\\sim$-PID when applying to high dimensional neural data. I think this point may determine whether the proposed contribution is timely and impactful or a solid technical upgrade without practical consequences in neuroscience.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses:**\n\n> No analysis is done to show the stability of delta-PID over increasing dimensionality. According to the paper, there are two differences between delta-PID and ~-PID. The first one is the \"additivity\" property, which is clearly shown by Examples 8-9 in Section 4. The second one is that ~-PID uses an exact upper bound, while delta-PID uses an approximate upper bound.\n>\n> However, no analysis is performed to show how the second difference will affect the performance. In other words, people may be curious about the stability of delta-PID over increasing dimensionality if we use examples with no \"additivity\" property (delta-PID and ~-PID both agress with the ground truth).\n>\n> Currently, it's hard to see the differences between delta-PID and ~-PID when applying to high dimensional neural data. I think this point may determine whether the proposed contribution is timely and impactful or a solid technical upgrade without practical consequences in neuroscience.\n\nWe thank the reviewer for raising this point. We have added a comparison with the delta-PID for a larger dimensionality (please see Fig. 6,7 in the attached PDF).\n\nJust as the delta-PID fails at a dimensionality of 2, it also fails at 64 (i.e., the error due to additivity extends to higher dimensions). Furthermore, we found that the delta-PID method begins to take a very large amount of time to compute (> 2 hours) at $d_M = d_X = d_Y = 128$, and completely exceeds the memory capacity of a 48GiB workstation at a dimensionality of 256. We have included a timing analysis (please see Fig.7 in the attached PDF) to show the difference in speed of the two methods; our method runs over 1000 times faster at a dimensionality of 64. We will add this analysis to the paper.\n\nIt is much harder to measure the difference between the delta and ~-PIDs due to the difference caused by the delta-PID recovering an approximate upper bound. The cases where the delta-PID agrees with the ~-PID are points where the solution is known in closed form due to a theorem from [12]. In other words, these are “easy” cases. All non-trivial cases are the ones that are constructed using the additivity property. The effective result of this analysis is that the degree to which the ~-PID outperforms the delta-PID due to the fact that the latter’s upper bound is approximate is still unknown. We will mention this point in the revised paper.\n\nHowever, we still believe our method will be impactful and have practical consequences for neuroscience, because it works accurately at much higher dimensions, runs much faster, and due to the importance of the additivity property (as described below).\n\n**Questions:**\n\n> 1. Does delta-PID also agree with the ground truth even when distributions are not Gaussian?\n\nTo the extent that we consider Banerjee et al.’s (ISIT 2018) method to provide ground truth, we show in the paper in Fig. 5 that the delta-PID does not agree with the ground truth to the same degree as the ~-PID for a non-Gaussian distribution. However, it should be noted that Banerjee et al.’s method uses the ~-PID definition. Thus any difference between the delta_G-PID (dashed line) and Banerjee et al.’s method (“ground truth”) could be a result of the difference in the two definitions. Unless we have a more accurate way to compute the delta-PID for discrete distributions, we would not be able to say whether the difference was purely due to the difference in definitions, or due to a difference in accuracy. For additional clarity, we will add this discussion to the paper.\n\n> 2. Could you please discuss how the \"additivity\" property would affect the analysis of communication among brain regions?\n\nAdditivity is an extremely fundamental property: effectively, it states that the PID values of an isolated system should not depend on the PID values of another isolated system. Without additivity, it is not possible to examine systems in isolation, since broadening your view to include a different isolated system could change the PID values of the first system.\n\nHypothetically, if we take two separate individuals receiving completely independent stimuli, and examine the PID between the activity in brain regions M, X and Y in each of their brains, then the effective unique information that X1 and X2 have about M1 and M2 with respect to Y1 and Y2 should be equal to the sum of the unique information in each of the individuals taken separately.\n\nAs another idea, suppose we are trying to examine visual information flow and auditory information flow in a multi-sensory integration task. We may want to understand the degree to which the activity in the sub-regions of the visual and auditory systems depend on each other. A reasonable null model of independence between the two systems would be that the PID values of the joint system will be equal to the sum of the PID values in the two individual systems. Then, measuring the actual degree to which the joint PID value is not equal to the sum will be a meaningful measure of dependence between the systems. This would only be possible with a PID definition that _guaranteed_ additivity of independent sub-systems. However, since the delta-PID does not satisfy the additivity property, we cannot guarantee that the aforementioned null would be the correct null model, and we would not be able to perform such an analysis. We will add this example to the supplementary material."}, {"review_id": "Qr1GgB2Thg", "reviewer": "Reviewer_YrjF", "summary": "The paper proposed a new method for partial information decomposition (PID) on multivariate Gaussian distributions. The issue of bias was discussed, and a correction method was provided. The method was tested on synthetic canonical examples and real data.", "questions": "1. Would non-modeled confounding factors affect the results?\n\nSummary after rebuttal\nThe added simulations for testing higher dimensional data further justify the high score I gave.", "limitations": "The paper does not have negative societal impacts.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The introduction clearly lays out the problem.\n2. The paper is extremely well written with notations clearly defined.\n3. The method extends prior work on PID with new properties. \n4. The method is rigorously tested on simulated data.\n5. The method solves an important problem, namely “the extent to which one region’s activity uniquely explains that of another, while excluding information corresponding to spontaneous behaviors” as stimulus could make two regions seem correlated.", "weaknesses": "1. The paper does not provide enough real data to show effectiveness in real applications.\n2. Testing on higher dimensional settings would be important.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Weaknesses**\n\n> 1. The paper does not provide enough real data to show effectiveness in real applications.\n\nWe thank the reviewer for raising this point. We have tried to ameliorate this concern by demonstrating how the method works with a larger number of neurons (i.e., PCA components) in the same dataset (Fig. 8 and 9 in the attached PDF). We have also demonstrated our method on more simulated non-Gaussian examples, to understand the extent to which our method applies in non-Gaussian settings (Figs. 1-4 in the attached PDF). The focus of our paper was showing that our method works well on examples where the ground truth was known, as well as showing a proof of concept on real data. A more extensive evaluation of our method on real datasets is certainly necessary, but we feel it is beyond the scope of this paper.\n\n> 2. Testing on higher dimensional settings would be important.\n\nWe have increased the dimensionality of the simulated high-dimensional example with known ground truth (Example 10, and we find that our method begins to depart from the ground truth at a dimensionality of 512 (please see the overall rebuttal, and Fig. 5 in the attached PDF)\n\n**Questions**\n\n> Would non-modeled confounding factors affect the results?\n\nOur Gaussian PID depends only on the joint covariance matrix between M, X and Y. Different models of external confounding that result in the same covariance matrix will not affect the PID values obtained. However, if one or more of these confounding variables are included in M, X or Y, and the PID is then re-computed, that may change the PID values completely compared to the case where the confounding variable is not included.\n\nConfounders are also important to consider when interpreting PID values. Our ~_G-PID method is based on the covariance matrix between M, X and Y, and is thus a correlational quantity. We cannot make causal claims based on observed PID values. Many of the same caveats that apply to interpreting correlations would also apply to interpreting observed PID values. A detailed analysis of the connections between different structural causal models (e.g., see Peters et al. 2017, “Elements of Causal Inference”) and their PID profiles has not been undertaken in the literature, to our knowledge, but could be the subject of future work (which is also better enabled by this paper)."}]

Author rebuttal:
We sincerely thank all of the reviewers for taking the time and effort to provide thoughtful and constructive reviews of our work.

We use this space to address a few comments that were common across reviewers, and to explain the additional analyses we perform in response to these comments. We also reiterate the scope and the central contributions of our paper, and describe how these are laid out in the paper.

**Summary of new analyses (and figure numbers in the attached PDF)**

1. Increased dimensionality in Example 10 (Fig 5)
2. Increased number of PCA components in neural data analysis (Fig 8,9). Please see response to reviewer wULj for details.
3. Comparison with delta-PID at higher dimensions (Fig 6,7)
4. More non-Gaussian examples: Binomial, and zero-inflated Binomial and Poisson (Figs. 1-4)

**Scope of the paper**

Our paper provides a new method to compute and estimate the ~-PID definition for multivariate, jointly Gaussian distributions P_MXY. This is done by assuming Gaussian optimality in the optimization problem used to compute the PID, which helps reduce the number of optimization variables. To justify the Gaussian optimality assumption, we consider a number of Gaussian examples with known ground truth, and show that the output of our method agrees with ground truth, even at high dimensionalities. The aforementioned are the subject of Sections 1-5, which deal strictly with Gaussian P_MXY, and endeavor to demonstrate the efficacy of our method on Gaussian distributions. However, since most real data is non-Gaussian, in Section 6, we showed using a low-dimensional simulation, that our method also gives reasonable results on non-Gaussian (Poisson) data. We then apply it to a publicly available neuroscientific dataset to investigate the amount of redundancy between different visual cortical brain regions.

**Extending our method to higher dimensionalities**

In the submitted paper, we presented evidence of agreement with ground truth only up to the dimensions of M, X and Y each being 128 (i.e., a total of 384 dimensions for a covariance matrix of size 384x384). We have now extended this analysis to higher dimensions (Fig. 5 in the attached PDF), showing that our method continues to agree with ground truth up to M, X and Y each being 256-dimensional, i.e., a total of 768 dimensions, for a covariance matrix of size 768x768 (see new figures in attached PDF). Our method begins to deviate from the ground truth at 512 (at which point the covariance matrix has 2.36 million elements). We will investigate the root issue of this problem and provide an update within a week.

It should be noted that our method is the first to provide PID estimates that agree with ground truth at such high dimensions, and that most other works do not come close (see comparisons below).

**More non-Gaussian examples**

In the submitted paper, we only tested our method on a single non-Gaussian, multivariate Poisson example. We have now extended this to include a multivariate Binomial example, as well as zero-inflated versions of the multivariate Poisson and Binomial setups (Fig 1-4 in the attached PDF). We find that our method performs well on the original Poisson and Binomial versions, agreeing with the discrete PID estimator of Banerjee et al. (ISIT 2018).

In the zero-inflated versions, when we have a bimodal distribution with a point-mass at zero, our Gaussian PID estimator no longer recovers the same _absolute_ PID values as the estimator by Banerjee et al.. However, when normalized by the total mutual information, the _relative_ PID-values are still close to the values computed by the Banerjee et al. method and follow the correct trends in every case. Thus, for such a distribution, an improved estimate of the mutual information (e.g., Kraskov et al. 2004; Belghazi et al. ICML 2018) could potentially be used to correct the absolute PID values, after a more thorough evaluation.

The primary challenge with testing non-Gaussian distributions is the absence of good ground truth. We currently rely on a method from Banerjee et al. (2018), designed for discrete distributions, as ground truth. We cannot be sure if the errors seen in the ultimate plot are a result of errors in our method, or errors in computing the ground truth itself.

**Improvements over the delta-PID & other PID measures**

Many reviewers asked us to expand upon how our work improves on prior art. We provide a summary below, which we will include in the Related Work section of our paper. Also, we have now compared our method with the delta-PID at higher dimensions (Fig 6,7 in the attached PDF).

We divide previous works into 3 categories:
1. Compared to Venkatesh et al. (ISIT 2022), here we have:
- Better ground truth validation
- Ability to compute at higher dimensionality (256, rather than 64 or 128)
- Faster computation (>1000x faster at d_M=d_X=d_Y=64)
- Satisfy additivity (which is important; see response to reviewer DDWz)
- Correct for bias in estimates
2. Compared to discrete PID estimators (Banerjee et al., ISIT 2018; Makkeh et al., Entropy 2018), we have
- Much higher dimensionality: both discrete PID estimators were demonstrated for discrete M, X and Y with support sizes of at most 18 in Makkeh et al., which corresponds to a total support of 18^3=5832 for Q_MXY. If we assume a single neuron can have at most 10 spikes (i.e., support of 10), then these methods can handle ~4 neurons, since log_10(5832) < 4.
- We discuss estimation and how bias can be corrected
3. Compared to the continuous PID estimators of Pakman et al. [17] and Liang et al. [21], we have:
- An estimator that matches ground truth at high dimensionalities. [17] shows ground truth validation only in scalar cases, while [21] matches ground truth only for scalar binary cases and deviates from ground truth in other examples.
- We demonstrate our method on higher dimensionalities compared to [17].
- We show bias correction

Meta-review: This paper makes significant contributions to an important topic: that of estimating partial information decompositions, which is a complex task in high dimension.
Compared to previous works, the author's method scales better with dimension and is bias-corrected.
The theory presented in this paper is sound.
All reviewers noted the clarity of the writing of the paper and that of the experiments.

This paper has a positive consensus, as all five reviewers argue for acceptance of the paper.

One reviewer initially argued for the rejection of the paper, yet the authors managed to convince them of the merits of their approach with a very detailed rebuttal. The added experiments run during the discussion period will enhance an already great paper.

In short, the consensus is that this paper will have a practical impact on brain signal processing.
Decision: Accept (spotlight)