Estimating the Unique Information of Continuous Variables
Year: 2021

Abstract:
The integration and transfer of information from multiple sources to multiple targets is a core motive of neural systems. The emerging field of partial information decomposition (PID) provides a novel information-theoretic lens into these mechanisms by identifying synergistic, redundant, and unique contributions to the mutual information between one and several variables. While many works have studied aspects of PID for Gaussian and discrete distributions, the case of general continuous distributions is still uncharted territory. In this work we present a method for estimating the unique information in continuous distributions, for the case of one versus two variables. Our method solves the associated optimization problem over the space of distributions with fixed bivariate marginals by combining copula decompositions and techniques developed to optimize variational autoencoders. We obtain excellent agreement with known analytic results for Gaussians, and illustrate the power of our new approach in several brain-inspired neural models. Our method is capable of recovering the effective connectivity of a chaotic network of rate neurons, and uncovers a complex trade-off between redundancy, synergy and unique information in recurrent networks trained to solve a generalized XOR task.

text: The integration and transfer of information from multiple sources to multiple targets is a core motive of neural systems. The emerging field of partial information decomposition (PID) provides a novel information-theoretic lens into these mechanisms by identifying synergistic, redundant, and unique contributions to the mutual information between one and several variables. While many works have studied aspects of PID for Gaussian and discrete distributions, the case of general continuous distributions is still uncharted territory. In this work we present a method for estimating the unique information in continuous distributions, for the case of one versus two variables. Our method solves the associated optimization problem over the space of distributions with fixed bivariate marginals by combining copula decompositions and techniques developed to optimize variational autoencoders. We obtain excellent agreement with known analytic results for Gaussians, and illustrate the power of our new approach in several brain-inspired neural models. Our method is capable of recovering the effective connectivity of a chaotic network of rate neurons, and uncovers a complex trade-off between redundancy, synergy and unique information in recurrent networks trained to solve a generalized XOR task.

heading: 1 Introduction and background
text: In neural systems, often multiple neurons are driven by one external event or stimulus; conversely multiple neural inputs can converge onto a single neuron. A natural question in both cases is how multiple variables hold information about the singleton variable. In their seminal work [1], Williams and Beer proposed an axiomatic extension of classic information theory to decompose the mutual information between multiple source variables and a single target variable in a meaningful way. For the case of two sources X1, X2, their partial information decomposition (PID) amounts to expressing the mutual information of X1, X2 with a target Y as a sum of four non-negative terms,
I(Y : (X1, X2)) = U(Y :X1\X2) + U(Y :X2\X1) +R(Y : (X1, X2)) + S(Y : (X1, X2)) , (1.1)
corresponding to unique (U1, U2), redundant (R) and synergistic (S) contributions, respectively. These terms should also obey the consistency equations
I(Y :X1) = R(Y : (X1, X2)) + U(Y :X1\X2) , (1.2) I(Y :X2) = R(Y : (X1, X2)) + U(Y :X2\X1) . (1.3)
The PID has proved useful in understanding information processing by distributed systems in a diverse array of fields including machine learning [2, 3], earth science [4] and cellular automata [5],
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
and particularly in neuroscience [6–10], where notions of synergy and redundancy, traditionally considered mutually exclusive and distinguished by the sign of
∆ = I(Y : (X1, X2))− I(Y :X1)− I(Y :X2) , = S(Y : (X1, X2))−R(Y : (X1, X2)) , (1.4)
have long played a central role in the quest to understand how neural circuits integrate information from multiple sources [11–14]. The novelty of the PID framework here is in separating the measures of synergy and redundancy in (1.4).
The above abstract formulation of PID provides three equations for four unknowns, and only becomes operational once one of U1, U2, R, or S is defined. This has been done in [15] via a definition of the unique information:
Definition 1 (BROJA [15]). Given three random variables (Y,X1, X2) with joint probability density p(y, x1, x2), the unique information U1 of X1 with respect to Y is
U(Y :X1\X2) = min q∈Q Iq(Y :X1|X2) , (1.5)
= min q∈Q
∫ dydx1dx2 q(y, x1, x2) log ( q(y, x1|x2)
q(y|x2)q(x1|x2)
) , (1.6)
where Q = {q(y, x1, x2) | q(y, xi) = p(y, xi), i = 1, 2} . (1.7)
In words, we minimize the conditional mutual information I(Y :X1|X2) over the space of density functions that preserve the marginal densities p(y, x1) and p(y, x2). The above definition implies, along with (1.2)-(1.3), that the unique and redundant information only depend on the marginals p(y, x1), p(y, x2), and that the synergy can only be estimated from the full p(y, x1, x2).
The original definition in [15] was limited to discrete random variables. Here, we show that the extension to continuous variables is well-defined and can be practically estimated.
Motivation from decision theory [15]. Consider for simplicity discrete variables. A decision maker DM1 can choose an action a from a finite set A, and receives a reward u(a, y) based on the selected action and the state y, which occurs with probability p(y). Notably, DM1 has no knowledge of y, but observes instead a random signal x1 sampled from p(x1|y). Choosing the action maximizing the expected reward for each x1, his maximal expected reward is
R1 = ∑ x1 p(x1) max a|x1 ∑ y p(y|x1)u(a, y) . (1.8)
DM1 is said to have no unique information about y w.r.t. another decision maker DM2 that observes x2 ∼ p(x2|y) – if R2 ≥ R1 for any set A, any distribution p(y), and any reward function u(a, y). A celebrated theorem by Blackwell [16, 17] states that such a generic advantage by DM2 occurs iff there exist a stochastic matrix q(x1|x2) which satisfies
p(x1|y) = ∑ x2 p(x2|y)q(x1|x2) . (1.9)
But this occurs precisely when the unique information (1.5) vanishes, since then there exists a joint distribution q(y, x1, x2) in Q for which y ⊥ x1|x2, which implies q(x1|x2, y) = q(x1|x2), and thus (1.9) holds. Similar results exist for continuous variables [18, 19]. Thus the unique information from Definition 1 quantifies a departure from Blackwell’s relation (1.9).
In this work we present a definition and a method to estimate the BROJA unique information for generic continuous probability densities. Our approach is based on the observation that the constraints (1.7) can be satisfied with an appropriate copula parametrization, and makes use of techniques developed to optimize variational autoencoders. We only consider one-dimensional Y,X1, X2 for simplicity, but the method can be naturally extended to higher dimensional cases. In Section 2 we review related works, in Section 3 we present our method and Section 4 contains several illustrative examples.

heading: 2 Related works
text: Partial information decomposition offers a solution to a repeated question that was not addressed by ‘classical’ information theory regarding the relations between two sources and a target [1]. From a mathematical perspective a ’functional definition’ has to be made, meaning that such a definition should align with our intuitive notions. Yet, as shown in [20], not all intuitively desirable properties of a PID can be realized simultaneously. Thus, different desirable properties are chosen for distinct application scenarios. Thus, various proposals for decomposition measures are not seen as conflicting but as having different operational interpretations. For example, the BROJA approach used here builds on desiderata from decision theory, while other approaches appeal to game theory [21] or the framework of Kelly gambling [22]. Yet other approaches use arguments from information geometry [23]. Other approaches assume agents receiving potentially conflicting or incomplete information about the source variables for the purpose of inference or decryption (see e.g. [24, 25]). In [26] the authors separate the specific operational interpretations of PID measures from the general structure of information decomposition.
The actual computation of the BROJA unique information is non-trivial, even for discrete variables. Optimization methods exist for the latter case [27–29], and analytic solutions are only known when all the variables are univariate binary [30]. For continuous probability densities, an earlier definition aligned with the BROJA measure was made by Barret [31], but only applies to Gaussian variables. For Barret’s measure, an analytic solution is known when p(y, x1, x2) is a three-dimensional Gaussian density [31], but does not generalize to higher dimensional Gaussians [32].

heading: 3 Bounding and estimating the unique information
text: We proceed in two steps. We first introduce a parametrization of the optimization space Q in (1.7) and then introduce and optimize an upper bound on the unique information.

heading: 3.1 Parametrizing the optimization space with copulas
text: To characterize the optimization space Q in (1.5)-(1.7), it is convenient to recall that according to Sklar’s theorem [33], any n-variate probability density can be expressed as
p(x1 . . . xn) = p(x1) . . . p(xn)c(u1 . . . un) , (3.1)
where p(xi) is the marginal and ui = F (xi) is the CDF of each variable. The dependency structure among the variables is encoded in the function c:[0, 1]n → [0, 1]. This is a copula density, a probability density on the unit hypercube with uniform marginals [34],∫
[0,1]n−1
n∏ j=1,j 6=i duj c(u1 . . . un) = 1 ∀i . (3.2)
Note that under univariate reparametrizations z′i = g(zi), the ui’s and the copula c remain invariant. For an overview of copulas in machine learning, see [35]. Proposition 1. Under the BROJA Definition 1 of unique information, all the terms of the partial information decomposition in (1.1)-(1.3) are independent of the univariate marginals p(x1), p(x2), p(y), and only depend on the copula c(uy, u1, u2).
Proof. Expressing q(y, x1, x2), q(x1, x2), q(y, x2) via copula decompositions (3.1), and changing variables as duy = q(y)dy, etc., the objective function in (1.6) becomes
Iq(Y :X1|X2) = ∫
[0,1]3
duydu1du2 c(uy, u1, u2) log
( c(uy, u1, u2)
c(uy, u2)c(u1, u2)
) . (3.3)
Note that the copula of any marginal distribution is the marginal of the copula:
c(uy, u2) = ∫ [0,1] du1 c(uy, u1, u2) , c(u1, u2) = ∫ [0,1] duy c(uy, u1, u2) . (3.4)
Thus the optimization objective and the unique information are independent of the univariate marginals. A similar result holds for the mutual information terms in the l.h.s. of (1.1)-(1.3).1 It follows that none of the PID terms in (1.1)-(1.3) depend on the univariate marginals, and therefore all the PID terms are invariant under univariate reparametrizations of (y, x1, x2).
In order to parametrize the optimization space Q in (1.7) using copulas, consider the factorization
p(y, x1, x2) = p(x1)p(y|x1)p(x2|y, x1) . (3.5)
Using the copula decomposition (3.1) for n = 2, the last two factors in (3.5) can be expressed as
p(y|x1) = p(y, x1) p(x1) = p(y)p(x1)c(y, x1) p(x1) = c(uy, u1)p(y) , (3.6)
and similarly
p(x2|y, x1) = p(x1, x2|y) p(x1|y) , (3.7)
= c1,2|y(u1|y, u2|y)p(x2|y) , (3.8) = c1,2|y(u1|y, u2|y)c(uy, x2)p(x2) , (3.9)
where we defined the conditional CDFs,
ui|y = F (ui|uy) = ∂C(uy, ui)
∂uy i = 1, 2 (3.10)
and C(uy, ui) is the CDF of c(uy, ui). Note that the function c1,2|y(u1|y, u2|y) in (3.8) is not the conditional copula c(u1, u2|uy), but rather the copula of the conditional p(x1, x2|y). Using expressions (3.6) and (3.9), the full density (3.5) becomes
p(y, x1, x2) = p(y)p(x1)p(x2)c(uy, u1, u2) , (3.11)
where
c(uy, u1, u2) = c(uy, u1) c(uy, u2)c1,2|y(u1|y, u2|y) . (3.12)
This is a simple case of the pair-copula construction of multivariate distributions [38–40], which allows to expand any n-variate copula as a product of (conditional) bivariate copulas.
Proposition 2. The copula of the conditional, c1,2|y(·, ·), parametrizes the space Q in (1.7).
Proof. Since q(y, xi) = p(y, xi) (i = 1, 2), the copula factors in
p(y, xi) = p(y)p(xi) c(uy, ui) , i = 1, 2 (3.13)
are fixed in Q. Therefore, in the copula decomposition (3.12) for q(y, x1, x2) ∈ Q, only the last factor can vary in Q. Let us denote by θ the parameters of a generic parametrization for the copula c1,2|y(u1|y, u2|y). Since the latter is conditioned on uy, the parameters can be taken as a function θ(uy). It follows that the copula of q necessarily has the form
cθ(uy, u1, u2) = c(uy, u1) c(uy, u2) c1,2|θ(uy)(u1|y, u2|y) , (3.14)
and the parameters of the function θ(uy) are the optimization variables.2
1The connection between mutual information and copulas was discussed in [36, 37]. 2We note that in multivariate pair-copula expansions it is common to assume constant conditioning parame-
ters θ [41], but we do not make such a simplifying assumption.

heading: 3.2 Optimizing an upper bound
text: Inserting now the expression (3.14) into the objective function (3.3) we get I[θ] = Ecθ(uy,u1,u2) log [ c(uy, u1)c1,2|θ(uy)(u1|y, u2|y) ] − Ecθ(u1,u2) log cθ(u1, u2) , (3.15)
which is our objective function and satisfies the marginal constraints (1.7). Note that apart from the optimization parameters θ, it depends on the bivariate copulas c(uy, u1) and c(uy, u2) which should be estimated from the observed data. Given D observations {y(i), x(i)1 , x (i) 2 )}Di=1, we map each value to [0, 1] via the empirical CDFs of each coordinate (y, x1, x2). Computing the latter has a O(D logD) cost from sorting each coordinate and yields a data set {u(i)y , u(i)1 , u (i) 2 )}Di=1. The latter set is used to estimate copula densities c(uy, u1) and c(uy, u2) by fitting several parametric and non-parametric copula models [42], and choosing the best pair of models using the AIC criterion.3 From the learned copulas we also get the conditional CDF functions ui|y = F (ui|uy) that appear in the arguments of the first term in (3.15).
A variational upper bound. Minimizing (3.15) directly w.r.t. θ is challenging because the second term depends on the copula marginal cθ(u1, u2) which has no closed form, as it requires integrating (3.14) w.r.t. uy. We introduce instead an inference distribution rφ(uy|u1, u2), with parameters φ, that approximates the conditional copula cθ(uy|u1, u2), and consider the bound
log cθ(u1, u2) = log ∫ du′y cθ(u ′ y, u1, u2) ≥ ∫ du′y rφ(u ′ y|u1, u2) log cθ(u ′ y, u1, u2)
rφ(u′y|u1, u2) , (3.16)
which follows from Jensen’s inequality and is tight when rφ(u′y|u1, u2) = cθ(u′y|u1, u2). This expression gives an upper bound on Iq[θ], which can be minimized jointly w.r.t. (θ, φ).
A disadvantage of the bound (3.16) is that its tightness depends strongly on the expressiveness of the inference distribution rφ(u′y|u1, u2). This situation can be improved by considering a multiple-sample generalization proposed by [44],
log cθ(u1, u2) ≥ DA,θ,φ(u1, u2) ≡ Ep(u(1)y ...u(A)y ) log
[ 1
A A∑ a=1 cθ(u (a) y , u1, u2) rφ(u (a) y |u1, u2)
] , (3.17)
where the expectation is w.r.t. A independent samples of rφ(u′y|u1, u2). DA,θ,φ(u1, u2) coincides with the lower bound in (3.16) for A = 1 and satisfies [44]
DA+1,θ,φ(u1, u2) ≥ DA,θ,φ(u1, u2), (3.18) lim A→∞ DA,θ,φ(u1, u2) = log cθ(u1, u2) . (3.19)
Thus, even when rφ(u′y|u1, u2) 6= cθ(u′y|u1, u2), the bound can be made arbitrarily tight for large enough A. Inserting (3.17) in (3.15), we get finally
Iq[θ] ≤ B1[θ] +B2[θ, φ] , (3.20)
where
B1[θ] = Ecθ(uy,u1,u2) log [ c(uy, u1)c1,2|θ(uy)(u1|y, u2|y) ] , (3.21)
B2[θ, φ] = −Ecθ(u1,u2)DA,θ,φ(u1, u2) , (3.22)
and we minimize the r.h.s. of (3.20) w.r.t. (θ, φ). Low-variance estimates of the gradients to perform the minimization can be obtained with the reparametrization trick [45, 46], as discussed in detail in Appendix A. In our examples below we use for c1,2|θ(uy) a bivariate Gaussian copula (reviewed in Appendix B). Such a copula has just one parameter θ ∈ [−1,+1], and thus the optimization is done over the space of functions θ(uy):[0, 1]→ [−1,+1], which we parametrize with a two-layer neural network. Similarly, we parametrize rφ(uy|u1, u2) with a two-layer neural network. Details of these networks are in Appendix D.
3For this fitting/model selection step, we used the pyvinecopulib python package [43].
While the term B2 in our bound is similar to the negative of the ELBO bound in importance weighted autoencoders (IWAEs) [44], there are some differences between the two settings, the most important being that we are interested in the precise value of the bound at the minimum, rather than the learned functions cθ, rφ. Note also that our latent variables u (k) y are one-dimensional, as opposed to the usual higher dimensional latent distributions of variational autoencoders, and that the empirical expectation over data observations in IWAEs is replaced in B2 by the expectation over cθ(u1, u2), whose parameters are also optimized.
Estimating the other PID terms In the following we adopt the minimal value taken by the upper bound (3.20) as our estimate of U1. The other terms in the partial information decomposition are obtained from the consistency relations (1.1)-(1.3), after estimating the mutual informations I(Y : (X1, X2)), I(Y : X1), I(Y : X2). There are several methods for the latter. In our examples, we use the observed data to fit additional copulas c(u1, u2) and c12|θ(uy) and estimate I(Y :X1) ' 1 D ∑D i=1 log c(u (i) y , u (i) 1 ) and similarly for the other terms. Note that all our estimates have sources of potential bias. Firstly, the estimation of the parametric copulas is subject to model or parameter misspecification, which can be ameliorated by more refined model selection strategies. Secondly, the optimized bound might not saturate, biasing the estimate upwards. This can be improved using higher A values and improving the gradient-based optimizer used.

heading: 4 Examples
text: Comparison with exact results for Gaussians. Consider a three-dimensional Gaussian with correlations ρy,xi between y, xi for i = 1, 2. The exact solution to (1.5) in this case is [31]
U(Y :X1\X2) = 1
2 log ( 1− ρ2y,x2 1− ρ2y,x1 ) 1 [ρy,x2 < ρy,x1 ] . (4.1)
Fig. 1 compares the above expression with estimates from our method. Here we know that cy,1 and cy,2 are Gaussian copulas, with parameters ρy,x1 , ρy,x2 , and we assumed a Gaussian copula for c1,2|y,θ(u1|y, u2|y) as well. For each pair of values ρy,x1 , ρy,x2 . In this and the rest of the experiments, we optimized the parameters (θ, φ) using the ADAM algorithm [47] with a fixed learning rate 10−2 during 1200 iterations, and using A = 50. The results reported correspond to the mean of the bound in the last 100 iterations. The comparison in Fig. 1 shows excellent agreement.
Model systems of three neurons. The nature of information processing of neural systems is a prominent area of application of the PID framework, since synergy has been proposed as natural measure of information modification [7, 48]. We consider two models:
M1 (X1, X2) ∼ N (0, ρ212),
Y = tanh(w1X1 + w2X2).
M2 (X1, X2) ∼ N (0, ρ212), Y = X21/ ( 0.1 + w1X 2 1 + w2X 2 2 ) .
(4.2)
Both models are parameterized by the correlation ρ12 and weights w1, w2. Model 1 is a particularly simple neural network. The tanh activation does not affect its copula, and even for a linear activation function the variables are not jointly Gaussian since Y is deterministic on (X1, X2). Model 2 is inspired by a normalization operation widely believed to be canonical in neural systems [49] and plays a role in common learned image compression methods [50]. The results, presented in Figure 2. are obtained from 3000 samples from each model
Computational aspects of connectivity in recurrent neural circuits. We apply our continuous variable PID to understand computational aspects of the information processing between recurrently coupled neurons (Fig. 3). A large amount of work has been devoted to applying information theoretic measures for quantifying directed pairwise information transfer between nodes in dynamic networks and neural circuits [51]. However, classical information theory only allows for the quantification of information transfer, whereas the framework of PID enables further decomposition of information processing into transfer, storage, and modification, providing further insights into the computation within a recurrent system [52]. Transfer entropy (TE) [53] is a popular measure to estimate the directed transfer of information between pairs of neurons [54, 55], and is sometimes approximated by linear Granger causality. Intuitively, TE between a process X and a process Y measures how much the past of X , X−, can help to predict the future of Y , Y +, accounting for its past Y −. Although TE quantifies how much information is transferred between neurons, it does not shed light on the computation emerging from the interaction of X− and Y −. Simply put, the information transferred from X− could enter Y +, independently of the past state Y −, or it could be fused in a non-trivial way with the information in the state in Y −[52, 56]. PID decomposes the TE into modified transfer (quantified by S(Y +:X−, Y −)) and unique transfer (quantified by U(Y +:X− \ Y −)) terms (see the Appendix for a proof):
TE(X → Y ) = I(Y +:X−|Y −) = U(Y +:X− \ Y −) + S(Y +:X−, Y −) .
Furthermore, the information kept by the system through time can be quantified by the unique storage (given by U(Y +:Y − \X−)) and redundant storage (given by R(Y +:X−, Y −)) in PID [48]. This perspective is a new step towards understanding how the information is processed in recurrent systems beyond merely detecting the direction functional interactions estimated by traditional TE methods (see Appendix G, for details). To explore these ideas, we simulated chaotic networks of rate neurons with an a-priori causal structure consisting of two sub-networks X and Y (Fig. 3a, see [57] for more details on causal analyses of this network model). The sub-network X is a Rossler attractor of three neurons obeying the dynamical equations:
Ẋ1 = −X2 −X3 Ẋ2 = X1 + αX2 Ẋ3 = β +X3(X1 − γ)
(4.3)
where {α, β, γ} = {0.2, 0.2, 5.7}. There are 100 neurons in the sub-network Y from which we chose the first three, Y1:3, to simulate the effect of unobserved nodes. Neurons within the sub-network Y obey the dynamical equations
Ẏ = −λY + 10 tanh(JY XX + JY Y Y ) (4.4) where JY X ∈ R100×3 has all its entries equal to 0.1, and JY Y is the recurrent weight matrix of the Y sub-network, sampled as zero-mean, independent Gaussian variables with standard deviation g = 4. No projections exist from the downstream sub-network Y to the upstream sub-network X. We simulated time series from this network (exhibiting chaotic dynamics, see Fig. 3a) and estimated the PID as unique, redundant, and synergistic contribution of neuron i and neuron j at time t in shaping the future of neuron j at time t + 1. For each pair of neurons Zi, Zj ∈ {X1:3, Y1:3} we treated (Zti , Z t j , Z t+1 j ) T t=1 as iid samples
4 and ran PID on these triplets (i, j represent rows and columns in Fig. 3b-d). The PID uncovered the functional architecture of the network and further revealed non-trivial interactions between neurons belonging to the different sub-networks, encoded in four matrices: modified transfer S, unique transfer U1, redundant storageR, and unique storage U2 (details in Fig. 3d). The sum of the modified and unique transfer terms was found to be consistent with the TE (Fig. 3c, TE equal to S + U1, up to estimation bias). The TE itself captured the network effective connectivity, consistent with previous results [55, 57].
Uncovering a plurality of computational strategies in RNNs trained to solve complex tasks. A fundamental goal in neuroscience is to understand the computational mechanisms emerging from
4Note that the estimation of the PID from many samples of the triplets (Zti , Z t j , Z t+1 j ) is operationally the same whether such triplets are iid or, as in our case, temporally correlated. This is similar to estimating expectations w.r.t. the equilibrium distribution of a Markov chain by using temporally correlated successive values of the chain. In both cases, the temporal correlations do not introduce bias in the estimator but can increase the variance.
the collective interactions of recurrent neural circuits leading to cognitive function and behavior. Here, we show that PID opens a new window for assessing how specific computations arise from recurrent neural interactions. Unlike MI or TE, the PID quantifies the alternative ways in which a neuron determines the information in its output from its inputs, and thus can be a sensitive marker of different computational strategies. We here trained RNNs as models of cortical circuits [59] and used the PID to elucidate how the computations emerging from recurrent neural interactions contribute to task performance. We trained RNNs to solve a generalized version of the classic XOR classification problem with target labels corresponding to odd vs. even mixture components (Fig. 4a). Stimuli were presented for one time step (t = 0) and the network was trained to report the decision at t = 10. By tracking the temporal trajectories of the hidden layer activity we found that the network recurrent dynamics (represented as unfolded in time in Fig. 4b) progressively pulls the two input classes in opposite directions along the output weights (see Appendix). We used PID to dissect how a plurality of different strategies emerge from recurrent neural interactions in RNNs trained for solving a classification task. The computation emerged from the recurrent interaction between hidden neurons at different time steps. Do all successfully trained networks have a similar profile in terms of the PID terms? If so, this hints at a single computational strategy across these networks. If not, it is safe to assume that task performance is reached via different mechanisms, despite identical network architecture and training algorithm.
We found that on average across multiple networks S, R, and U rose with task difficulty (Fig. 4c), yet at all difficulties, individual networks differed strongly with respect to the ratio S/U , i.e. there were networks with larger average synergy across neuron pairs compared to the average unique information, and vice versa. For simple networks like the ones used here, one can inspect receptive fields to understand the reason for this differential behaviour (Fig. 4d-e). Indeed, networks with high average unique information displayed ’grandmother-cell’-like neurons, that would alone classify a large parts of the sample space, while in networks with higher average synergy such cells were absent (Fig. 4d). The emergence of these ’grandmother-cell’-like receptive fields is due to the recurrent dynamics. While in a feedforward architecture (WXX = 0) hidden layer receptive fields are captured by hyperplanes in input space, in the RNN the receptive fields are time dependent, where later times are interpreted as deeper layers (Fig. 4b) and thus can capture highly non-linear features in input space. The advantage of PID versus a manual inspection of receptive fields is twofold: First, the PID framework abstracts and generalizes descriptions of receptive fields as being e.g. ’grandmother-cell’-like; thus the concept of unique information stays relevant even in scenarios where the concept of a receptive field becomes meaningless, or inaccessible. Second, the quantitative outcomes of a PID rest only on information theory, not specfic assumptions about neural coding or computational strategies, and can be obtained for large numbers of neurons.
Comparison of our PID-based approach with the concept of neuronal selectivity used in neuroscience highlights interesting similarities and differences. Several kinds of selectivity (pure, mixed linear, and mixed non-linear) can be identified by performing regression analysis of neural responses vs. task variables [60]. In this framework, our grand-mother cells correspond to neurons with pure selectivity to the input class labels (a.k.a. "choice-selective" neurons). In the XOR task, [60] showed that non-linear mixed selectivity of neurons to the class labels is beneficial when solving the XOR task, by leading to a high-dimensional representation of the task variables. While selectivity profiles are a property of single neuron responses to task variables, our PID measures are a property of the combined activity of triplets of neurons and thus reveal emerging functional interactions between units and their computational algorithms (see also [7] and [52]). This allowed us to characterize a functional property of neural systems less studied than task variable selectivity: the computations that require functional mixing of the information from multiple units (measured by the average synergistic information) vs. the computations that rely on the output of individual neurons (measured by the unique information and described as grandmother cells). Concretely, by comparing PID and receptive fields we found that that in networks with high unique information, neurons typically have receptive fields with pure selectivity (grandmother cells, with large unique information to the class labels). In networks with high synergy, neurons show complex mixed selectivity to class labels.

heading: 5 Conclusions
text: We presented a partial information decomposition measure for continuous variables with arbitrary probability densities, thereby extending the popular BROJA PID measure for discrete variables. Extending PID measures to continuous variables drastically broadens the possible applications of the PID framework. This is important as the latter provides key insights into the way a complex system represents and modifies information in a computation – via asking which variables carry information about a target uniquely (such that it can only be obtained from that variable), redundantly, or only synergistically with other variables. Answering these questions is pivotal to understanding distributed computation in complex systems in general, and neural coding in particular. We believe that the methods presented here will allow PIDs to be extended efficiently in neuroscience for multiple continuous sources with potentially complex dependency structures, as would be common in cellular imaging data or activation properties of brain modules or areas in functional imaging. More generally, the approach we presented here would be relevant for other application domains such as machine learning, biomedical science, finance, and the physical sciences.

heading: Acknowledgments
text: We thank Thibault Vatter and Praveen Venkatesh for conversations. The work of AP is supported by the Simons Foundation, the DARPA NESD program, NSF NeuroNex Award DBI1707398 and The Gatsby Charitable Foundation. DG is supported by a Swartz Fellowship. AM and MW are supported by Volkswagenstiftung under the program ‘Big Data in den Lebenswissenschaften’ and by the Ministry for Science and Education of Lower Saxony and the Volkswagen Foundation through the ‘Niedersächsisches Vorab’. LM is supported by NINDS Grant NS118461 (BRAIN Initiative). ES is supported by the Simons Collaboration on the Global Brain (542997) as well as research support from Martin Kushner Schnur and Mr. and Mrs. Lawrence Feis, and is the Joseph and Bessie Feinberg Professorial Chair.

Review summaries:
(I am not an expert in information theory) The authors extend partial information decomposition to the case of continuous variables with arbitrary distributions.  My review will focus on the last section, where a recurrent neural network was trained on a generalized XOR task. By decomposing the unique, redundant and synergistic components of each neuron’s contribution to the output, the authors identify different strategies used by these networks.  

This paper extends the popular PID measure of Bertschinger et al. (2014) to continuous random variables, and provides a method to compute and estimate the PID in this case. The authors use a novel method based on copulas to reparameterize the optimization problem used to compute the PID, and estimate it by using an approximation based on a variational upper bound. They also provide several examples to show that their estimates match known PID quantities and that the PID can be useful in various neuroscientifically relevant settings.

Authors present a partial information decomposition (PID) method that can be applied to continuous variables thereby extending a popular approach for analysis of discrete variables. Authors validate their method  by leveraging known analytic results for Gaussian random variables. Additionally, the authors demonstrate that their approach is able to recover the effective connectivity of a neuronal network with chaotic dynamics. Finally, authors use their method to investigate a complex trade-off between redundant, synergetic, and unique information components in recurrent networks that have been trained to solve the XOR problem.

   This paper presents a novel method for computing the Partial    Information Decomposition (BROJA flavour) in the case of three    continuous random variables. The method uses a variational approach    on a lower bound derived from a copula representation of the joint    PDF of the three variables. The paper presents some application    examples, including two where the proposed method is used to    elucidate the properties of recurrent neural networks.