Title: DATA-EFFICIENT MUTUAL INFORMATION NEURAL ESTIMATOR FOR STATISTICAL DEPENDENCY TESTING Anonymous

Abstract: Measuring Mutual Information (MI) between high-dimensional, continuous, random variables from observed samples has wide theoretical and practical applications. Recent works have developed accurate MI estimators through provably lowbias approximations and tight variational lower bounds assuming abundant supply of samples, but require an unrealistic number of samples to guarantee statistical significance of the estimation. In this work, we focus on improving data efficiency and propose a Data-Efficient MINE Estimator (DEMINE) that can provide a tight lower confident interval of MI under limited data, through adding cross-validation to the MINE lower bound (Belghazi et al., 2018) . Hyperparameter search is employed and a novel meta-learning approach with task augmentation is developed to increase robustness to hyperparameters, reduce overfitting and improve accuracy. With improved data-efficiency, our DEMINE estimator enables statistical testing of dependency at practical dataset sizes. We demonstrate the effectiveness of DEMINE on synthetic benchmarks and real world fMRI data, with application of inter-subject correlation analysis.

1. INTRODUCTION:

Mutual Information (MI) is an important, theoretically grounded measure of similarity between random variables. MI captures general, non-linear, statistical dependencies between random variables. MI estimators that estimate MI from samples are important tools widely used in not only subjects such as physics and neuroscience, but also machine learning ranging from feature selection and representation learning to explaining decisions and analyzing generalization of neural networks.

Existing studies on MI estimation between general random variables focus on deriving asymptotic lower bounds and approximations to MI under infinite data, and techniques for reducing estimator bias such as bias correction, improved signal modeling with neural networks and tighter lower bounds. Widely used approaches include the k-NN-based KSG estimator (Kraskov et al., 2004) and the variational lower-bound-based Mutual Information Neural Estimator (MINE) family (Belghazi et al., 2018; Poole et al., 2018) .

Despite the empirical and asymptotic bias improvements, MI estimation has not seen wide adoption. The challenges are two-fold. First, the analysis of dependencies among variables -let alone any MI analyses for scientific studies -requires not only an MI estimate, but also confidence intervals (Holmes & Nemenman, 2019) around the estimate to quantify uncertainty and statistical significance. Existing MI estimators, however, do not provide confidence intervals. As low probability events may still carry a significant amount of information, the MI estimates could vary greatly given additional observations (Poole et al., 2018) . Towards providing upper and lower bounds of true MI under limited number of observations, existing MI lower bound techniques assume infinite data and would need further relaxations when a limited number of observations are provided. Closest to our work, Belghazi et al. (2018) studied the lower bound of the MINE estimator under limited data, but it involves bounds on generalization error of the signal model and would not yield useful confidence intervals for realistic datasets. Second, practical MI estimators should be insensitive to the choice of hyperparameters. An estimator should return a single MI estimate with its confidence interval irrespective of the type of the data and the number of observations. For learning-based approaches, this means that the model design and optimization hyperparameters need to not only be determined automatically but also taken into account when computing the confidence interval.

Towards addressing these challenges, our estimator, DEMINE, introduces a predictive MI lower bound for limited samples that enables statistical dependency testing under practical dataset sizes. Our estimator builds on top of the MINE estimator family, but performs cross-validation to remove the need to bound generalization error. This yields a much tighter lower bound agnostic to hyperparameter search. We automatically selected hyperparameters through hyperparameter search, and a new cross-validation meta-learning approach is developed, based upon few-shot meta-learning, to automatically decide initialization of model parameters. Meta-overfitting is strongly controlled through task augmentation, a new task generation approach for meta-learning. With these improvements, we show that DEMINE enables practical statistical testing of dependency for not only synthetic datasets but also for real world functional Magnetic Resonance Imaging (fMRI) data analysis capturing nonlinear and higher-order brain-to-brain coupling.

Our contributions are summarized as follows: 1) A data-efficient Mutual Information Neural Estimator (DEMINE) for statistical dependency testing; 2) A new formulation of meta-learning using Task Augmentation (Meta-DEMINE); 3) Application to real life, data-scarce applications (fMRI).

2. RELATED WORK:



2.1. MI ESTIMATION:

A widely used approach for estimating MI from samples is using k-NN estimates, notably the KSG estimator (Kraskov et al., 2004) . Gao et al. (2017) provided a comprehensive review and studied the consistency and of asymptotic confidence bound of the KSG estimator (Gao et al., 2018) . MI estimation can also be achieved by estimating individual entropy terms through kernel density estimation (Ahmad & Lin, 1976) or cross-entropy (McAllester & Statos, 2018) . Despite their good performance on random variables with few dimensions, MI estimation on high-dimensional random variables remains challenging for commonly used Gaussian kernels. Fundamentally, estimating MI requires accurately modeling the random variables, where high-capacity neural networks have shown excellent performance on complex high-dimensional signals such as text, image and audio.

Recent works on MI estimation have focused on developing tight asymptotic variational MI lower bounds where neural networks are used for signal modeling. The IM algorithm (Agakov, 2004) introduces a variational MI lower bound, where a neural network q(z|x) is learned as a variational approximation to the conditional distribution P (Z|X). The IM algorithm requires the entropy, H(Z), and E XZ log q(z|x) to be tractable, which applies to latent codes of Variational Autoencoders (VAEs) and Generative Adversarial Networks (GANs) as well as categorical variables. Belghazi et al. (2018) introduces MI lower bounds MINE and MINE-f which allow the modeling of general random variables and shows improved accuracy for high-dimensional random variables, with application to improving generative models. Poole et al. (2018) introduces a spectrum of energy-based MI estimators based on MINE and MINE-f lower bounds and a new TCPC estimator inspired by Contrastive Predictive Coding (Oord et al., 2018) for the case when multiple samples from P (Z|X) can be drawn.

Our work introduces cross-validation to the MINE-f estimator. We derive the lower bound of MINE-f under limited number of samples, and introduce meta-learning and hyperparameter search to enable practical statistical dependency testing.

2.2. GENERAL STATISTICAL DEPENDENCY TESTING:

Existing works in general statistical dependency testing (Bach & Jordan, 2002; Gretton et al., 2005a; Berrett & Samworth, 2019) have developed non-parametric independent criterions based on correlation and mutual information estimators equivalent to testing I(X; Z) = 0, followed by detailed bias and variance analyses. Our approach for independent testing suggest a different direction by harnessing the generalization power of neural networks and may improve test performance on complex signals. The p-values provided by our test do not involve approximated distributions and hold for small number of examples and arbitrary number of signal dimensions. As different statistical dependency testing approaches have explicit or implicit assumptions and biases that make them suitable in different situations, a fair comparison across different approaches is a challenging task. Instead, we focus on a self-contained presentation of our dependency test, and provide preliminary comparisons with a widely studied Hilbert-Schmidt independence criterion (HSIC) (Gretton et al., 2005a) in the appendix.

2.3. META LEARNING:

Meta-learning, or "learning to learn", seeks to improve the generalization capability of neural networks by searching for better hyperparameters (Maclaurin et al., 2015) , network architectures (Pham et al., 2018) , initialization (Finn et al., 2017a; 2018; Kim et al., 2018) and distance metrics (Vinyals et al., 2016; Snell et al., 2017) . Meta-learning approaches have shown significant performance improvements in applications such as automatic neural architecture search (Pham et al., 2018) , few-shot image recognition (Finn et al., 2017a) and imitation learning (Finn et al., 2017b) .

In particular, our estimator benefits from the Model-Agnostic Meta-Learning (MAML) (Finn et al., 2017a) framework which is designed to improve few-shot learning performance. A network initialization is learned to maximize its performance when fine-tuned on few-shot learning tasks. Applications include few-shot image classification and navigation.

We leverage the model-agnostic nature of MAML for MI estimation between generic random variable and adopt MAML for maximizing MI lower bounds. To construct a collection of diverse tasks for MAML learning from limited samples, inspired by MI's invariance to invertible transformations, we propose a task-augmentation protocol to automatically construct tasks by sampling random transformations to transform the samples. Results show reduced overfitting and improved generalization.

3. BACKGROUND:

In this section, we will provide the background necessary to understand our approach1 . We define X and Z to be two random variables, P (X, Z) is the joint distribution, and P (X) and P (Z) are the marginal distributions over X and Z respectively. Our goal is to estimate MI, I(X; Z) given independent and identically distributed (i.i.d.) sample pairs (x i , z i ), i = 1, 2 . . . n from P (X, Z). Let F = {T θ (x, z)} θ∈Θ be a class of scalar functions, where θ is the set of model parameters. Let q(x|z) = p(x) x,z) . Results from previous works (Belghazi et al., 2018; Poole et al., 2018) show that the following energy-based family of lower bounds of MI hold for any θ:

e T θ (x,z) E (x,z)∼P XZ e T θ (

I(X; Z) ≥ E (x,z)∼P XZ log q(x|z) p(x) = E (x,z)∼P XZ T θ (x, z) -Ex∼P X log Ez∼P Z e T θ (x,z) IEB1 ≥ E (x,z)∼P XZ T θ (x, z) -log Ex∼P X ,z∼P Z e T θ (x,z) IMINE ≥ E (x,z)∼P XZ T θ (x, z) -Ex∼P X ,z∼P Z e T θ (x,z) + 1 IMINE-f,IEB (1)

where, E is the expectation over the given distribution. Based on I MINE , the MINE estimator I(X, Z) n is defined as in Eq.2. Estimators for I EB1 , I MINE-f and I EB can be defined similarly.

EQUATION

With infinite samples to approximate expectation, Eq.2 converges to the lower bound I(X, Z) ∞ = sup θ∈Θ I MINE . Note that the number of samples n needs to be substantially more than the number of model parameters d = |θ| to guarantee that T θ (X, Y ) does not overfit to the samples (x i , z i ), i = 1, 2 . . . n and overestimate MI. Formally, the sample complexity of MINE is defined as the minimum number of samples n in order to achieve Eq.3,

EQUATION

Specifically, MINE proves that under the following assumptions: 1)

T θ (X, Z) is L-Lipschitz; 2) T θ (X, Z) ∈ [-M, M ], 3) {θ i ∈ [-K, K],

∀i ∈ 1, . . . , d}, the sample complexity of MINE is given by Eq.4.

EQUATION

)

For example, a neural network with dimension d = 10, 000, M = 1, K = 0.1 and L = 1, achieving a confidence interval of = 0.1 with 95% confidence (δ = 0.05) would require n ≥ 18, 756, 256 samples. This is achievable for synthetic example generated by GANs like that studied in Belghazi et al. (2018) . For real data, however, the cost of data acquisition for reaching statistically significant estimation can be prohibitively expensive. Our approach instead uses the MI lower bounds specified in Eq.1 from a prediction perspective, inspired by cross-validation. Our estimator, DEMINE, improves sample complexity by disentangling data for lower bound estimation from data for learning a generalizable T θ (X, Z). DEMINE enables high-confidence MI estimation on small datasets.

4. APPROACH:

Section 4.1 specifies DEMINE for predictive MI estimation and derives the confidence interval; Section 4.2 formulates Meta-DEMINE, explains task augmentation, and defines the optimization algorithms.

4.1. PREDICTIVE MUTUAL INFORMATION ESTIMATION:

In DEMINE, we interpret the estimation of MINE-f lower bound2 Eq.1 as a learning problem. The goal is given a limited number of samples, infer the optimal network T θ * (X, Z) with parameters θ * defined as follows:

θ * = arg max θ∈Θ E P XZ T θ (X, Z) -E P X E P Z e T θ (X,Z) + 1.

Specifically, samples from P (X, Z) are subdivided into a training set {(x i , z i ) train , i = 1, . . . , m} and a validation set {(x i , z i ) val , i = 1, . . . , n}. The training set is used for learning a network θ as an approximation to θ * whereas the validation set is used for computing the DEMINE estimation I(X, Z) n, θ defined as in Eq.5.

EQUATION

We propose an approach to learn θ, DEMINE. DEMINE learns θ by maximizing the MI lower bound on the training set as follows:

θ = arg min θ∈Θ L({(x, z)} train , θ), where, L({(x, z)}B, θ) = - 1 |B| |B| i=1 T θ (xi, zi)B + 1 |B| 2 |B| i=1 |B| j=1 e T θ (x i ,z j ) B -1. (6)

The DEMINE algorithm is shown in Algorithm 2 in appendix.

Sample complexity analysis.:

Because θ is learned independently of validation samples {(x i , z i ) val , i = 1, . . . , n}, the sample complexity of the DEMINE estimator does not involve the model class F and the sample complexity is greatly reduced compared to MINE-f. DEMINE estimates I(X, Z) ∞, θ when infinite number of samples are provided, defined as:

EQUATION

We now derive the sample complexity of DEMINE defined as the number of samples n required for I(X, Z) n, θ to be a good approximation to I(X, Z) ∞, θ in Theorem 1.

Theorem 1. For T θ (X, Z) bounded by [L, U ], given any accuracy and confidence δ, we have:

Pr(|I(X, Z) n, θ -I(X, Z) ∞, θ | ≤ ) ≥ 1 -δ

when the number of validation samples n satisfies:

n ≥ n * , s.t. f (n * ) ≡ min 0≤ξ≤ 2e -2ξ 2 n * (U -L) 2 + 4e - ( -ξ) 2 n * 2(e U -e L ) 2 = δ (8) Proof. Since T θ (X, Z) is bounded by [L, U ],

applying the Hoeffding inequality to the first half of Eq.5 yields:

Pr(| 1 n n i=1 T θ (x i , z i ) -E P XZ T θ (X, Z)| ≥ ξ) ≤ 2e -2ξ 2 n (U -L) 2

As e T θ (X,Z) is bounded by [e L , e U ], applying the Hoeffding inequality twice to the second half of Eq.5:

Pr(|EP X EP Z e T θ (X,Z) -1 n n i=1 EP Z e T θ (x i ,z) | ≥ ζ) ≤ 2e - 2ζ 2 n (e U -e L ) 2 Pr(|EP Z 1 n n i=1 e T θ (x i ,z) -1 n n j=1 1 n n i=1 e T θ (x i ,z j ) | ≥ ζ) ≤ 2e - 2ζ 2 n (e U -e L ) 2

Combining the above bounds results in:

Pr(|I(X, Z) n, θ -I(X, Z) ∞, θ | ≤ ξ + 2ζ) ≥ 1 -2e -2ξ 2 n (U -L) 2 -4e -2ζ 2 n (e U -e L ) 2

By solving ξ to minimize n according to Eq.8 we have:

Pr(|I(X, Z) n, θ -I(X, Z) ∞, θ | ≤ ) ≥ 1 -δ.

Theorem 1 also implies the following MI lower confidence interval under limited number of samples

Pr(I(X; Z) ≥ I(X, Z) n, θ -) ≥ 1 -δ

Compared to MINE, as per the example shown in Section 3, for M = 1 (i.e. L = -1 and U = 1), δ = 0.05, = 0.1, our estimator requires n = 10, 742 compared to MINE requiring n = 18, 756, 256 i.i.d validation samples to estimate a lower bound, which makes MI-based dependency analysis feasible for domains where data collection is prohibitively expensive, e.g. fMRI scans. In practice, sample complexity can be further optimized by optimizing hyperparameters U and L.

Note that unlike Eq.3, Theorem 1 bounds the closeness of the DEMINE estimate, I(X, Z) n, θ , not towards the MI lower bound sup θ∈Θ I MINE-f , but towards the MI lower bound I(X, Z) ∞, θ . Therefore, the sample complexity of DEMINE as in Eq.8 makes fair comparison with the sample complexity of MINE as in Eq.4. MINE's higher sample complexity stems from the need to bound the generalization error of T θ (X, Z) on unseen {(x, z)}. Existing generalization bounds are known to be overly loose, as over-parameterized neural networks have been shown to generalize well in classification and regression tasks (Zhang et al., 2016) . By using a learning-based formulation, DEMINE not only avoids the need to bound generalization error, but also allows further generalization improvements by learning θ through meta-learning.

In the following section, we present a meta-learning formulation, Meta-DEMINE, that learns θ for generalization given the same model class and training samples.

4.2. META-LEARNING:

Given training data {(x i , z i ) train , i = 1, . . . m}, Meta-DEMINE first generates MI estimation tasks each consisting of a meta-training split A and a meta-val split B through a novel task augmentation process. And then a parameter initialization θ init is then learned to maximize MI estimation performance on the generated tasks using initialization θ init as shown in Eq.9.

EQUATION

)

Here θ (t) = MetaTrain (x, z) A , θ (0) is the meta-training process of starting from an initialization θ (0) and applying Stochastic Gradient Descent (SGD)3 over t steps to learn θ where in every meta training iteration we have:

θ (t) ← θ (t-1) -γ∇L((x, z) A , θ (t-1) ).

Finally, θ is learned using the entire training set {(x i , z i ) train , i = 1, . . . , m} with θ init as initialization: θ = MetaTrain (x, z) train , θ init .

Task Augmentation: Meta-DEMINE adapts MAML (Finn et al., 2017a) for MI lower bound maximization. MAML has been shown to improve generalization performance in N -class K-shot image classification. MI estimation, however, does not come with predefined classes and tasks. A naive approach to produce tasks would be through cross-validation -partitioning training data into meta-training and meta-validation splits. However, merely using cross-validation tasks is prone to overfitting -a θ init , which memorizes all training samples would as a result have memorized all metavalidation splits. Instead, Meta-DEMINE generates tasks by augmenting the cross-validation tasks through task augmentation. Training samples are first split into meta-training and meta-validation splits, and then transformed using the same random invertible transformation to increase task diversity. Meta-DEMINE generates invertible transformation by sequentially composing the following functions:

Mirror : m(x) = (2n -1)x, n ∼ Bernoulli( 1 2 ), Permute : P (x) = n P d , Permute dimensions. Offset : O(x) = x + , ∼ U (-0.1, 0.1), Gamma : G(x) = sign(x) |x| γ , γ ∼ U (0.5, 2),

Since the MI between two random variables is invariant to invertible transformations on each variable, MetaTrain(•, •) is expected to arrive at the same MI lower bound estimation regardless of the transformation applied. At the same time, memorization is greatly suppressed, as the same pair (x, z) can have different log p(x,z) p(x)p(z) under different transformations. More sophisticated invertible transformations (affine, piece-wise linear) can also be added. Task augmentation is an orthogonal approach to data augmentation. Using image classification as an example, data augmentation generates variations of the image, translated, or rotated images assuming that they are valid examples of the class. Task augmentation on the other hand, does not make such an assumption. Task augmentation requires the initial parameters θ init to be capable of recognizing the same class in a world where all images are translated and/or rotated, with the assumption that the optimal initialization should easily adapt to both the upright world and the translated and/or rotated world. Optimization: Solving θ init using the meta-learning formulation Eq.9 poses a challenging optimization problem. The commonly used approach is back propagation through time (BPTT) which computes second order gradients and directly back propagates gradients from MetaTrain((x, z) A , θ (0) ) to θ init . BPTT is very effective for a small number of optimization steps, but is vulnerable to exploding gradients and is memory intensive. In addition to BPTT, we find that stochastic finite difference algorithms such as Evolution Strategies (ES) (Salimans et al., 2017) and Parameter-Exploring Policy Gradients (PEPG) (Sehnke et al., 2010) can sometimes improve optimization robustness. In practice, we switch betwen BPTT and PEPG depending on the number of meta-training iterations. Meta-DEMINE algorithm is specified in Algorithm 1.

5. EVALUATION ON SYNTHETIC DATASETS:

Dataset. We evaluate our approaches DEMINE and Meta-DEMINE against baselines and state-ofthe-art approaches on 3 synthetic datasets: 1D Gaussian, 20D Gaussian and sine wave. For 1D and 20D Gaussian datasets, following Belghazi et al. ( 2018), we define two k-dimensional multivariate Gaussian random variables X and Z which have component-wise correlation corr(X i , Z j ) = δ ij ρ, where ρ ∈ (-1, 1) and δ ij is Kronecker's delta. Mutual information I(X; Z) has a closed form solution I(X; Z) = -k ln(1 -ρ 2 ). For the sine wave dataset, we define two random variables X and Z, where X ∼ U(-1, 1), Z = sin(aX + π

2 ) + 0.05 , and ∼ N (0, 1). Estimating mutual information accurately given few pairs of (X, Z) requires the ability to extrapolate the sine wave given few examples. Ground truth MI for sine wave dataset is approximated by running the the KSG Estimator (Kraskov et al., 2004 ) on 1, 000, 000 samples. Implementation. We compare our estimators, DEMINE and Meta-DEMINE, against the KSG estimator (Kraskov et al., 2004) MI-KSG and MINE-f (Belghazi et al., 2018) . For both DEMINE and Meta-DEMINE, we study variance reduction mode, referred to as -vr, where hyperparameters are selected by optimizing 95% confident estimation mean (µ -2σ µ ) and statistical significance mode, referred to as -sig, where hyperparameters are selected by optimizing 95% confident MI Algorithm 1 Meta-DEMINE Input Data: {(x, z) train , (x, z) val } Parameters: batch B, Meta Learning Iterations N M , Task Augmentation Iterations N T , Optimization Iterations N O , Ratio r, Learning rate η, Meta Learning Rate η meta Output: MI, T θinit (X, Z), T θ (X, Z) 1: for i = 1 : N M do 2:

for j = 1 : N T do 3:

A = r × train, B = train -A 4: Split (x, z) train into (x, z) A and (x, z) B 5: Transformation R x for x, R x (•) = m(P(O(G(•)))) 6: Transformation R z for z, R z (•) = m(P(O(G(•)))) 7: θ (0) meta ← θ init 8: for k = 1 : N O do 9: Sample a batch of (x, z) B ∼ (x, z) A 10: Compute L (R x (x), R z (z)) B , θ (k) meta 11: Compute ∇ θ (k) meta L -gradient for θ meta 12:

Update θ meta using Adam Kingma & Ba (2014) with η 13: end for

14: Compute L meta (R x (x), R z (z)) B , θ (N O ) meta 15:

Compute ∇ θ0 L meta -gradient to θ init using BPTT lower bound (µ -). Samples (x, z) are split 50%-50% into (x, z) train and (x, z) val . We use a separable network architecture T θ (x, z) = M tanh(w cos f (x), g(z) + b) -t . f and g are MLP encoders that embed signals x and z into vector embeddings. Hyperparameters t ∈ [-1, 1] and M control upper and lower bounds T θ (x, z) ∈ [-M (1 + t), M (1 -t)]. Parameters w and b are learnable parameters. MLP design and optimization hyperparameters are selected using Bayesian hyperparameter optimization (Bergstra et al., 2013) described below.

Hyperparameter search on DEMINE-vr and DEMINE-sig was conducted using the hyperopt package4 . Seven hyperparameters were involved in hyperparameter search: 1) number of encoder layers Meta-learning hyperparameters are chosen as outer loop N M = 3, 000 iterations, task augmentation N T = 1 iterations, r = 0.8, η meta = η 3 , with task augmentation mode m(P (O(•))). N O was capped at 30 iterations for 1D and 20D Gaussian datasets due to memory limit. For the sine wave datasets with large N O , we used PEPG (Sehnke et al., 2010) rather than BPTT.

For MI-KSG, we use off-the-shelf implementation by Gao et al. (2017) with default number of nearest neighbors k = 3. MI-KSG does not provide any confidence interval. For MINE-f, we use the same network architecture same as DEMINE-vr. we implement both the original formulation which optimizes T θ on (x, z) till convergence (10k iters), as well as our own implementation MINE-f-ES with early stopping, where optimization is stopped after the same number of iterations as DEMINEvr to control overfitting.

Results:

. Figure 1 (a) shows MI estimation performance on 20D Gaussian datasets with varying ρ ∈ {0, 0.1, 0.2, 0.3, 0.4, 0.5} using N = 300 samples. Results are averaged over 5 runs to compare estimator bias, variance and confidence. Note that Meta-DEMINE-sig detects the highest p < 0.05 confidence MI, outperforming DEMINE-sig which is a close second. Both detect p < 0.05 statistically significant dependency starting ρ = 0.3, whereas estimations of all other approaches are low confidence. It shows that in contrary to common belief, estimating the variational lower bounds with high confidence can be challenging under limited data. MINE-f estimates MI > 3.0 and MINE-f-ES estimates positive MI when ρ = 0, both due to overfitting, despite MINE-f-ES having the lowest empirical bias. DEMINE variants have relatively high empirical bias but low variance due to tight upper and lower bound control, which provides a different angle to understand bias-variance trade off in MI estimation (Poole et al., 2018) .

Figure 1(b, c, d ) shows MI estimation performance on 1D, 20D Gaussian and sine wave datasets with fixed ρ = 0.8, 0.3 and a = 8π respectively, with varying N ∈ {30, 100, 300, 1000, 3000} number of samples. More samples asymptotically improves empirical bias across all estimators. As opposed to 1D Gaussian datasets which are well solved by N = 300 samples, higher-dimensional 20D Gaussian and higher-complexity sine wave datasets are much more challenging and are not solved using N = 3000 samples with a signal-agnostic MLP architecture. DEMINE-sig and Meta-DEMINE-sig detect p < 0.05 statistically significant dependency on not only 1D and 20D Gaussian datasets where x and z have non-zero correlation, but also on the sine wave datasets where correlation between x and z is 0. This means that DEMINE-sig and Meta-DEMINE-sig can be used for nonlinear dependency testing to complement linear correlation testing.

We study the effect of cross-validation meta-learning and task augmentation on 20D Gaussian with ρ = 0.3 and N = 300. )) prevent overfitting, but do not provide performance benefits, possibly because their complexity is insufficient or excessive for 20 adaptation steps. Further more, task augmentation with no adaptation (N O = 0) falls back to data augmentation, where samples from transformed distributions are directly used to learn T θ (x, z). Data augmentation with O(•) outperforms no augmentation, but is unable to outperform baseline and suffers from overfitting. It shows that task augmentation provides improvements orthogonal to data augmentation.

6. APPLICATION: FMRI INTER-SUBJECT CORRELATION (ISC) ANALYSIS:

Humans use language to effectively transmit brain representations among conspecifics. For example, after witnessing an event in the world, a speaker may use verbal communication to evoke neural representations reflecting that event in a listener's brain (Hasson et al., 2012) . The efficacy of this transmission, in terms of listener comprehension, is predicted by speaker-listener neural synchrony and synchrony among listeners (Stephens et al., 2010) . To date, most work has measured brainto-brain synchrony by locating statistically significant inter-subject correlation (ISC); quantified as the Pearson product-moment correlation coefficient between response time series for corresponding voxels or regions of interest (ROIs) across individuals (Hasson et al., 2004; Schippers et al., 2010; Silbert et al., 2014; Nastase et al., 2019) . Using DEMINE and Meta-DEMINE for statistical dependency testing, we can extend ISC analysis to capture nonlinear and higher-order interactions in continuous fMRI responses. Specifically, given synchronized fMRI response frames in two brain 42.8 28.0 32.8 35.9 0.637 16.5 7.9 11.6 12.0 0.035 Meta-DEMINE-vr 1vR 47.2 32.5 39.9 41.0 0.752 13.7 7.9 8.2 8.9 0.031 Abbreviations: P: Pieman; F: Forgot; Br: Bronx; Bk: Black, MI: Mutual Information.

*Note that all the results are averaging over the subjects.

regions X and Z across K subjects X i , Z i , i = 1, . . . , K as random variables. We model the conditional mutual information I(X i ; Z j |i = j) as the MI form of pair-wise ISC analysis. By definition, I(X i ; Z j |i = j) first computes MI between activations X i and Z j from subjects i and j respectively, and then average across pairs of subjects i = j. It can be lower bounded using Eq. 7 by learning a T θ (x, z) shared across all subject pairs.

Dataset. We study MI-based and correlation-based ISC on a fMRI story comprehension dataset by Nastase et al. ( 2019) with 40 participants listening to four spoken stories. Average story duration is 11 minutes. An fMRI frame with full brain coverage is captured at repetition time 1 TR =1.5 seconds with 2.5mm isotropic spatial resolution. We restricted our analysis to subsets of voxels defined using independent data from previous studies: functionally-defined masks of high ISC voxels (ISC; 3,800 voxels) and dorsal Default-Mode Network voxels (dDMN; 3,940 voxels) from Simony et al. (2016) Results. We first examine, for the fine grained HCM-MMP1 brain regions, which have p < 0.05 statistically significant MI and Pearson's correlation. Table 1 shows the result. Overall, more regions have statistically significant correlation than dependency. This is expected because correlation requires less data to detect. But Meta-DEMINE is able to find 6 brain regions that have statistically significant dependency but lacks significant correlation. This shows that MI analysis can be used to complement correlation-based ISC analysis.

By considering temporal ISC over time, fMRI signals can be modeled with improved accuracy. In Table 2 we apply DEMINE and Meta-DEMINE with L = 10TRs (15s) sliding windows as random variables to study amount of information that can be extracted from ISC and dDMN masks. We use between-subject time-segment classification (BSC) for evaluation (Haxby et al., 2011; Guntupalli et al., 2016) . Each fMRI scan is divided into K non-overlapping L = 10 TRs time segments. The BSC task is one versus rest retrieval: retrieve the corresponding time segment z of an individual given a group of time segments x excluding that individual, measured by top-1 accuracy. For retrieval score, T θ (X, Z) is used for DEMINE and Meta-DEMINE and ρ(X, Z) is used for Pearson's correlation as a simple baseline. With CNN as encoder, DEMINE and Meta-DEMINE model the signal better and achieve higher accuracy. Also. Meta-DEMINE is able to extract 0.75 nats of MI from the ISC mask over 10 TRs or 15s, which could potentially be improved by more samples.

7. CONCLUSION:

We illustrated that a predictive view of the MI lower bounds coupled with meta-learning results in data-efficient variational MI estimators, DEMINE and Meta-DEMINE, that are capable of performing statistical test of dependency. We also showed that our proposed task augmentation reduces overfitting and improves generalization in meta-learning. We successfully applied MI estimation to real world, data-scarce, fMRI datasets. Our results suggest a greater avenue of using neural networks and meta-learning to improve MI analysis and applying neural network-based information theory tools to enhance the analysis of information processing in the brain. Model-agnostic, high-confidence, MI lower bound estimation approaches -including MINE, DEMINE and Meta-DEMINE-are limited to estimating small MI lower bounds up to O(log N ) as pointed out in (McAllester & Statos, 2018) , where N is the number of samples. In real fMRI datasets, however, strong dependency is rare and existing MI estimation tools are limited more by their ability to accurately characterize the dependency. Nevertheless, when quantitatively measuring strong dependency, cross-entropy (McAllester & Statos, 2018) or model-based quantities, alternatives to MI, such as correlation or CCA, may be measured with high confidence.

A APPENDIX:

A.1 THE DEMINE ALGORITHM (Glorot & Bengio, 2010) 2: for i = 1 :

Algorithm 2 DEMINE Input Data: {(x, z) train , (x, z) val } Parameters: Batch B, Iterations N O , Learning rate η Output: MI, T θ (X, Z) 1: θ (0) ← Xavier Initialization

N O do 3: Sample a batch of (x i , z i ) B ∼ (x, z) train 4: Compute L (x i , z i ) B , θ (i-1) 5: Compute ∇ (i) θ L -gradient for θ 6:

Update θ (i) using Adam (Kingma & Ba, 2014) with η 7: end for

8: MI = I(X, Z) n,θ (N O ) 9: return MI, θ (N O )

A.2 ADDITIONAL DETAILS OF THE FMRI DATASET:

The dataset we used contains 40 participants (mean age = 23.3 years, standard deviation = 8.9, range: 1853; 27 female) recruited to listen to four spoken stories 56 . The stories were renditions of "Pie Man" and "Running from the Bronx" by Jim OGrady (O'Grady, 2018b;a), "The Man Who Forgot Ray Bradbury" by Neil Gaiman (Gaiman, 2018) , and "I Knew You Were Black" by Carol Daniel (Daniel, 2018) ; story durations were 7, 9, 14, and 13 minutes, respectively. After scanning, participants completed a questionnaire comprising 25-30 questions per story intended to measure narrative comprehension. The questionnaires included multiple choice, True/False, and fill-in-theblank questions, as well as four additional subjective ratings per story. Functional and structural images were acquired using a 3T Siemens Prisma with a 64-channel head coil. Briefly, functional images were acquired in an interleaved fashion using gradient-echo echo-planar imaging with a multiband acceleration factor of 3 (TR/TE = 1500/31 ms where TE stands for "echo time", resolution = 2.5 mm isotropic voxels, full brain coverage). All fMRI data were formatted according to the Brain Imaging Data Structure (BIDS) standard (Gorgolewski et al., 2016) and preprocessed using the fMRIPrep library (Esteban et al., 2018) . Functional data were corrected for slice timing, head motion, and susceptibility distortion, and normalized to MNI space using nonlinear registration. Nuisance variables comprising head motion parameters, framewise displacement, linear and quadratic trends, sine/cosine bases for high-pass filtering (0.007 Hz), and six principal component time series from cerebrospinal fluid (CSF) and white matter (WM) were regressed out of the signal using the Analysis of Functional NeuroImages (AFNI) software suite (Cox, 1996) .

The fMRI data comprise X ∈ R Vi×T for each subject, where V i represents the flattened and masked voxel space and T represents the number of samples (in TRs) during auditory stimulus presentation.

Additional Details on Dataset Collection Functional and structural images were acquired using a 3T Siemens Magnetom Prisma with a 64-channel head coil. Functional, blood-oxygenation-leveldependent (BOLD) images were acquired in an interleaved fashion using gradient-echo echo-planar imaging with pre-scan normalization, fat suppression, a multiband acceleration factor of 3, and no in-plane acceleration: TR/TE = 1500/31 ms, flip angle = 67 • , bandwidth = 2480 hz per pixel, resolution = 2.5 mm 3 isotropic voxels, matrix size = 96 x 96, Field of view (FoV) = 240 x 240 mm, 48 axial slices with roughly full brain coverage and no gap, anteriorposterior phase encoding. At the beginning of each scanning session, a T1-weighted structural scan (where T1 stands for "longitudinal relaxation time"), was acquired using a high-resolution single-shot Magnetization-Prepared 180 degrees radio-frequency pulses and RApid Gradient-Echo (MPRAGE) sequence with an in-plane acceleration factor of 2 using GeneRalized Autocalibrating Partial Parallel Acquisition (GRAPPA): TR/TE/TI = 2530/3.3/1100 ms where TI stands for inversion time, flip angle = 7 • , resolution = 1.0 x 1.0 x 1.0 mm voxels, matrix size = 256 x 256, FoV = 256 x 256 x 176 mm, 176 sagittal slices, ascending acquisition, anteriorposterior phase encoding, no fat suppression, 5 min 53 s total acquisition time. At the end of each scanning session a T2-weighted (where T2 stands for "transverse relaxation time") structural scan was acquired using the same acquisition parameters and geometry as the T1-weighted structural image: TR/TE = 3200/428 ms, 4 minutes 40 seconds total acquisition time. A field map was acquired at the beginning of each scanning session, but was not used in subsequent analyses.

Additional Details on Dataset Preprocessing Preprocessing was performed using the fMRIPrep library 7 Esteban et al. (2018), a Nipype library8 (Gorgolewski et al., 2011) based tool. T1-weighted images were corrected for intensity non-uniformity using the N4 bias field correction algorithm (Tustison et al., 2010) and skull-stripped using Advanced Normalization Tools (ANTs) (Avants et al., 2008) . Nonlinear spatial normalization to the International Consortium for Brain Mapping (ICBM) 152 Nonlinear Asymmetrical template version 2009c (Fonov et al., 2009) was performed using ANTs. Brain tissue segmentation cerebrospinal fluid, white matter, and gray matter was was performed using FSL library's9 FAST tool Zhang et al. (2001) . Functional images were slice timing corrected using AFNI software's 3dTshift (Cox, 1996) and corrected for head motion using FSL library's MCFLIRT tool (Jenkinson et al., 2002) . "Fieldmap-less" distortion correction was performed by co-registering each subject's functional image to that subject's intensity-inverted T1-weighted image (Wang et al., 2017) constrained with an average field map template (Treiber et al., 2016) . This was followed by co-registration to the corresponding T1-weighted image using FreeSurfer software's10 boundary-based registration (Greve & Fischl, 2009) with 9 degrees of freedom. Motion correcting transformations, field distortion correcting warp, BOLD-to-T1 transformation and T1-to-template (MNI) warp were concatenated and applied in a single step with Lanczos interpolation using ANTs. Physiological noise regressors were extracted applying "a Component Based Noise Correction Method" aCompCor (Behzadi et al., 2007) . Six principal component time series were calculated within the intersection of the subcortical mask and the union of CSF and WM masks calculated in T1w (T1 weighted) space, after their projection to the native space of each functional run. Framewise displacement (Power et al., 2014) was calculated for each functional run. Functional images were downsampled to 3 mm resolution. Nuisance variables comprising six head motion parameters (and their derivatives), framewise displacement, linear and quadratic trends, sine/cosine bases for high-pass filtering (0.007 Hz cutoff), and six principal component time series from an anatomically-defined mask of cerebrospinal fluid and white matter were regressed out of the signal using AFNI's 3dTproject (Cox, 1996) . Functional response time series were z-scored for each voxel.

A.3 COMPARISON BETWEEN HSIC AND DEMINE:

We first review the Hilbert-Schmidt independence criterion (HSIC), a widely-studied correlationbased independence criterion and discuss its connections with the MINE family of mutual information lower bound methods, and then study DEMINE and a spectral HSIC implementation on the synthetic datasets.

The HSIC approach (Gretton et al., 2005b; a) is based on a necessary and sufficient condition of independence: two random variables X and Z are independent if and only if for all bounded or positive functions f and g R

d → R, E XZ f (X)g(Z) -E X f (X)E Z g(Z) = 0, or equivalently E XZ (f (X)-E X f (X))(g(Z)-E Z g(Z)) = 0.

A proof can be constructed by showing equivalence to the definition of independence, P (X, Z) = P (X)P (Z).

To construct an independence test, existing approaches (Gretton et al., 2005b; a) use Reproducing Kernel Hilbert Spaces (RKHS) for f and g, a function space that not only covers all functions between [0, 1], but also allows computationally efficient estimation or bounding of

COCO(X, Z) = sup f,g E XZ f (X)g(Z) -E X f (X)E Z g(Z)

given samples, and test COCO(X, Z) = 0. Confidence intervals are derived through McDiarmid's inequality, or using closed-form distributions to approximate the test statistics to a certain order of moments, and compute the confidence interval from the closed-form distribution.

The COCO(X, Z) used by HSIC estimators bears great resemblance to the MINE family of mutual information estimators. In fact, it can be shown that

EQUATION

It means that within a family of decomposable functions where T θ (X, Z) = f (X)g(Z), COCO(X,Z) is an upperbound to the MINE estimates. In addition, the equivalence of COCO(X,Z) = 0 and I(X, Z) = 0 seems to suggest a form of mutual information bound. On the other hand, MINE allows the use of non-decomposable T θ (X, Z). Existing results on MINE (Poole et al., 2018) seem to suggest that a non-decomposable T θ (X, Z) gives superior empirical mutual information estimation performance over a decomposable T θ (X, Z). The necessity of non-decomposable T θ (X, Z) designs and mutual information lower bounds under decomposable designs of T θ (X, Z) may be subjects of further research.

Similar to the MINE estimators, HSIC-based estimators tend to have loose confidence intervals due to the need to bound generalization error of kernels f and g on unseen data points. We expect a cross-validation-based approach like DEMINE to also improve the performance of the HSIC-based estimators.

Comparison between DEMINE and HSIC on synthetic benchmarks. We compare Canonical Correlation Analysis (CCA), DEMINE, DEMINE-meta and HSIC for independent testing on our 4 synthetic Gaussian and sine wave benchmarks presented in Section 5. Results for a single random seed is reported for a compact presentation, but we have ran experiments using multiple random seeds and find the result of a single random seed representative enough.

For CCA, we compute p-value using the χ 2 test. For HSIC, we report p-value using a publicly available implementation for a spectral HSIC test (Zhang et al., 2018) 11 . The default kernel is used. Hyperparameters are set to recommended setting when available. For DEMINE and DEMINE-meta, the setup is identical to Section 5. A 2-sided 95% confidence interval is reported, but showing only the lower side.

Experiment results are compiled in Table 3 . Statistically significant dependence detections with p < 0.05 are bolded. Results show that spectral HSIC requires less data to test dependency for the simple Gaussians dataset. But on the more challenging sine wave dataset, DEMINE-sig and DEMINEmeta-sig perform better. Overall, we find DEMINE more complementary to linear correlations for dependency testing on complex signals. Note that Gaussian kernels are used for spectral HSIC. More complex kernels have potential to improve results.

B SANITY CHECK ON STATISTICAL DEPENDENCY TESTING:

We performed sanity check of our approach, as well as several statistical dependency testing implementations that we compare against. We run different statistical dependency testing implementations on our 1D Gaussian ρ = 0.0, N = 30 samples dataset where X and Z are independent. A large number of runs with different random seeds are performed. False positive rate of p < 0.05 statistical 

:

We follow the same notations inBelghazi et al. (2018). We encourage the review of Belghazi et al. (2018);Poole et al. (2018) for a detailed understanding of IMINE, IEB1, and IEB.

MINE lower bound can also be interpreted in the predictive way, but will result in a higher sample complexity than MINE-f lower bound. We choose MINE-f in favor of a lower sample complexity over bound tightness.

In practice, the Adam optimizer(Kingma & Ba, 2014) is used for faster optimization. The Adam optimizer uses first and second order momentums of the gradient to speed up optimization. Illustrating SGD for simplicity.

Hyperopt package: https://github.com/hyperopt/hyperopt.

Two of the stories were told by a professional storyteller undergoing an fMRI scan; however, fMRI data for the speaker were not analyzed for the present work due to the head motion induced by speech production.

The study was conducted in compliance with the Institutional Review Board of the University

https://github.com/poldracklab/fmriprep

https://github.com/nipy/nipype

https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL

https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferWiki

https://github.com/oxmlcs/kerpy. We also experimented with classic HSIC with gamma approximationGretton et al. (2005a) https://github.com/amber0309/HSIC and a block HSIC implementation(Zhang et al., 2018) from https://github.com/oxmlcs/kerpy, but find that they both report significantly more than 5% false positives for independent 1D and 20D gaussians ρ = 0 at N = {30, 100} across 100,000-1,000,000 random seeds, indicating errors in confidence interval

calculations.12 This is a false positive case for CCA, because for this sine wave data ground truth correlation is 0.

Figures and Tables:

θ init using Adam Kingma & Ba (2014) with η meta 18: end for 19: θ (0) ← θ init 20: for i = 1 : N O do 21: Sample a batch of (x, z) B ∼ (x, z) train 22: Compute L (x, z) B , θ (i) 23: Compute gradient ∇ θ L 24: Update θ using Adam with η 25: end for 26: Compute MI = L (x, z) val , θ (N O ) 27: return MI, θ init , θ (N O )

[1, 5], 2) encoder hidden size[8, 256], 3) learning rate η [10 -4 , 3 × 10 -1 ] in log scale, 4) number of optimization iterations N O[5, 200]  (sine wave [5, 5000]) in log scale, 5) batch size B [256, 1024], 6) M , [10 -3 , 5] in log scale, 7) t, [-1, 1]. Mean µ and sample standard deviation σ of MI estiamte computed over 3-fold cross-validation on (x, z) train . DEMINE-vr maximizes two sigma low µ-2σ µ where σ µ = 1 √ 3 σ due to 3-fold cross-validation. DEMINE-sig maximizes statistical significance µwhere is two-sided 95% confidence interval of MI. Meta-DEMINE-vr and Meta-DEMINEsig subsequently reuse these hyperparameters as DEMINE-vr and DEMINE-sig.

Figure 2 plots performance of Meta-DEMINE-vr over N M = 3000 meta iterations under combinations of task augmentations modes and number of adaptation iterations N O ∈ {0, 20}. Overall, task augmentation modes which involve axis flipping m(•) and permutation P (•) are the most successful. With N O = 20 steps of adaptation, task augmentation modes P (•), m(P (•)) and m(P (O(•))) prevent overfitting and improves performance. The performance improvements of task augmentation is not simply from change in batch size, learning rate or number of optimization iterations, because meta-learning without task augmentation for both N O = 0 and 20 could not outperform baseline. Meta-learning without task augmentation and with task augmentation but using only O(•) or G(•) result in overfitting. Task augmentation with m(•) or m(P (O(G(•))

Figure 1: Comparing MI Estimation performance of DEMINE and Meta-DEMINE with the KSG estimator Kraskov et al. (2004) and MINE-f Belghazi et al. (2018) on different datasets using varying number of samples. The bars show estimator mean and standard deviation averaged over 5 runs with different seeds. The error bars show 95% confidence interval (not available for MI-KSG). The statistical significance focused variants DEMINE-sig and Meta-DEMINE-sig achieves the highest 95% confident MI estimation. Meta-DEMINE improves over DEMINE most of the time.

(a) Meta-DEMINE-vr NO = 0. (b) Meta-DEMINE-vr NO = 10. (c) Meta-DEMINE-vr NO = 20.

Figure 2: To study the effect of task augmentation and number of adaptation steps, we run Meta-DEMINE-vr with different task augmentation modes and vary number of adaptation iterations N O ∈ {0, 10, 20} on Gaussian 20D, ρ = 0.3 dataset. Combinations of permutation and mirroring operations are effective in reducing overfitting and improving performance.Table 1: Number of HCP-MMP1 regions with significant correlation (r) and MI (DEMINE, Meta-DEMINE) during listening.

Number

Table content: <table><tr><td/><td colspan="3">of HCP-MMP1</td></tr><tr><td colspan="4">regions with significant correla-</td></tr><tr><td colspan="4">tion (r) and MI (DEMINE, Meta-</td></tr><tr><td colspan="3">DEMINE) during listening.</td><td/></tr><tr><td>No. shared</td><td colspan="2">r DEMINE</td><td>Meta</td></tr><tr><td/><td/><td>-sig</td><td>-DEMINE</td></tr><tr><td/><td/><td/><td>-sig</td></tr><tr><td>r</td><td>37</td><td>24</td><td>23</td></tr><tr><td>DEMINE-sig</td><td>24</td><td>28</td><td>26</td></tr><tr><td colspan="2">Meta-DEMINE-sig 23</td><td>26</td><td>29</td></tr></table>

Segment classification accuracy for NeuralMI versusPearson's correlation in 1-vs-rest*.

Table content: <table><tr><td>Classification</td><td/><td/><td colspan="2">ISC Mask</td><td/><td/><td colspan="3">dDMN Mask</td><td/></tr><tr><td>Accuracy (%)</td><td>P</td><td>F</td><td>Br</td><td>Bk</td><td>MI</td><td>P</td><td>F</td><td>Br</td><td>Bk</td><td>MI</td></tr><tr><td>Chance</td><td colspan="10">3.7 1.8 2.6 1.9 N/A 3.7 1.8 2.6 1.9 N/A</td></tr><tr><td>Pearson's r 1vR</td><td colspan="10">35.0 20.4 25.8 31.5 N/A 14.8 6.4 11.8 9.9 N/A</td></tr><tr><td>DEMINE-vr 1vR</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

as well as 180 HCP-MMP1 multimodal cortex parcels fromGlasser et al. (2016). All masks were defined in MNI space. Implementation. We compare MI-based ISC using DEMINE and Meta-DEMINE with correlationbased ISC using Pearson's correlation. DEMINE and Meta-DEMINE setup follows Section 5. The fMRI data were partitioned by subject into a train set of 20 subjects and a validation set of 20 different subjects. Residual 1D CNN is used instead of MLP as the encoder for studying temporal dependency. For Pearson's correlation, high-dimensional signals are reshaped to 1D for correlation analysis. Effective sample size for confidence interval calculation is the number of unique nonoverlapping fMRI samples.

Table content: <table/>

Statistical dependency testing comparison between CCA, spectral HSIC and DEMINE algorithms on synthetic datasets. Statistically significant dependency with p < 0.05 detections are bolded. See text for detailed experimental setups.

Table content: <table><tr><td>Problem</td><td>Samples</td><td>CCA</td><td>Spectral HSIC</td><td>DEMINE-sig</td><td>DEMINE-meta-sig</td></tr><tr><td>1D Gaussian ρ = 0.8</td><td>30</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.006-0.017</td><td>I ≥ 0.006-0.017</td></tr><tr><td>1D Gaussian ρ = 0.8</td><td>100</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.006-0.009</td><td>I ≥ 0.006-0.009</td></tr><tr><td>1D Gaussian ρ = 0.8</td><td>300</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.143-0.132</td><td>I ≥ 0.146-0.132</td></tr><tr><td>1D Gaussian ρ = 0.8</td><td>1000</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.278-0.168</td><td>I ≥ 0.292-0.168</td></tr><tr><td>1D Gaussian ρ = 0.8</td><td>3000</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.365-0.146</td><td>I ≥ 0.344-0.146</td></tr><tr><td>20D Gaussian ρ = 0.3</td><td>30</td><td>p=0.015</td><td>p=0.019</td><td>I ≥ -0.003-0.017</td><td>I ≥ 0.001-0.017</td></tr><tr><td>20D Gaussian ρ = 0.3</td><td>100</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.000-0.009</td><td>I ≥ 0.001-0.009</td></tr><tr><td>20D Gaussian ρ = 0.3</td><td>300</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.005-0.005</td><td>I ≥ 0.007-0.005</td></tr><tr><td>20D Gaussian ρ = 0.3</td><td>1000</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.322-0.170</td><td>I ≥ 0.376-0.170</td></tr><tr><td>20D Gaussian ρ = 0.3</td><td>3000</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.632-0.253</td><td>I ≥ 0.689-0.253</td></tr><tr><td>20D Gaussian ρ = 0.0</td><td>300</td><td>p=0.624</td><td>p=0.498</td><td>I ≥ 0.000-0.005</td><td>I ≥ 0.000-0.005</td></tr><tr><td>20D Gaussian ρ = 0.1</td><td>300</td><td>p=0.000</td><td>p=0.014</td><td>I ≥ 0.000-0.005</td><td>I ≥ 0.000-0.005</td></tr><tr><td>20D Gaussian ρ = 0.2</td><td>300</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.002-0.005</td><td>I ≥ 0.003-0.005</td></tr><tr><td>20D Gaussian ρ = 0.3</td><td>300</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.005-0.005</td><td>I ≥ 0.007-0.005</td></tr><tr><td>20D Gaussian ρ = 0.4</td><td>300</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.191-0.140</td><td>I ≥ 0.260-0.140</td></tr><tr><td>20D Gaussian ρ = 0.5</td><td>300</td><td>p=0.000</td><td>p=0.001</td><td>I ≥ 0.621-0.388</td><td>I ≥ 0.815-0.388</td></tr><tr><td>Sine wave a = 8π</td><td>30</td><td>p=0.826</td><td>p=0.856</td><td>I ≥ 0.000-0.005</td><td>I ≥ -0.002-0.017</td></tr><tr><td>Sine wave a = 8π</td><td>100</td><td>p=0.962</td><td>p=0.913</td><td>I ≥ 0.006-0.009</td><td>I ≥ 0.006-0.009</td></tr><tr><td>Sine wave a = 8π</td><td>300</td><td>p=0.093</td><td>p=0.498</td><td>I ≥ 0.448-0.306</td><td>I ≥ 0.511-0.306</td></tr><tr><td>Sine wave a = 8π</td><td>1000</td><td>p=0.047 12</td><td>p=0.111</td><td>I ≥ 1.351-0.450</td><td>I ≥ 1.316-0.450</td></tr><tr><td>Sine wave a = 8π</td><td>3000</td><td>p=0.166</td><td>p=0.094</td><td>I ≥ 1.711-0.442</td><td>I ≥ 1.756-0.442</td></tr></table>

Sanity check of several statistical dependency testing implementations. See text for details.

Table content: <table/>

Peer review (iclr2020):
Second update: We thank the authors for updating their paper. The work is now improving, and is on the right track for publication at a future conference. There are a few comments on the new results, and suggestions for further improvement: * The issue of possible false positives due to sample dependence in fMRI data has again been ignored in the rebuttal. Without careful treatment of this effect, these results are vacuous. * There is still no comparison with competing nonparametric tests on the fMRI data. * the results linking the COCO and MINE estimators are interesting. Some statements don*t make sense, however, eg. *HSIC-based estimators tend to have loose confidence intervals due to the need to bound generalization error of kernels f and g on unseen data points.* First, f and g are functions, not kernels. Second, testing with HSIC or COCO does not require generalisation to unseen data points: this is why testing is an easier problem than regression. For HSIC testing, I am surprised to read in the footnote that the Gamma approximation report significantly more than 5% errors, especially given the Table 4 results that show the correct level. In any case, I very strongly suggest using a permutation approach to obtain the test threshold for HSIC, which is by far the most robust and reliable method. The Gamma approximation has no statistical guarantees, as stated explicitly in the HSIC testing paper of NeurIPS 2007. Permutation gives a guarantee of the correct level. See the NeurIPS paper for details. Once you have verified the correct false positive rate for the permutation threshold, then you can compute the p-value on the alternative. While the comparison with HSIC is a helpful one, it is also required to compare with the competing method closest to yours, i.e. the Berrett and Samworth test. In addition, HSIC is a non-adaptive test, but your test is adaptive, so a fairer comparison would be to a modern adaptive test such as *An Adaptive Test of Independence with Analytic Kernel Embeddings.* * The false positive rate in the sanity check is far below the design level of 0.05. This is as I expected, given the use of the Hoeffding bound. This should be stated clearly in the main text, and not disclosed in the final sentence of the final page of the appendix. ==================== Update: thank you for your rebuttal. *the necessary and sufficient condition of dependency, is more general and is complementary to other techniques that make stronger assumptions about the data. * The alternative tests listed are nonparametric. That is to say, unlike Pearson correlation, they do not make specific parametric assumptions on the nature of the dependence. Rather they make generic smoothness assumptions (your test also makes such assumptions by the choice of neural network architecture). Thus, comparison with the prior work in Statistics and Machine Learning is relevant, since these tests have the same aims and scope as your tests. The cited Berrett and Samworth MI test uses a permutation approach to obtaining the test threshold, not an asymptotic approach (see the results of Section 4 of that paper). Several of the other cited tests also use a permutation approach for the test threhsold. These tests are therefore relevant prior work. * If things change very little from one second to the next, the signals could be very similar and may not really be, intuitively, independent samples and may bias result of the study. However, which independence assumptions to use is not in scope for our paper, because our fMRI study is trying to show that dependency testing works * In the work cited by Chwialkowski et al, failure to account for the dependence between samples results in excessive false positives. This is because, for dependent data, the effective sample size is reduced, and the tests must be made more conservative to correct for this effect. It is therefore the case that the fMRI results may be false positives. Re level:*This proof could be experimentally verified ...* This should be verified. In particular, Hoeffding can be very loose in practice, which is likely to be observed in experiments. ====== The authors propose a procedure for improving neural mutual information estimates, via a combination of data augmentation and cross validation. They then use these estimates in hypothesis testing, on low dimensional toy datasets and on high dimensional real-world fMRI data. The improved training procedures for MI estimation are of interest, however the hypothesis testing parts of the paper could still be improved. In hypothesis testing, it is important to verify that the test has the correct level (false positive rate). This is all the more essential when the estimate has required optimisation over parameters. It is not clear from the presentation that this has been confirmed. There are a number of prior approaches to testing for multivariate statistical dependence in the machine learning and statistics literature (including a 2017 paper which uses mutual information). A small selection is given below, although a literature search will reveal many more papers. In the absence of citation or comparison with any of the prior work on multivariate statistical dependence testing, the current submission is not suitable for publication. In statistics: --------------- https://arxiv.org/abs/1711.06642 Nonparametric independence testing via mutual information Thomas B. Berrett, Richard J. Samworth 2017 Measuring and testing dependence by correlation of distances Gábor J. Székely, Maria L. Rizzo, and Nail K. Bakirov Ann. Statist. Volume 35, Number 6 (2007), 2769-2794. Large-scale kernel methods for independence testing Qinyi ZhangEmail Sarah Filippi, Arthur Gretton, Dino Sejdinovic Statistics and Computing January 2018, Volume 28, Issue 1, pp 113–130| Cite as In machine learning: --------------------- Multivariate tests of association based on univariate tests Heller, Ruth and Heller, Yair Advances in Neural Information Processing Systems 29 2016 A Kernel Statistical Test of Independence Gretton, Arthur and Fukumizu, Kenji and Choon H. Teo and Song, Le and Sch*{o}lkopf, Bernhard and Alex J. Smola Advances in Neural Information Processing Systems 20 2008 http://papers.nips.cc/paper/3201-a-kernel-statistical-test-of-independence.pdf http://proceedings.mlr.press/v70/jitkrittum17a/jitkrittum17a.pdf An Adaptive Test of Independence with Analytic Kernel Embeddings Wittawat Jitkrittum, Zoltán Szabó, Arthur Gretton ; ICML 2017, PMLR 70:1742-1751 Time dependence ---------------- It is also the case that if the variables have time dependence, then appropriate corrections must be made for the test threshold, to avoid excessive false positives. Does the fMRI data exhibit time dependence? For the case of multivariate statistical dependence testing, such corrections are described e.g. in: https://papers.nips.cc/paper/5452-a-wild-bootstrap-for-degenerate-kernel-tests.pdf A Wild Bootstrap for Degenerate Kernel Tests Kacper Chwialkowski, Dino Sejdinovic, Arthur Gretton NeurIPS 2014