Nonparanormal Information Estimation
Shashank Singh Barnabás Pøczos


[
Nonparanormal Information Estimation





Shashank Singhsss1@andrew.cmu.edu
Carnegie Mellon University,
             5000 Forbes Ave., Pittsburgh, PA 15213 USA
Barnabás Póczosbapoczos@cs.cmu.edu
Carnegie Mellon University,
             5000 Forbes Ave., Pittsburgh, PA 15213 USA



mutual information,entropy,nonparanormal,Gaussian copula,Spearman,Kendall
0.3in]


We study the problem of using i.i.d. samples from an unknown multivariate probability distribution p to estimate the mutual information of p. This problem has recently received attention in two settings: (1) where p is assumed to be Gaussian and (2) where p is assumed only to lie in a large nonparametric smoothness class. Estimators proposed for the Gaussian case converge in high dimensions when the Gaussian assumption holds, but are brittle, failing dramatically when p is not Gaussian. Estimators proposed for the nonparametric case fail to converge with realistic sample sizes except in very low dimensions. As a result, there is a lack of robust mutual information estimators for many realistic data. To address this, we propose estimators for mutual information when p is assumed to be a nonparanormal (a.k.a., Gaussian copula) model, a semiparametric compromise between Gaussian and nonparametric extremes. Using theoretical bounds and experiments, we show these estimators strike a practical balance between robustness and scaling with dimensionality.





§ INTRODUCTION



This paper is concerned with the problem of estimating entropy or mutual information of an unknown probability density p over ^D, given n i.i.d. samples from p. Entropy and mutual information are fundamental information theoretic quantities, and consistent estimators for these quantities have a host of applications within machine learning, statistics, and signal processing. For example, entropy estimators have been used for goodness-of-fit testing <cit.>,  parameter estimation in semi-parametric models <cit.>, texture classification and image registration <cit.>,  change point detection <cit.>, and anomaly detection in networks <cit.>. Mutual information is a popular nonparametric measure of dependence, whose estimators have been used in feature selection <cit.>, clustering <cit.>, learning graphical models <cit.>, fMRI data processing <cit.>, prediction of protein structures <cit.>, boosting and facial expression recognition <cit.>, and fitting deep nonlinear models <cit.>. Estimators for both entropy and mutual information have been used in independent component and subspace analysis <cit.>.

Motivated by these and other applications, several very recent lines of work (discussed in [sec:related_work]Section <ref>) have studied information estimation,[We will collectively call the closely related problems of entropy and mutual information estimation information estimation.] focusing largely on two settings:


  * Gaussian Setting: If p is known to be Gaussian, there exist information estimators with mean squared error (MSE) at most 2log(1 - D/n) and an (almost matching) minimax lower bound of 2D/n <cit.>.

  * Nonparametric Setting: If p is assumed to lie in a nonparametric smoothness class, such an s-order[Here, s encodes the degree of smoothness, roughly corresponding to the number of continuous derivatives of p.] Hólder or Sobolev class, then the minimax MSE is of asymptotic order ≍max{ n, n^-8s/4s + D} <cit.>.


In the Gaussian setting, consistent estimation is tractable even in the high-dimensional case where D increases fairly quickly with n, as long as D/n → 0. However, optimal estimators for the Gaussian setting rely heavily on the assumption of joint Gaussianity, and their performance can degrade quickly when the data deviate from Gaussian. Especially in high dimensions, it is unlikely that data are jointly Gaussian, making these estimators brittle in practice. In the nonparametric setting, the theoretical convergence rate decays exponentially with D, and, it has been found empirically that information estimators for this setting fail to converge at realistic sample sizes in all but very low dimensions. Also, most nonparametric estimators are sensitive to tuning of bandwidth parameters, which is challenging for information estimation, since no empirical error estimate is available for cross-validation.

Given these factors, though the Gaussian and nonparametric cases are fairly well understood in theory, there remains a lack of practical information estimators for the common case where data are neither exactly Gaussian nor very low dimensional. The main goal of this paper is to fill the gap between these two extreme settings by studying information estimation in a semiparametric compromise between the two, known as the “nonparanormal” (a.k.a. “Gaussian copula”) model (see [def:nonparanormal]Definition <ref>). The nonparanormal model, analogous to the additive model popular in regression <cit.>, limits complexity of interactions among variables but makes minimal assumptions on the marginal distribution of each variable. The result scales better with dimension than nonparametric models, while being more robust than Gaussian models.

Paper Organization: [sec:problem_statement]Section <ref> gives definitions and notation to formalize the nonparanormal information estimation problem. [sec:related_work]Section <ref> discusses the history of the nonparanormal model and prior work on information estimation, motivating our contributions. [sec:estimators]Section <ref> proposes three estimators, while [sec:main_results]Section <ref> presents our theoretical error bounds, proven in the Appendix. 
[sec:empirical]Section <ref> provides simulation results. While most of the paper discusses mutual information estimation, [sec:entropy]Section <ref> discusses additional considerations arising in entropy estimation.  [sec:conc_and_future]Section <ref> presents some concluding thoughts and avenues for future work.



§ PROBLEM STATEMENT AND NOTATION


There are a number of distinct generalizations of mutual information to more than two variables. The definition we consider is simply the difference between the sum of marginal entropies and the joint entropy:

(Multivariate mutual information)
Let X_1,…,X_D be -valued random variables with a joint probability density
p : ^D → [0, ∞) and marginal densities p_1,...,p_D : → [0, ∞). The multivariate mutual information
I(X) of X = (X_1,…,X_D) is defined by

    I(X)
        := _X ∼ p[
            log(
                 p(X)/∏_j = 1^D p_j(X_j))
          ] 
        = ∑_j = 1^D H(X_j) - H(X),

where H(X) = -_X ∼ p [log p(X)] denotes entropy of X.


This notion of multivariate mutual information, originally due to <cit.> (who called it “total correlation”) measures total dependency, or redundancy, within a set of D random variables. It has also been called the “multivariate constraint” <cit.> and “multi-information” <cit.>. Many related information theoretic quantities can be expressed in terms of I(X), and can thus be estimated using estimators of I(X). Examples include pairwise mutual information I(X,Y) = I((X,Y)) - I(X) - I(Y),
which measures dependence between (potentially multivariate) random variables X and Y, conditional mutual information

    I(X|Z) = I((X,Z)) - ∑_j = 1^D I((X_j,Z)),

which is useful for characterizing how much dependence within X can be explained by a latent variable Z <cit.>, and transfer entropy (a.k.a. directed information) T_X → Y, which measures predictive power of one time series X on the future of another time series Y. I(X) is also related to entropy via Eq. (<ref>), but, unlike the above quantities, this relationship depends on the marginal distributions of X, and hence involves some additional considerations, as discussed in [sec:entropy]Section <ref>.

We now define the class of nonparanormal distributions, from which we assume our data are drawn.

(Nonparanormal distribution, a.k.a. Gaussian copula model)
A random vector X = (X_1,…,X_D)^T is said to have a nonparanormal
distribution (denoted X ∼𝒩𝒫𝒩(Σ; f)) if there exist functions {f_j}_j = 1^D such that each
f_j : → is a diffeomorphism
[A diffeomorphism is a continuously differentiable bijection g : → R ⊆ such that g is continuously differentiable.

]
and f(X) ∼𝒩(0, Σ), for some (strictly) positive definite Σ∈^D × D with 1's on the diagonal (i.e., each σ_j = Σ_j,j = 1).
[Setting [ f(X) ] = 0 and each σ_j = 1 ensures model identifiability, but does not reduce the model space, since these parameters can be absorbed into the marginal transformation f.] Σ is called the latent covariance of X and f is called the marginal transformation of X.



The nonparanormal family relaxes many constraints of the Gaussian family. Nonparanormal distributions can be multi-modal or heavy-tailed, can encode noisy nonlinear dependencies amongst variables, and need not be supported on ^D. Assumptions made by a nonparanormal model on the marginals are minimal; any desired continuously differentiable marginal cumulative distribution function (CDF) F_i of the variable X_i corresponds to the marginal transformation f_i(x) = Φ(F_i(x)) (where Φ is the standard normal CDF). As examples, for a Gaussian variable Z, the 2-dimensional case, X_1∼𝒩(0,1), and X_2 = T(X_1 + Z) is completely captured by a Gaussian copula when T(x) = x^3, T = tanh, T = Φ, or any other diffeomorphism. On the other hand, the limits of the Gaussian copula appear, for example, when T(x) = x^2, which is not bijective; then, if [Z] = 0, the Gaussian copula approximation of (X_1,X_2) will model X_1 and X_2 as independent.

We are now ready to formally state our problem:

Formal Problem Statement:
Given n i.i.d. samples X_1,...,X_n ∼𝒩𝒫𝒩(Σ;f), where Σ and f are both unknown, we would like to estimate I(X).

Other notation: D denotes the dimension of the data (i.e., Σ∈^D × D and f : ^D →^D). For a positive integer k, [k] = {1,...,k} denotes the set of positive integers less than k (inclusive). For consistency, where possible, we use i ∈ [n] to index samples and j ∈ [D] to index dimensions (so that, e.g., X_i,j denotes the j^th dimension of the i^th sample). Given a data matrix X ∈^n × D, our estimators depend on the empirical rank matrix

    R ∈ [n]^n × D   with   
    R_i,j := ∑_k = 1^n 1_{X_i,j≥ X_k,j}.


For a square matrix A ∈^k × k, |A| denotes the determinant of A, A^T denotes the transpose of A, and





    A_2 := max_x ∈^k
    x_2 = 1Ax_2
       and   A_F := √(∑_i, j ∈ [k] A_i,j^2)

denote the spectral and Frobenius norms of A, respectively. When A is symmetric, λ_1(A) ≥λ_2(A) ≥⋯≥λ_D(A) are its eigenvalues.



§ RELATED WORK AND OUR CONTRIBUTIONS




 §.§ The Nonparanormal


Nonparanormal models have been used for modeling dependencies among high-dimensional data in a number of fields, such as graphical modeling of gene expression data <cit.>, of neural data <cit.>, and of financial time series <cit.>, extreme value analysis in hydrology <cit.>, and informative data compression <cit.>.

Besides being more robust generalizations of Gaussians, nonparanormal distributions are also theoretically motivated in certain contexts. For example, the output Z of a neuron is often modeled by feeding a weighted linear combination Y = ∑_k = 1^N w_k X_k of inputs into a nonlinear transformation Z = f(Y). When the components of X are independent, the central limit theorem suggests Y is approximately normally distributed, and hence Z is approximately nonparanormally distributed <cit.>.

With one recent exception <cit.>, previous information estimators for the nonparanormal case <cit.>, rely on fully nonparametric information estimators as subroutines, and hence suffer strongly from the curse of dimensionality. Very recently, <cit.> proposed what we believe is the first mutual information estimator tailored specifically to the nonparanormal case; their estimator is equivalent to one of the estimators (I_G, described in Section <ref>) we study. However, they focused on its applications to neuroimaging data analysis, and did not study its performance theoretically or empirically.



 §.§ Information Estimation

Our motivation for studying the nonparanormal family comes from trying to bridge two recent approaches to information estimation. The first has studied fully non-parametric entropy estimation, assuming only that data are drawn from a smooth probability density p, where smoothness is typically quantified by a Hölder or Sobolev exponent s ∈ (0, ∞), roughly corresponding to the continuous differentiability of s.
In this setting, the minimax optimal MSE rate has been shown by <cit.> to be O ( max{ n, n^-8s/4s + D}).
This rate slows exponentially with the dimension D, and, while many estimators have been proposed <cit.> for this setting, their practical use is limited to a few dimensions[“Few” depends on s and n, but <cit.> suggest nonparametric estimators should only be used with D at most 4-6. <cit.> tried using several nonparametric information estimators on the Communities and Crime UCI data set (n = 2195, D = 10), but found all too unstable to be useful.].

The second area is in the setting where data are assumed to be drawn from a truly Gaussian distribution. Here the high-dimensional case is far more optimistic. While this case had been studied previously <cit.>, <cit.> recently provided a precise finite-sample analysis based on deriving the exact probability law of the log-determinant log|Σ̂| of the scatter matrix Σ̂. From this, they derived a deterministic bias correction, giving an estimator for which they prove an MSE upper bound of 2log( 1 - D/n) and a high-dimensional central limit theorem for the case D →∞ as n →∞ (but D < n).

<cit.> also prove a minimax lower bound of 2D/n on MSE, with several interesting consequences. First, consistent information estimation is possible only if D/n → 0. Second, since, for small x, log(1 - x) ≈ x, this lower bound essentially matches the above upper bound when D/n is small. Third, they show this lower bound holds even when restricted to diagonal covariance matrices. Since the upper bound for the general case and the lower bound for the diagonal case essentially match, it follows that Gaussian information estimation is not made easier by structural assumptions such as Σ being bandable, sparse, or Toeplitz, as is common in, for example, stationary Gaussian process models <cit.>.

This 2D/n lower bound extends to our more general nonparanormal setting. However, we provide a minimax lower bound suggesting that the nonparanormal setting is strictly harder, in that optimal rates depend on Σ. Our results imply nonparanormal information estimation does become easier if Σ is assumed to be bandable or Toeplitz.

A closely related point is that known convergence rates for the fully nonparametric case require the density p to be bounded away from 0 or have particular tail behavior, due to singularity of the logarithm near 0 and resulting sensitivity of Shannon information-theoretic functionals to regions of low but non-zero probability. In contrast, <cit.> need no lower-bound-type assumptions in the Gaussian case. In the nonparanormal case, we show some such condition is needed to prove a uniform rate, but a weaker condition, a positive lower bound on λ_D(Σ), suffices.



The main contributions of this paper are the following:


  * We propose three estimators, Î_G, Î_ρ, and Î_τ,[<cit.> proposed Î_G for use in neuroimaging data analysis. To the best of our knowledge, Î_ρ and Î_τ are novel.] for the mutual information of a nonparanormal distribution.

  * We prove upper bounds, of order O(D^2/(λ_D^2(Σ)n)) on the mean squared error of Î_ρ, providing the first upper bounds for a nonparanormal information estimator. This bound suggests nonparanormal estimators scale far better with D than nonparametric estimators.

  * We prove a minimax lower bound suggesting that, unlike the Gaussian case, difficulty of nonparanormal information estimation depends on the true Σ.

  * We give simulations comparing our proposed estimators to Gaussian and nonparametric estimators. Besides confirming and augmenting our theoretical predictions, these help characterize the settings in which each nonparanormal estimator works best.

  * We present entropy estimators based on Î_G, Î_ρ, and Î_τ. Though nonparanormal entropy estimation requires somewhat different assumptions from mutual information estimation, we show that entropy can also be estimated at the rate O(D^2/(λ_D^2(Σ)n)).















§ NONPARANORMAL INFORMATION ESTIMATORS



In this section, we present three different estimators, I_G, I_ρ, and I_τ, for the mutual information of a nonparanormal distribution. We begin with a lemma providing common motivation for all three estimators.

Since mutual information is invariant to diffeomorphisms of individual variables, it is easy to see that the mutual information of a nonparanormal random variable is the same as that of the latent Gaussian random variable. Specifically:

(Nonparanormal mutual information): Suppose X ∼𝒩𝒫𝒩(Σ; f). Then,

    I(X) = -1/2log|Σ|.




Lemma <ref> shows that mutual information of a nonparanormal random variable depends only the latent covariance Σ; the marginal transformations are nuisance parameters, allowing us to avoid difficult nonparametric estimation; the estimators we propose all plug different estimates of Σ into Eq. (<ref>), after a regularization step described in Section <ref>.



 §.§ Estimating Σ by Gaussianization













The first estimator Σ̂_G of Σ proceeds in two steps. First, the data are transformed to have approximately standard normal marginal distributions, a process <cit.> referred to as “Gaussianization”. By the nonparanormal assumption, the Gaussianized data are approximately jointly Gaussian. Then, the latent covariance matrix is estimated by the empirical covariance of the Gaussianized data.

More specifically, letting Φ denote the quantile function of the standard normal distribution and recalling the rank matrix R defined in (<ref>), the Gaussianized data

    X̃_i,j
      := Φ( R_i,j/n + 1)
         (for  i ∈ [n], j ∈ [D])

are obtained by transforming the empirical CDF of the each dimension to approximate Φ. Then, we estimate Σ by the empirical covariance
Σ̂_G
  := 1/n∑_i = 1^n X̃_i X̃_i^T.

















































 §.§ Estimating Σ by rank correlation


The second estimator actually has two variants, I_ρ and I_τ, respectively based on relating the latent covariance to two classic rank-based dependence measures, Spearman's ρ and Kendall's τ. For two random variables X and Y with CDFs F_X,F_Y : → [0, 1], ρ and τ are defined by

    ρ(X, Y)     := (F_X(X),F_Y(Y)) 
    and   τ(X, Y)     := ((X - X'), (Y - Y')),

respectively, where 
    (X, Y) = [(X - [X])(Y - [Y])]/√([X][Y])
 denotes the standard Pearson correlation operator and (X',Y') is an IID copy of (X,Y). ρ and τ generalize to the D-dimensional setting in the form of rank correlation matrices ρ, τ∈ [-1,1]^D × D with ρ_i,j = ρ(X_i, X_j) and τ_i,j = τ(X_i, X_j) for each i ∈ [n],j ∈ [D].

I_ρ and I_τ are based on a classical result relating the correlation and rank-correlation of a bivariate Gaussian:

<cit.>: Suppose (X,Y) has a Gaussian joint distribution with covariance Σ. Then,

    (X, Y)
     = 2sin(π/6ρ(X, Y) )
     = sin( π/2τ(X, Y) ).

 


ρ and τ are often preferred to Pearson correlation for their relative robustness to outliers and applicability to non-numerical ordinal data. While these are strengths here as well, the main reason for their relevance is that they are invariant to marginal transformations (i.e., for diffeomorphisms f, g : →, ρ(f(X),g(Y)) = ±ρ(X, Y) and τ(f(X),g(Y)) = ±τ(X,Y)). As a consequence, the identity provided in Theorem <ref> extends unchanged to the case (X,Y) ∼𝒩𝒫𝒩(Σ;f). This suggests an estimate for Σ based on estimating ρ or τ and plugging this element-wise into the transform x ↦ 2sin( π/6 x ) or x ↦sin( π/2 x ), respectively. Specifically, Σ_ρ is defined by

    Σ̂_ρ := 2 sin( π/6ρ̂),
         where   ρ̂= (R)

is the empirical correlation of the rank matrix R, and sine is applied element-wise. Similarly, Σ̂_τ := sin( π/2τ̂),
where

    τ̂_j,k := 1/n2∑_i ≠ℓ∈ [n](X_i,j - X_ℓ,j)(X_i,k - X_ℓ,k).




 §.§ Regularization and estimating I


Unfortunately, unlike usual empirical correlation matrices, none of Σ̂_G, Σ̂_ρ, or Σ̂_τ is almost surely strictly positive definite. As a result, directly plugging into the mutual information functional (<ref>) may give ∞ or even be undefined.

To correct for this, we propose a regularization step, in which we project each estimated latent covariance matrix onto the (closed) cone 𝒮(z) of symmetric matrices with minimum eigenvalue z > 0. Specifically, for any z > 0, let

    𝒮(z) := { A ∈^D × D : A = A^T, λ_D(A) ≥ z }.

For any symmetric matrix A ∈^D × D with eigendecomposition Σ̂= Q Λ Q (i.e., QQ^T = Q^TQ = I_D and Λ is diagonal), the projection A_z of A onto 𝒮(z) is defined as A_z := Q Λ_z Q, where Λ_z is the diagonal matrix with j^th nonzero entry ( Λ_z )_j,j = max{ z, Λ_j,j}. We call this a “projection” because A_z is precisely the Frobenius norm projection of A onto 𝒮(z) (see, e.g., <cit.>):
A_z = min_B ∈^D × DA - B_F.

Applying this regularization to Σ̂_G, Σ̂_ρ, or Σ̂_τ gives a strictly positive definite estimate Σ̂_G,z, Σ̂_ρ,z, or Σ̂_τ,z, respectively, of Σ. We can then estimate I by plugging this into Equation (<ref>), giving our three estimators:

    Î_G,z
    := -1/2log| Σ̂_G,z|,
        Î_ρ,z
    := -1/2log| Σ̂_ρ,z| 
     and   Î_τ,z
    := -1/2log| Σ̂_τ,z|.




§ UPPER BOUNDS ON THE ERROR OF Î_Ρ,Z







Here, we provide finite-sample upper bounds on the error of the estimator Î_ρ based on Spearman's ρ. Proofs are given in the Appendix.

We first bound the bias of the estimator:

Suppose X_1,...,X_n i.i.d.∼𝒩𝒫𝒩(Σ;f). Then, there exists a constant C > 0 such that, for any z > 0, the bias of Î_ρ,z is at most

    | [ Î_ρ,z] - I |
    ≤ C ( D/z√(n) + log|Σ_z|/|Σ|),

where Σ_z is the projection of Σ onto 𝒮(z).



The first term of the bias stems from nonlinearity of the log-determinant function in Equation <ref>, which we analyze via Taylor expansion. The second term,

    log|Σ_z|/|Σ|
      = ∑_λ_j(Σ) < zlog( z/λ_j(Σ)),

is due to the regularization step and is actually, but is difficult to simplify or bound without further assumptions on the spectrum of Σ and a choice of z, which we discuss later. We now turn to bounding the variance of Î_ρ,z. We first provide an exponential concentration inequality for Î_ρ,z around its expectation, based on McDiarmid's inequality:

Suppose X_1,...,X_n i.i.d.∼𝒩𝒫𝒩(Σ;f). Then, for any z, > 0,

    [
        |
          Î_ρ,z - [ Î_ρ,z]
        | > ]
      = 2 exp(
        - n z^2^2/18 π^2 D^2).




Such exponential concentration bounds are useful when one wants to simultaneously bound the error of multiple uses of an estimator, and hence we present it separately as it may be independently useful. However, for the purpose of understanding convergence rates, we are more interested in the variance bound that follows as an easy corollary:


Suppose X_1,...,X_n i.i.d.∼𝒩𝒫𝒩(Σ;f). Then, for any z > 0, the variance of Î_ρ,z is at most

    [ Î_ρ,z]
      ≤36 π^2 D^2/z^2 n.




Given these bias and variance bounds, a bound on the MSE of Î_ρ,z follows via the usual bias-variance decomposition:

Suppose X ∼𝒩𝒫𝒩(Σ;f). Then, there exists a constant C such that






    [ ( Î_ρ,z - I )^2 ]
    ≤ C ( D^2/z^2n + log^2 |Σ_z|/|Σ|).



A natural question is now how to optimally select the regularization parameter z. While the bound (<ref>) is clearly convex in z, it depends crucially on the unknown spectrum of Σ, and, in particular, on the smallest eigenvalues of Σ. As a result, it is difficult to choose z optimally in general, but we we can do so for certain common subclasses of covariance matrices. For example, if Σ is Toeplitz or bandable (i.e., for some c ∈ (0,1), all |Σ_i,j| ≤ c^|i - j|), then the smallest eigenvalue of Σ can be bounded below <cit.>. When Σ is bandable, as we show in the Appendix, this bound can be independent of D. In these cases, the following somewhat simpler MSE bound can be used:

Suppose X ∼𝒩𝒫𝒩(Σ;f), and suppose z ≤λ_D(Σ). Then, there exists a constant C > 0 such that

    [ ( Î_ρ,z - I )^2 ] ≤CD^2/z^2n.



















§ LOWER BOUNDS IN TERMS OF Σ


When the data X_1,...,X_n i.i.d∼𝒩(0,Σ) are truly Gaussian, using the plug-in estimator

    Î = -1/2log| Σ̂|
        Σ̂= 1/n∑_i = 1^n X_i X_i^T

is the empirical covariance matrix), <cit.> showed that the distribution of Î - I is independent of the true correlation matrix Σ. This follows from  the “stability” of Gaussians (i.e., that nonsingular linear transformations of Gaussian random variables are Gaussian). In particular,

    Î - I = log|Σ̂| - log|Σ| = log|Σ^-1/2Σ̂Σ^-1/2|,
 and Σ^-1/2Σ̂Σ^-1/2 has the same distribution as logΣ̂ does in the special case that Σ = I_D is the identity. This property is both somewhat surprising, given that I →∞ as |Σ| → 0, and useful, leading to a tight analysis of the error of Î and confidence intervals that do not depend on Σ. 

It would be convenient if any nonparanormal information estimators satisfied this property. Unfortunately, the main result of this section is a negative one, showing that this property is unlikely to hold without additional assumptions:


Consider the 2-dimensional case

    X_1,...,X_n i.i.d∼𝒩(0,Σ),    with   Σ = [ 1 σ; σ 1 ],

and let σ_* ∈ (0,1). Suppose an estimator Î = Î(R) of I_σ = -1/2log(1 - σ^2) is a function of the empirical rank matrix R ∈^n × 2 of X. Then, there exists a constant C > 0, depending only n, such that the worst-case MSE of Î over σ∈ (0,σ_*) satisfies

    sup_σ∈ (0,σ^*)[ ( Î(R) - I_σ)^2 ]
       ≥1/64( C - log(1 - σ_*^2) )^2



Clearly, this lower bound tends to ∞ as σ→ 1.
As written, this result lower bounds the error of rank-based estimators in the Gaussian case when σ≈ 1.

However, to the best of our knowledge, all methods for estimating Σ in the nonparanormal case are functions of R, and prior work <cit.> has shown that the rank matrix R is a generalized sufficient statistic for Σ (and hence for I) in the nonparanormal model. Thus, it is reasonable to think of lower bounds for rank-based estimators in the Gaussian case as lower bounds for any estimator in the nonparanormal case.

The proof of this result is based on the simple observation that the rank matrix can take only finitely many values. Hence, as σ→ 1, R tends to be perfectly correlated, providing little information about σ, whereas the dependence of the estimand I_σ on σ increases sharply. This is intuition is formalized in the Appendix
using Le Cam's lemma for lower bounds in two-point parameter estimation problems.



§ EMPIRICAL RESULTS


We compare 5 mutual information estimators:


  * Î: Gaussian plug-in estimator with bias-correction (see <cit.>).

  * Î_G: Nonparanormal estimator using Gaussianization.

  * Î_ρ: Nonparanormal estimator using Spearman's ρ.

  * Î_τ: Nonparanormal estimator using Kendall's τ.

  * Î_kNN: Nonparametric estimator using k-nearest neighbor (kNN) statistics.


For I_ρ and I_τ, we used a regularization constant z = 10^-3. We did not regularize for I_G. Although this implies [I_G = ∞] > 0, this is extremely unlikely for even moderate values of n and never occurred during our experiment, which all use n ≥ 32. We will thus omit denoting dependence on z. For I_kNN, except as noted in Experiment 3, k = 2, based on recent analysis <cit.> suggesting that small values of k are best for estimation.

Sufficient details to reproduce experiments are given in the Appendix,

and MATLAB source code is available at [Omitted for anonymity]. We report MSE based on 1000 i.i.d. trials of each condition. 95% confidence intervals were consistently smaller than plot markers and hence omitted to avoid cluttering plots. Except as specified otherwise, each experiment had the following basic structure: In each trial, a correlation matrix Σ was drawn by normalizing a random covariance matrix from a Wishart distribution, and data X_1,...,X_n i.i.d.∼𝒩(0, Σ) drawn. All 5 estimators were computed from X_1,...,X_n and squared error from true mutual information (computed from Σ) was recorded. Unless specified otherwise, n = 100 and D = 25.

Since our nonparanormal information estimators are functions of ranks of the data, neither the true mutual information nor our non-paranormal estimators depend on the marginal transformations. Thus, except in Experiment 2, where we show the effects of transforming marginals, and Experiment 3, where we add outliers to the data, we perform all experiments on truly Gaussian data, with the understanding that this setting favors the Gaussian estimator.

All experimental results are displayed in Figure <ref>.











Experiment 1 (Dependence on n):
We first show nonparanormal estimators have “parametric” O(n) dependence on n, unlike Î_kNN, which converges far more slowly. For large n, MSEs of Î_G, Î_ρ, and Î_τ are close to that of Î.


Experiment 2 (Non-Gaussian Marginals):






Next, we show nonparanormal estimators are robust to non-Gaussianity of the marginals, unlike Î. We applied a nonlinear transformation f to a fraction α∈ [0, 1] of dimensions of Gaussian data. That is, we drew Z_1,...,Z_n i.i.d.∼𝒩(0,Σ) and then used data X_1,...,X_n, where

    X_i,j = {[ T(Z_i,j)  j < α D;    Z_i,j   j ≥α D ].,   ∀ i ∈ [n], j ∈ [D],

for a diffeomorphism T. Here, we use T(z) = e^z. The Appendix shows similar results for several other T. Î performs poorly even when α is quite small. Poor performance of Î_kNN may be due to discontinuity of the density at x = 0.

Experiment 3 (Outliers):





We now show that nonparanormal estimators are far more robust to the presence of outliers than Î or Î_kNN. To do this, we added outliers to the data according to the method of <cit.>. After drawing Gaussian data, we independently select ⌊β n ⌋ samples in each dimension, and replace each i.i.d. uniformly at random from {-5,+5}. Performance of Î degrades rapidly even for small β. Î_kNN can fail for atomic distributions, Î_kNN = ∞ whenever at least k samples are identical. This mitigate this, we increased k to 20 and ignored trials where Î_kNN = ∞, but Î_kNN ceased to give any finite estimates when β was sufficiently large.

For small values of β, nonparanormal estimators surprisingly improve. We hypothesize this is due to convexity of the mutual information functional Eq. (<ref>) in Σ. By Jensen's inequality, estimators which plug-in an approximately unbiased estimate Σ̂ of Σ are biased towards overestimating I. Adding random (uncorrelated) noise reduces estimated dependence, moving the estimate closer to the true value.

If this nonlinearity is indeed a major source of bias, it may be possible to derive a von Mises-type bias correction (see <cit.>) accounting for higher-order terms in the Taylor expansion of the log-determinant.








Experiment 4 (Dependence on Σ):
Here, we verify our results in [sec:Sigma_lower_bound]Section <ref> showing that MSE of rank-based estimators approaches ∞ as |Σ| → 0, while MSE of Î is independent of Σ. Here, we set D = 2 and Σ as in Eq. (<ref>), varying σ∈ [0,1]. Indeed, the MSE of Î does not change, while the MSEs of Î_G, Î_ρ, and Î_τ all increase as σ→ 1. This increase seems mild in practice, with performance worse than of Î only when σ > 0.99. Î_τ appears to perform far better than Î_G and Î_ρ in this regime.

Performance of I_kNN degrades far more quickly as σ→ 1. This phenomenon is explored by <cit.>, who lower bound error of I_kNN in the presence of strong dependencies, and proposed a correction to improve performance in this case.

It is also interesting that errors of Î_ρ and Î_τ drop as σ→ 0. This is likely because, in this regime, the main source of error is the variance of ρ̂ and τ̂ (as -log(1 - σ^2) ≈σ^2 when σ≈ 0). When n →∞ and D is fixed, both 2sin(πρ̂/6) and sin(πτ̂/2) are asymptotically normal estimates of σ, with asymptotic variances proportional to (1 - σ^2)^2 <cit.>. By the delta method, since dI/dσ = σ/1 - σ^2, Î_ρ and Î_τ are asymptotically normal estimates of I, with asymptotic variances proportional to σ^2 and hence vanishing as σ→ 0.














§ ESTIMATING ENTROPY


Thus far, we have discussed estimation of mutual information I(X). Mutual information is convenient because it is invariant under marginal transformation, and hence I(X) = I(f(X)) depends only on Σ. While the entropy H(X) does depend on the marginal transform f, fortunately, by Eq. (<ref>), H(X) differs from I(X) only by a sum of univariate entropies. Univariate nonparametric estimation of entropy in has been studied extensively, and there exist several estimators (e.g., based on sample spacings <cit.>, kernel density estimates <cit.> or k-nearest neighbor methods <cit.>) that can estimate H(X_j) at the rate ≍ n in MSE under relatively mild conditions on the marginal density p_j. While the precise assumptions vary with the choice of estimator, they are mainly (a) that p_j be lower bounded on its support or have particular (e.g., exponential) tail behavior, and (b) that p_j be smooth, typically quantified by a Hölder or Sobolev condition. Details of these assumptions are in the Appendix. 

Under these conditions, since there exist estimators Ĥ_1,...,Ĥ_D and a constant C > 0 such that

    [(Ĥ_j - H(X_j))^2] ≤ C/n,
      ∀ j ∈ [D].

Combining these estimators with an estimator, say Î_ρ,z, of mutual information gives an estimator of entropy:

    Ĥ_ρ,z := ∑_j = 1^D Ĥ_j - Î_ρ,z.

If we assume z = λ_D(Σ) is bounded below by a positive constant, combining inequality (<ref>) with Corollary <ref> gives

    [ ( Ĥ_ρ,z - H(X) )^2 ] ≤CD^2/n,

where the constant C may differ from in (<ref>) but is independent of n and D.










§ CONCLUSIONS AND FUTURE WORK



This paper we suggests nonparanormal information estimation as a practical compromise between the difficult nonparametric case and the restrictive Gaussian case. We proposed three estimators for this problem, and provided the first upper bounds for nonparanormal information estimation. We also provided lower bounds showing how dependence on Σ differs from the Gaussian case, and we demonstrated empirically that nonparanormal estimators are more robust than Gaussian estimators, even when dimension is too high for fully nonparametric estimators.

Collectively, these results suggest that, by scaling to moderate or high dimensionality without relying on Gaussianity, nonparanormal information estimators may be effective tools with a number of machine learning applications. While the best choice of information estimator inevitably depends on context, as a crude off-the-shelf guide for practitioners, the estimators we might suggest, in order of preference, are:


  * fully nonparametric if D < 6, n > max{100,10^D}.

  * Î_ρ if D^2/n is small and data may have outliers.

  * Î_τ if D^2/n is small and dependencies may be strong.

  * Î_G otherwise.

  * Î only given strong belief that data are nearly Gaussian.












There are many natural open questions in this line of work. First, in the nonparanormal model, we focused on estimating mutual information I(X), which does not depend on marginal transforms f, and entropy, which decomposes into I(X) and 1-dimensional entropies. In both cases, additional structure imposed by the nonparanormal model allows estimation in higher dimensions than fully nonparametric models. Can nonparanormal assumptions lead to higher dimensional estimators for the many other useful nonlinear functionals of densities (e.g., L_p norms/distances and more general (e.g., Rényi or Tsallis) entropies, mutual informations, and divergences) that do not decompose?

Second, there is a gap between our upper bound rate of Σ_2^2 D^2/n and the only known lower bound of 2D/n (from the Gaussian case), thought we also showed that bounds for rank-based estimators depend on Σ. Is quadratic dependence on D optimal? How much do rates improve under structural assumptions on Σ? Upper bounds should be derived for other estimators, such as Î_G and Î_τ. The 2D/n lower bound proof of <cit.> for the Gaussian case, based on the Cramer-Rao inequality <cit.>, is unlikely to tighten in the nonparanormal case, since Fisher information is invariant to diffeomorphisms of the data. Hence, a new approach is needed if the lower bound in the nonparanormal case is to be raised.





Finally, our work also applies to estimating the log-determinant log|Σ| of the latent correlation matrix in a nonparanormal model. In addition to information estimation, the work of <cit.> on estimating log|Σ| in the Gaussian setting was motivated by the use of log|Σ| in several other multivariate statistical tools, including

quadratic discriminant analysis (QDA) and

MANOVA <cit.>. Can our estimators lead to more robust nonparanormal versions of these procedures?


icml2017





§ LEMMAS

Our proofs rely on the following lemmas.

(Convexity of the inverse operator norm):
The function A ↦A_2 is convex over A ≻ 0.



For A, B ≻ 0, let C := τ A + (1 - τ) B. Then,

    Ĉ_2
        = 1/inf_x ∈^D x^T C x
        = 1/inf_x ∈^Dτ x^T A x + (1 - τ) x^T B x
       ≤1/τinf_x ∈^D x^T A x + (1 - τ) inf_x ∈^D x^T B x
       ≤τ1/inf_x ∈^D x^T A x
      + (1 - τ) 1/inf_x ∈^D x^T B x
        = τ A_2 + (1 - τ)  B_2

via convexity of the function x ↦ 1/x on (0, ∞).



(Mean-Value Bound on the Log-Determinant):
Matrix derivative of log-determinant. Suppose A, B ≻ 0. Then, for λ := min{λ_D(A), λ_D(B)},

    | log |A| - log |B| | ≤1/λA - B_F.




Proof: First recall that the log-determinant is continuously differentiable over the strict positive definite cone, with ∇_X log |X| = X for any X ≻ 0. Hence, by the matrix-valued version of the mean value theorem,

    log |A| - log |B| = tr(C(A - B)),

where C = τ A + (1 - τ) B for some τ∈ (0, 1). Since for positive definite matrices, the inner product can be bounded by the product of the operator and Frobenius norms, and clearly C ≻ 0, we have

    | log |A| - log |B| | = C_2A - B_F.

Finally, it follows by Lemma <ref> that

    | log |A| - log |B| | ≤1/λA - B_F.





§ PROOFS OF MAIN RESULTS





























































































Here, we give proofs of our main theoretical results, beginning with upper bounds on the MSE of Î_ρ and proceeding to minimax lower bounds in terms of Σ.



§ UPPER BOUNDS ON THE MSE OF Î_Ρ




    | [ log |Σ̂_z| ] - log|Σ| | 
       ≤ C ( Σ_2^2D/z^2n + ( ∑_λ_j(Σ) < zlog( z/λ_j(Σ)) )^2 ).




By the triangle inequality,

    | [ log |Σ̂_z| ] - log|Σ| |
       ≤| [ log |Σ̂_z| ] - log|Σ_z| | 
          + | log |Σ_z| - log |Σ| |

For the first term, applying the matrix mean value theorem (Lemma <ref>) and the inequality
A_F ≤√(D)A_2

    | [ log| Σ̂_z | ] - log |Σ_z| |
       ≤[ | log| Σ̂_z | - log |Σ_z| | ] 
       ≤1/z[ Σ̂_z - Σ_z _F ] 
       ≤√(D)/z[ Σ̂_z - Σ_z _2 ] 
       ≤C_MZΣ_2 D/z√(n),

where we used Theorem 1 of <cit.>, which gives a
constant C_MZ such that

    [ Σ̂_z - Σ_z _2 ] ≤ C_MZΣ_2 √(D/n).

Via the bound Σ_2 ≤√(D)Σ_∞, this reduces to

    [ Σ̂_z - Σ_z _2 ]
      ≤ C_MZD/√(n).














    [ Î]
      ≤36 π^2 D^2/z^2 n.





By the Efron-Stein inequality, since X_1,…,X_n are independent and identically distributed,

    [ Î]
       ≤1/2∑_i = 1^n
              [
                (
                  log |Σ̂_z| - log |Σ̂_z^(i)|
                )^2
              ] 
        = n/2[
                (
                  log |Σ̂_z| - log |Σ̂_z^(1)|
                )^2
              ],

where Σ̂_z^(1) is our estimator after independently re-sampling the
first sample X_1. Applying the multivariate mean-value theorem (Lemma <ref>), we have

    | log |Σ̂_z| - log |Σ̂_z^(1)| |
      ≤1/zΣ̂_z - Σ̂_z^(1)_F.

Σ̂_τ_2 ≤1/z.
Since 𝒮(z) is convex and the Frobenius norm is supported by an inner product, the operation of projecting onto 𝒮(z) is a contraction. In particular,
( Σ̂_z - Σ̂_z^(1)) _F
  ≤( Σ̂- Σ̂^(1)) _F
Applying the mean value theorem to the function x ↦ 2sin( π/6 x ),

    ( Σ̂- Σ̂^(1)) _F^2
        = ∑_j,k = 1^D ( Σ̂- Σ̂^(1))_j,k^2 
       ≤π^2/9∑_j,k = 1^D ( ρ̂_j,k - ρ̂_j,k^(1))^2 
        = π^2/9ρ̂- ρ̂^(1)_F^2.

From the formula

    ρ̂_j,k = 1 - 6 ∑_i = 1^n d_i,j,k^2/n(n^2 - 1),

(where d_i,j,k denotes the difference in ranks of X_i,j and X_i,k in X_1,j,...,X_n,j and X_1,k,...,Y_n,k, respectively), one can see, since |d_1,j,k - d_1,j,k'| ≤ n and, for i ≠ 1, |d_i,j,k - d_i,j,k'| ≤ 1, that

    | ρ̂_j,k - ρ̂_j,k^(1)| ≤18/n,

and hence that

    ρ̂- ρ̂^(1)_F ≤18D/n.

It follows from inequality (<ref>) that

    Σ̂_z - Σ̂_z^(1)_F ≤6 π D/n.

Altogether, this gives

    | log |Σ̂_z| - log |Σ̂_z^(1)| |
      ≤6 π D/z n.

Then, McDiarmid's Inequality gives, for all > 0,

    [ | Î - [ Î] | > ]
      = 2 exp(
        - n z^2^2/18 π^2 D^2).

This translates to a variance bound of

    [ Î]
      ≤36 π^2 D^2/z^2 n.





 §.§ Lower bound for rank-based estimators in terms of Σ

One (perhaps surprising) result of <cit.> is that, as long as D/n → 0, the convergence rate of the estimator is independent of the true correlation structure Σ. Here, we show that this desirable property does not hold in the nonparanormal case.


Consider the 2-dimensional case

    X_1,...,X_n i.i.d∼𝒩(0,Σ),    with   Σ = [ 1 σ; σ 1 ],

and let σ_* ∈ (0,1). Suppose an estimator Î = Î(R) of I_σ = -1/2log(1 - σ^2) is a function of the empirical rank matrix R ∈^n × 2 of X (as defined in (<ref>)). Then, there exists a constant C > 0, depending only n, such that the worst-case MSE of Î over σ∈ (0,σ_*) satisfies

    sup_σ∈ (0,σ^*][ ( Î(R) - I_σ)^2 ]
       ≥1/64( C - log(1 - σ_*^2) )^2 
       →∞   as   σ_* → 1.




Note that the rank matrix R can take only finitely many values. Let ℛ be the set of all (n!)^D possible rank matrices and let ℛ_1 ⊆ℛ be the set of rank matrices that are perfectly correlated. Then, as σ→ 1, [R ∈ℛ_1] → 1, so, in particular, we can pick σ_0 (depending only on n) such that, for all σ≥σ_0, [R ∈ℛ_1] ≥1/2. Since the data are i.i.d., all rank matrices in ℛ_1 have equal probability. It follows that

    D_TV(_0||_1)
     = 1/2_0 - _1_1
     ≤1/2.

Finally, by Le Cam's Lemma (see, e.g., Section 2.3 of <cit.>),

    inf_Îsup_σ∈{σ_0,σ_1}[ ( Î - I_σ)^2 ] 
       ≥(I_σ_* - I_σ_0)^2/8( 1 - D_TV(P_σ_0,P_σ_1) ) 
       ≥(log(1 - σ_0^2) - log(1 - σ_*^2))^2/64





§ DETAILS OF EXPERIMENTAL METHODS


Here, we present details needed to reproduce our numerical simulations. Note that MATLAB source code for these experiments is available at [Omitted for anonymity.], including a single runnable script that performs all experiments and generates all figures presented in this paper.
Specific details needed to reproduce experiments are given in the Appendix,


In short, experiments report empirical mean squared errors based on 100 i.i.d. trials of each condition. We initially computed 95% confidence intervals, but these intervals were consistently smaller than marker sizes, so we omitted them to avoid cluttering plots. Except as specified otherwise, each experiment followed the same basic structure, as follows: In each trial, a random correlation matrix Σ∈ [-1,1]^D × D was drawn by normalizing a covariance matrix from a Wishart distribution W(I_D,D) with identity scale matrix and D degrees of freedom. Data X_1,...,X_n were then drawn i.i.d. from 𝒩(0, Σ). All estimators were applied to the same data. Unless specified otherwise, n = 100 and D = 25.





 §.§ Computational Considerations

In general, the running time of all the nonparanormal estimators considered is O(D n log n + D^2 n + D^3) (i.e., O(Dn log n) to rank or Gaussianize the variables in each dimension, D^2n to compute the covariance matrix, and O(D^3) to compute the log-determinant). All log-determinants log|Σ| were computed by summing the logarithms of the diagonal of the Cholesky decomposition of Σ, as this is widely considered to be a fast and numerically stable approach. Note however that faster (O(D)-time) randomized algorithms <cit.> have been proposed to approximate the log-determinant).



§ ADDITIONAL EXPERIMENTAL RESULTS


Here, we present variants on the experiments presented in the main paper, which support but are not necessary for illustrating our conclusions.


 §.§ Effects of Other Marginal Transformations


In [subsec:hat_I_nonrobust]Section <ref>, we showed that the Gaussian estimator Î is highly sensitive to failure of the Gaussian assumption for even a small fraction of marginals. Figure <ref>, illustrates this for the transformation x ↦exp(x), but we show here that this is not specific to the exponential transformation. As shown in Figures <ref> nearly identical results hold when the marginal transformation f is the hyperbolic tangent function x ↦tanh(x), the cubic function x ↦ x^3, sigmoid function x ↦1/1 + e^-x, or standard normal CDF.












































































































































































































































































§ SPECIFIC ASSUMPTIONS FOR ESTIMATING H(X)


As shown in the main paper, to estimate the entropy of a nonparanormal distribution at the rate O(D^2/n), it suffices to the univariate entropy of each variable X_j at the rate O(1/n). To do this, additional assumptions are required on the marginal densities p_j. Here, we give detailed sufficient conditions for this.

Letting S_j ⊆ denote the support of p_j, the two key assumptions can be roughly classified as follows:



  * 1/2-order smoothness[This is stronger than the 1/4-order smoothness mandated by the minimax rate for entropy estimation <cit.>, but appears necessary for most practical entropy estimators. See Section 4 of <cit.> for further details.]; e.g., a Hölder condition:

    0pt
    sup_x ≠ y ∈ S_j|p_j(x) - p_j(y)|/|x - y|^1/2
     < L,

or a (slightly weaker) Sobolev condition:

    0pt
    ∫_S_j p_j^2(x)   dx < ∞  and ∫_S_j( |ξ|^1/2 |ℱ[ p_j ](ξ)| )^2 dξ < L,

(where ℱ[ p_j ](ξ) denotes the Fourier transform of p_j evaluated at ξ) for some constant L > 0.

  * absolute bounds p_j(x) ∈ [κ_1,κ_2] for all x ∈ S_j or (a_j,b_j)-exponential tail bounds

    f(x)/exp(-a_j x^b_j)∈ [κ_1,κ_2]    for all  x ∈ S_j

for some κ_1, κ_2 ∈ (0, ∞).


Under these assumptions, there are a variety of nonparametric univariate entropy estimators that have been shown to converge at the rate O(1/n) <cit.>.



§ LOWER BOUNDING THE EIGENVALUES OF A BANDABLE MATRIX

Recall that, for c ∈ (0,1), a matrix Σ∈^D × D is called c-bandable if there exists a constant c ∈ (0,1) such that, for all i,j ∈ D, |Σ_i,j| ≤ c^|i - j|.

Here, we show simple bounds on the eigenvalues of a bandable correlation matrix Σ. While this result is fairly straightforward, a brief search the literature turned up no comparable results. <cit.>, who originally introduced the class of bandable covariance matrices, separately assumed the existence of lower and upper bounds on the eigenvalues to prove their results. In the context of information estimation, this results of particular interest because, when c < 1/3 it implies a dimension-free positive lower bound on the minimum eigenvalue of Σ, hence complementing our upper bound in Theorem <ref>.


Suppose a symmetric matrix Σ∈^D × D is c-bandable and has identical diagonal entries Σ_j,j = 1. Then, the eigenvalues λ_1(Σ),...,λ_D(Σ) of Σ can be bounded as

    1 - 3c/1 - c≤λ_1(Σ),...,λ_D(Σ) ≤1 + c/1 - c.

In particular, when c < 1/3, we have

    0 < 1 - 3c/1 - c≤λ_D(Σ).




The proof is based on the Gershgorin circle theorem <cit.>. In the case of a real symmetric matrix Σ, this states that the eigenvalues of Σ lie within a union of intervals

    {λ_1(Σ),...,λ_D(Σ) }⊆⋃_j = 1^D [ Σ_j,j - R_j, Σ_j,j + R_j ],

where R_j := ∑_k ≠ j |Σ_j,k| is the sum of the absolute values of the non-diagonal entries of the j^th row of Σ. In our case, since the diagonal entries of Σ are all Σ_j,j = 1, we simply have to bound

    max_j ∈ [D] R_j ≤∑_k ≠ j c^|k - j|.

This geometric sum is maximized when j = ⌈ D/2 ⌉, giving

    R_j
     ≤ 2∑_δ = 1^⌊ D/2 ⌋ c^δ
     = 2c1 - c^⌊ D/2 ⌋/1 - c≤2c/1 - c.

Finally, the inclusion (<ref>) gives

    λ_D(Σ) ≥ 1 - 2c/1 - c = 1 - 3c/1 - c > 0

when c < 1/3.
1 + 2c/1 - c = 1 + c/1 - c.










