Minimum Mutual Information and Non-Gaussianity Through the Maximum Entropy Method: Theory and Properties

Abstract

The application of the Maximum Entropy (ME) principle leads to a minimum of the Mutual Information (MI), I(X,Y), between random variables X,Y, which is compatible with prescribed joint expectations and given ME marginal distributions. A sequence of sets of joint constraints leads to a hierarchy of lower MI bounds increasingly approaching the true MI. In particular, using standard bivariate Gaussian marginal distributions, it allows for the MI decomposition into two positive terms: the Gaussian MI (I_g), depending upon the Gaussian correlation or the correlation between ‘Gaussianized variables’, and a non‑Gaussian MI (I_ng), coinciding with joint negentropy and depending upon nonlinear correlations. Joint moments of a prescribed total order p are bounded within a compact set defined by Schwarz-like inequalities, where I_ng grows from zero at the ‘Gaussian manifold’ where moments are those of Gaussian distributions, towards infinity at the set’s boundary where a deterministic relationship holds. Sources of joint non-Gaussianity have been systematized by estimating I_ng between the input and output from a nonlinear synthetic channel contaminated by multiplicative and non-Gaussian additive noises for a full range of signal-to-noise ratio (snr) variances. We have studied the effect of varying snr on I_g and I_ng under several signal/noise scenarios.

1. Introduction

One of the most commonly used information theoretic measures is the mutual information (MI) [1], measuring the total amount of probabilistic dependence among random variables (RVs)—see [2] for a unifying perspective and axiomatic review. MI is a positive quantity vanishing iff RVs are independent.

MI is an effective tool for multiple purposes, namely and among others: (a) Blind signal separation [3] and Independent Component Analysis (ICA) [4], both of which look for transformed and/or lagged stochastic time series [5] which minimize MI; (b) Predictability studies, Predictable Component Analysis [6] and Forecast Utility [7], all of which are focused on the analysis and decomposition of the MI between probabilistic forecasts and observed states.

Analytical expressions of MI are known for a few number of parametric joint distributions [8,9]. Alternatively, it can be numerically estimated by different methods, such as Maximum Likelihood estimators, Edgeworth expansion, Bayesian methods, equiprobable and equidistant histograms, kernel‑based probability distribution functions (PDFs), K-nearest neighbors technique—see [10,11] and references therein for a survey of estimation methods and scoring comparison studies.

In the bivariate case, treated here, the MI $I(X,Y)$, between RVs $X,Y$is the Kullback–Leibler (KL) divergence: $I(X,Y)\equiv D\left(p_{XY}\left|\right|p_X{p}_Y\right)=E_{p_{XY}}\left(\log (p_{XY}/(p_X{p}_Y))\right)\ge 0$, between the joint probability function $p_{XY}$ and the product of marginal probability distributions $p_X{p}_Y$ where $E_{p_{XY}}$ is the expectation operator over the measure $p_{XY}$. The MI is invariant for smooth invertible transformations of $X,Y$.

The goal is the determination of theoretical lower MI bounds under certain conditions or, in other words, the minimum mutual information (MinMI) [12] between two RVs $X,Y$, consistent, both with imposed marginal distributions and cross-expectations assessing their linear and nonlinear covariability. Those lower bounds can be obtained due to the application of the Maximum Entropy (ME) method to distributions [13] and to the inequality: $D\left(p_{ME}\left|\right|q\right)\le D\left(p\left|\right|q\right)$. Here, $p_{ME}$ is a Maximum Entropy probability distribution (MEPD) with respect to q, say $p_{ME}=\arg {\min }_{p\in \mathrm{\Omega }}\left(D(p||q)\right)$, where Ω is a PDF class verifying a given set of constraints [1]. Therefore, by using the joint MEPD $p=p_{XY-ME}$ and $q=p_X{p}_Y$, the lower MI bound is obtained: $D\left(p_{XY-ME}\left|\right|p_X{p}_Y\right)\le D\left(p_{XY}\in \mathrm{\Omega }\left|\right|p_X{p}_Y\right)\le I(X,Y)$. Finding of the bivariate $p_{XY-ME}$ is straightforward if the marginal distributions are themselves well defined ME distributions. We solve that by transforming the single variables $X,Y$ into others with imposed ME probability mass distributions through the so called ME-anamorphoses [14].

The joint ME probability distribution $p_{XY-ME}$ is derived from the minimum of a functional in terms of a certain number of Lagrange parameters. The properties of multivariate ME distributions have been studied for various ME constraints, namely: (a) imposed marginals and covariance matrix [15]; (b) generic joint moments [16]. Abramov [17,18,19] has developed efficient and stable numerical iterative algorithms for computing ME distributions forced by sets of polynomial expectations. Here we use a bivariate version of the algorithm of [20], already tested in [21].

By taking a sequence of encapsulated sets of joint ME constraints we obtain an increasing hierarchy of lower MI bounds converging towards the total MI.

We particularize this methodology to the case where $X,Y$ are standard Gaussians, issued from single homeomorphisms of an original pair of variables $\widehat{X},\widehat{Y}$ by Gaussian anamorphosis [14]. Then we get the MI $I(X,Y)=I(\widehat{X},\widehat{Y})$, which is decomposed into two generic positive quantities, a Gaussian MI I_g and a non-Gaussian MI I_ng [21], vanishing under bivariate Gaussianty. The Gaussian MI is given by $I_g(c_g)\equiv -1/2\log (1-c_g^2)$, where $c_g(\widehat{X},\widehat{Y})\equiv cor(X,Y)$ is the Pearson linear correlation between the Gaussianized variables. As for the non-Gaussian MI term I_ng it relies upon imposed nonlinear correlations between the single Gaussian variables. The MI reduces to $I_g(c_g)$ when only moments of order one and two are imposed as ME constraints. We will note that, for certain extreme non-Gaussian marginals, $I_g(c=cor(\widehat{X},\widehat{Y}))>I(\widehat{X},\widehat{Y})$, thus showing that $I_g(c)$ [22] is not a proper MI lower bound in general.

The correlation $c_g$, hereafter called Gaussian correlation [21] is a nonlinear concordance measure like the Spearman rank correlation [23] and the Kendall τ, which, by definition are invariant for monotonically growing smooth marginal transformations. These measures have the good property of being expressed as functionals of the copula density functions [24], uniquely dependent on the cross dependency between variables.

The non-Gaussian MI _ng holds some interesting characteristics. It coincides with the joint negentropy (deficit of entropy with respect to that of the Gaussian PDF with the same mean, variance and covariance) in the space of ‘Gaussianized’ variables, which is invariant for any orthogonal or oblique rotation of them. In particular for uncorrelated rotated variables, it coincides with the ‘compactness’, which measures the concentration of the joint distribution around a lower-dimensional manifold [25] which is given by $D\left(p_{XY}\left|\right|p_{SG}\right)$, i.e., the KL divergence with respect to the spherical Gaussian $p_{SG}$ (Gaussian with an isotropic covariance matrix with the same trace as that of $p_{XY}$, say the total variance).

We also show that I_ng comprises a series of positive terms associated to a p-sequence of imposed monomial expectations of total even order p (2,4,6,8…). The higher the number of independent constraints, the higher the order of terms that are retained in that series and the more information is extracted from the joint PDF.

We have shown that the possible values of the cross-expectations lie within a bounded set obtained by Schwarz-like inequalities. We illustrate the range of I_ng values within those sets as function of third and fourth-order cross moments. Near the set’s boundary, I_ng tends to infinity where a deterministic relationship holds and the ME problem functional is ill-conditioned.

In order to better understand the possible sources of joint non-Gaussianity and non-Gaussian MI, we have used the preceding method for computing I_g and I_ng between the input and the output of a nonlinear channel contaminated by multiplicative and non-Gaussian noise for a full range of the signal-to-noise (snr) variance ratio. We put in evidence that sources of non-Gaussian MI arise from the nonlinearity of the transfer function, multiplicative noise and additive non-Gaussian noise [26].

Many of the results of the paper are straightforwardly generalized to the multivariate case with three or more random variables.

The paper is then organized as follows: Section 2 formalizes the Minimum Mutual Information (MinMI) principle from maximum entropy distributions, while Section 3 particularizes that principle to the MI decomposition into Gaussian and non-Gaussian MI parts. Section 4 addresses the non‑Gaussianity in a nonlinear non-Gaussian channel. The paper ends with conclusions and appendices with theoretical proofs and the numerical algorithm for solving the ME problem. This paper is followed by a companion one [27] on the estimation of non-Gaussian MI from finite samples with practical applications.

2. MI Estimation from Maximum Entropy PDFs

In this section we present the basis of the MI estimation between bivariate RVs, through the use of joint PDFs inferred by the maximum entropy (ME) method (ME-PDFs for short) in the space of transformed (anamorphed) marginals into specified ME distributions. We start with preliminary general concepts and definitions.

2.1. General Properties of Bivariate Mutual Information

Let $\left(X,Y\right)$ be continuous RVs with support given by the Cartesian product $S=S_X\otimes S_Y\subseteq {ℝ}^2$, and let the joint PDF ${\rho }_{X,Y}$ be taken absolutely continuous with respect to the product of the marginal PDFs ${\rho }_X$ and ${\rho }_Y$. The mutual information (in nat) between $X$ and $Y$ is a non-negative real functional of ${\rho }_{X,Y}$ expressed in different equivalent forms as:

$$I\left(X,Y\right)\equiv {\displaystyle \underset{S}{\mathrm{\int }}{\rho }_{X,Y}\left(x,y\right)\log \left(\frac{{\rho }_{X,Y}\left(x,y\right)}{{\rho }_X\left(x\right){\rho }_Y\left(y\right)}\right)dxdy}=D\left({\rho }_{XY}\left|\right|{\rho }_X{\rho }_Y\right)=H_{{\rho }_X}+H_{{\rho }_Y}-H_{{\rho }_{XY}}\ge 0,$$

(1)

in terms both of a KL divergence and Shannon entropies. The MI equals zero iff the RVs $X$ and $Y$ are statistically independent, i.e., iff ${\rho }_{X,Y}\left(x,y\right)={\rho }_X\left(x\right){\rho }_Y\left(y\right),\forall x\in S_X,y\in S_Y$, except possibly in a zero‑measure set in $S$. Under quite general regularity conditions of PDFs, statistical independence of $X$ and $Y$ is equivalent to the vanishing of the Pearson correlation between any pair of linear and/or nonlinear mapping functions $\tilde{X}(X)$ and $\tilde{Y}(Y)$. Then, linearly uncorrelated non-independent variables have at least a non-zero nonlinear correlation. The MI between transformed variables $\tilde{X}(X)$ and $\tilde{Y}(Y)$, differentiable almost everywhere, satisfies the ‘data processing inequality’ $I\left(\tilde{X}(X),\tilde{Y}(Y)\right)\le I\left(X,Y\right)\hspace{0.17em}$ [1] with equality occurring if both $\tilde{X}$ and $\tilde{Y}$ are smooth homeomorphisms of $X$ and $Y$ respectively.

2.2. Congruency between Information Moment Sets

• Definition 1: Following the notation of [28], we define the moment class ${\mathrm{\Omega }}_{T,\mathrm{\theta }}$ of bivariate PDFs ${\rho }_{XY}$ of $(X,Y)$, as:

  $${\mathrm{\Omega }}_{T,\mathrm{\theta }}=\left\{{\rho }_{XY}:E_{\rho }\left[T\right]=\mathrm{\theta }\right\}\hspace{0.17em}$$

  (2)

  where $E_{\rho }$ is the ρ-expectation operator, $T={\left(T_1,\mathrm{...},T_J\right)}^T$ is a J-dimensional vector composed of J absolutely $(X,Y)$ integrable functions with respect to ρXY, and $\mathrm{\theta }={\left({\mathrm{\theta }}_1,\mathrm{...},{\mathrm{\theta }}_J\right)}^T$ is the vector of function expectations of $T$. Here and henceforth, $\left(T,\mathrm{\theta }\right)$ is denoted as an information moment set.

The PDF $\rho \in {\mathrm{\Omega }}_{T,\mathrm{\theta }}$ that maximizes the Shannon entropy or the ME-PDF verifying the constraints associated to $\left(T,\mathrm{\theta }\right)$, hereby represented by ${\rho }_{T,\mathrm{\theta }}^{*}$, must exist since $H_{\rho }$ is a concave function in the non-empty convex set ${\mathrm{\Omega }}_{T,\mathrm{\theta }}$. The form and estimation of the ME-PDF is described in Appendix 1. The ME-PDF satisfies the following Lemma [1]:

• Lemma 1: Given two information encapsulated information moment sets: $\left(T,\mathrm{\theta }\right)\subseteq \left(T_1,{\mathrm{\theta }}_1\right)$ i.e., with $T_1$ including more constraining functions than $T$, the respective ME-PDFs ${\rho }_{{\rm T},\mathrm{\theta }}^{*},\hspace{0.17em}{\rho }_{T_1,{\mathrm{\theta }}_1}^{*}$, if they exist, satisfy the following conditions:

  $$E_{\rho }\left[-\log {\rho }_{T,\mathrm{\theta }}^{*}\right]=E_{{\rho }_{T1,\mathrm{\theta }1}^{*}}\left[-\log {\rho }_{T,\mathrm{\theta }}^{*}\right]=E_{{\rho }_{T,\mathrm{\theta }}^{*}}\left[-\log {\rho }_{T,\mathrm{\theta }}^{*}\right]=H_{{\rho }_{T,\mathrm{\theta }}^{*}}$$

  (3a)

  $$E_{\rho }\left[-\log {\rho }_{T1,\mathrm{\theta }1}^{*}\right]=E_{{\rho }_{T1,\mathrm{\theta }1}^{*}}\left[-\log {\rho }_{T1,\mathrm{\theta }1}^{*}\right]=H_{{\rho }_{T1,\mathrm{\theta }1}^{*}}$$

  (3b)

This means that any ME-PDF log-mean is unchanged when it is computed with respect to a more constraining ME-PDF. As a corollary we have $D\left({\rho }_{T_1,{\mathrm{\theta }}_1}^{*}\left|\right|{\rho }_{{\rm T},\mathrm{\theta }}^{*}\right)=H_{{\rho }_{T,\mathrm{\theta }}^{*}}-H_{{\rho }_{T_1,{\mathrm{\theta }}_1}^{*}}\ge 0$, i.e., the ME decreases with the increasing number of constraints. This can be understood bearing in mind that the entropy maximization is performed in a more reduced PDF class, since ${\mathrm{\Omega }}_{T_1,{\mathrm{\theta }}_1}\subseteq {\mathrm{\Omega }}_{T,\mathrm{\theta }}$.

• Definition 2: If two information moment sets $\left(T_1,{\mathrm{\theta }}_1\right)$, $\left(T_2,{\mathrm{\theta }}_2\right)$ are related by linear affine relationships, then the sets are referred to as “congruent”, a property hereby denoted as $\left(T_1,{\mathrm{\theta }}_1\right)\underset{PDF}{\cong }\left(T_2,{\mathrm{\theta }}_2\right)$, and consequently both PDF sets are equal, i.e., ${\mathrm{\Omega }}_{T_1,{\mathrm{\theta }}_1}={\mathrm{\Omega }}_{T_2,{\mathrm{\theta }}_2}$ [15]. A stronger condition than ‘congruency’ is the ME-congruency, denoted as $\left(T_1,{\mathrm{\theta }}_1\right)\underset{ME}{\cong }\left(T_2,{\mathrm{\theta }}_2\right)$, holding when both the associated ME-PDFs are equal. For example, both the univariate constraint sets $\left(T_1=X^2,{\mathrm{\theta }}_1=1\right)$ and $\left(T_2={(X^2,X^4)}^T,{\mathrm{\theta }}_2={(1,3)}^T\right)$ for $X\in ℝ$ lead to the same ME-PDF, the standard Gaussian N(0,1). Consequently both information moment sets are ME-congruent but not congruent since ${\mathrm{\Omega }}_{T_2,{\mathrm{\theta }}_2}\subset {\mathrm{\Omega }}_{T_1,{\mathrm{\theta }}_1}$. This is because the Lagrange multiplier of the ME functional (see Appendix 1) corresponding to the fourth moment is set to zero without any constraining effect. The congruency implies ME-congruency but not the converse.

2.3. MI Estimation from Maximum Entropy Anamorphoses

We are looking for a method of obtaining lower bound MI estimates from ME-PDFs. For that purpose we will decompose the information moment set as: $(T,\mathrm{\theta })=(T_{ind},{\mathrm{\theta }}_{ind})\cup (T_{cr},{\mathrm{\theta }}_{cr})$, say into a marginal or independent part $(T_{ind},{\mathrm{\theta }}_{ind})=(T_X,{\mathrm{\theta }}_X)\cup (T_Y,{\mathrm{\theta }}_Y)$ of single X and Y independent moments on $S_X,S_Y$, and a cross part $(T_{cr},{\mathrm{\theta }}_{cr})$ on $S_X\otimes S_Y$, made by moments of joint $\left(X,Y\right)$ functions not expressible as sums of single moments. For example, by taking $(T_X,{\mathrm{\theta }}_X)=\left((X,X^2),(0,1)\right)$, $(T_Y,{\mathrm{\theta }}_Y)=\left((Y,Y^2),(0,1)\right)$ and $(T_{cr},{\mathrm{\theta }}_{cr})=\left(XY,c\right)$ for $(X,Y)\in {ℝ}^2$ leads to a ME-PDF which is the bivariate Gaussian with correlation $c$ and standard Gaussians as marginal distributions.

The joint ME-PDF associated to the independent part $(T_{ind},{\mathrm{\theta }}_{ind})$ is the product of two independent ME-PDFs related to $(T_X,{\mathrm{\theta }}_X)$ and $(T_Y,{\mathrm{\theta }}_Y)$, with the joint entropy being the sum of marginal maximum entropies, as in the case of independent random variables [15]. The KL divergence between the ME-PDFs associated to $(T,\mathrm{\theta })$ and those associated to $(T_{ind},{\mathrm{\theta }}_{ind})$ is a proper MI lower bound or, expressed in other terms, a constrained mutual information. We denote it as:

$$I\left(X,Y:(T,\mathrm{\theta }),(T_{ind},{\mathrm{\theta }}_{ind})\right)\equiv D\left({\rho }_{T,\mathrm{\theta }}^{*}\left|\right|{\rho }_{{\rm T}ind,\mathrm{\theta }ind}^{*}\right)=H_{{\rho }_{Tind,\mathrm{\theta }ind}^{*}}-H_{{\rho }_{T,\mathrm{\theta }}^{*}}=H_{{\rho }_{TX,\mathrm{\theta }X}^{*}}+H_{{\rho }_{TY,\mathrm{\theta }Y}^{*}}-H_{{\rho }_{T,\mathrm{\theta }}^{*}}\ge 0$$

(4)

Its difference with respect to $I\left(X,Y\right)$ is given by:

$$I\left(X,Y\right)-I\left(X,Y:(T,\mathrm{\theta }),(T_{ind},{\mathrm{\theta }}_{ind})\right)=D\left({\rho }_{XY}\left|\right|{\rho }_{T,\mathrm{\theta }}^{*}\right)-\left[D\left({\rho }_X\left|\right|{\rho }_{{\rm T}X,\mathrm{\theta }X}^{*}\right)+D\left({\rho }_Y\left|\right|{\rho }_{{\rm T}Y,\mathrm{\theta }Y}^{*}\right)\right]\hspace{0.17em}$$

(5)

which can be negative. However, the positiveness is ensured when marginal distributions are set to the ME-PDFs constrained by ${\rho }_X(x)={\rho }_{{\rm T}X,\mathrm{\theta }X}^{*}(x),\forall x$ and${\rho }_Y(y)={\rho }_{{\rm T}Y,\mathrm{\theta }Y}^{*}(y),\forall y$, which results in KL divergences vanishing with respect to the marginal PDFs in Equation (5). This procedure is quite general because $(X,Y)$ can appropriately be obtained through an injective smooth maps $(X=X(\widehat{X}),Y=Y(\widehat{Y}))$ from an original pair of RVs $(\widehat{X},\widehat{Y})$preserving the MI, i.e., $I(\widehat{X},\widehat{Y})=I(X,Y)$. Those maps are monotonically growing homeomorphisms, the hereby called ME-anamorphoses, which are obtained by equaling mass probability functions of the original variable ${\rho }_{\widehat{X}}$ to those of the transformed variable ${\rho }_X$ (equally for ${\rho }_{\widehat{Y}}$ and ${\rho }_Y$) as:

$${\displaystyle {\mathrm{\int }}_{-\infty }^{X(\widehat{X})}{\rho }_{{\rm T}X,\mathrm{\theta }X}^{*}(u)du\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}}{\displaystyle {\mathrm{\int }}_{-\infty }^{\widehat{X}}{\rho }_{\widehat{X}}(u)du}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle {\mathrm{\int }}_{-\infty }^{Y(\widehat{Y})}{\rho }_{{\rm T}Y,\mathrm{\theta }Y}^{*}(u)du\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}}{\displaystyle {\mathrm{\int }}_{-\infty }^{\widehat{Y}}{\rho }_{\widehat{Y}}(u)du}$$

(6)

The moments in $(T_{ind},{\mathrm{\theta }}_{ind})$ are invariant under the transformation, $(\widehat{X},\widehat{Y})\to (X,Y)$, i.e., they are the same for both original and transformed variables: $E_{{\rho }_{\widehat{X}}}\left[T_X\right]=E_{{\rho }_X}\left[T_X\right]={\mathrm{\theta }}_X$ (idem for $Y$). Therefore, in the space of ME-anamorphed variables the ME-based MI (4) gives the minimum MI (MinMI), compatible with the cross moments $(T_{cr},{\mathrm{\theta }}_{cr})$.

Thanks to Lemma 1, a hierarchy of ME-based MI bounds is obtainable by considering successive supersets of the ME constraints on the ME-anamorphed variables, which is justified by the theorem below.

• Theorem 1: Let $(X,Y)$ be a pair of single random variables (RVs), distributed as the ME-PDF associated to the independent constraints $(T_{ind}=(T_X,T_Y),{\mathrm{\theta }}_{ind}=({\mathrm{\theta }}_X,{\mathrm{\theta }}_Y))$. Both variables can be obtained from previous ME-anamophosis. Let $(T_1,{\mathrm{\theta }}_1)=(T_{cr1}\cup T_{ind1},{\mathrm{\theta }}_{cr1}\cup {\mathrm{\theta }}_{ind1})$ be a subset of $(T_2,{\mathrm{\theta }}_2)=(T_{cr2}\cup T_{ind2},{\mathrm{\theta }}_{cr2}\cup {\mathrm{\theta }}_{ind2})$, i.e., $T_{cr1}\subseteq T_{cr2}$ and $T_{ind}\subseteq T_{ind1}\subseteq T_{ind2}$, such that all independent moment sets are ME-congruent (see Definition 2), i.e., $(T_{ind},{\mathrm{\theta }}_{ind})\underset{ME}{\cong }(T_{ind1},{\mathrm{\theta }}_{ind1})\underset{ME}{\cong }(T_{ind2},{\mathrm{\theta }}_{ind2})$, i.e., such that the independent extra moments in $(T_2,{\mathrm{\theta }}_2)$ are not further constraining the ME-PDF. Each marginal moment set is decomposed as $(T_{ind\hspace{0.17em}j}=(T_{Xj},T_{Yj}),{\mathrm{\theta }}_{ind\hspace{0.17em}j}=({\mathrm{\theta }}_{Xj},{\mathrm{\theta }}_{Yj}))\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}(j=1,2)$. For simplicity of notation, let us denote $I{\left(X,Y\right)}_j\equiv I\left(X;Y:(T_j,{\mathrm{\theta }}_j),(T_{ind\hspace{0.17em}j},{\mathrm{\theta }}_{ind\hspace{0.17em}j})\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(j=1,2)$. Then, the following inequalities between constrained mutual informations hold:

  $$I{\left(X,Y\right)}_j=I\left(X,Y\right)-D\left({\rho }_{XY}\left|\right|{\rho }_{T_j,{\mathrm{\theta }}_j}^{*}\right)\le I\left(X,Y\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(j=1,2)$$

  (7a)

  $$I{\left(X,Y\right)}_2-I{\left(X,Y\right)}_1=D\left({\rho }_{T_2,{\mathrm{\theta }}_2}^{*}\left|\right|{\rho }_{T_1,{\mathrm{\theta }}_1}^{*}\right)=H_{{\rho }_{T1,\mathrm{\theta }1}^{*}}-H_{{\rho }_{T2,\mathrm{\theta }2}^{*}}\ge 0$$

  (7b)

The proof is given in Appendix 2. The Theorem 1 justifies the possibility of building a monotonically growing sequence of lower bounds of $I\left(X,Y\right)$ from encapsulated sequences of cross-constraints $T_{cr\hspace{0.17em}1}\subseteq \mathrm{...}{T}_{cr\hspace{0.17em}j}\subseteq T_{cr\hspace{0.17em}j+1}\subseteq \mathrm{...}$ and independent constraints $T_{ind}\subseteq T_{ind\hspace{0.17em}1}\subseteq \mathrm{...}\subseteq T_{ind\hspace{0.17em}j}\subseteq T_{ind\hspace{0.17em}j+1}\subseteq \mathrm{...}$. In the sequence, the entropy associated to independent constraint sets is always constant due to the ME-congruency, while the entropy of the joint ME-PDF decreases, thus allowing the MinMIs to grow. Therefore, $I{\left(X,Y\right)}_j$ is the part of MI due to cross moments in $T_{cr\hspace{0.17em}j}$ and the positive difference $I{\left(X,Y\right)}_{j+1}-I{\left(X,Y\right)}_j$ is the increment of MI due to the additional cross moments in $T_{ind\hspace{0.17em}j+1}/T_{ind\hspace{0.17em}j}$, while marginals are kept as preset ME-PDFs (e.g., Gaussian, Gamma, Weibull).

3. MI Decomposition under Gaussian Marginals

3.1. Gaussian Anamorphosis and Gaussian Correlation

In this section we explain how to implement the sequence of MI estimators detailed in Section 2.3 for the particular case where $X$ and $Y$ are standard Gaussian RVs. Our aim is to estimate the MI between two original variables $(\widehat{X},\widehat{Y})$ of null mean and unit variance with real support. Those variables are then transformed through a homeomorphism, the Gaussian anamorphosis [14], into standard Gaussian RVs, respectively $X~N(0,1)$ and $Y~N(0,1)$, given by:

$$X(\widehat{x})={\mathrm{\Phi }}^{-1}\left({\displaystyle {\mathrm{\int }}_{-\infty }^{\widehat{x}}{\rho }_{\widehat{X}}(u)du}\right)=G_{\widehat{X}}(\widehat{x})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Y(\widehat{y})={\mathrm{\Phi }}^{-1}\left({\displaystyle {\mathrm{\int }}_{-\infty }^{\widehat{y}}{\rho }_{\widehat{Y}}(u)du}\right)=G_{\widehat{Y}}(\widehat{y})$$

(8)

where Φ is the mass distribution function for the standard Gaussian. If $(\widehat{X},\widehat{Y})$ are marginally non‑Gaussian, then the Gaussian anamorphoses are nonlinear transformations. In practice, Gaussian anamorphoses can be approximated empirically from finite data sets by equaling cumulated histograms. However, for certain cases, it is analytically possible to construct bivariate distributions with specific marginal distributions and the knowledge of the joint cumulative distribution function [29].

In the case of Gaussian anamorphosis, the information moment set $(T_{ind},{\mathrm{\theta }}_{ind})$ of Theorem 1 includes the first and second independent moments of each variable: $E[X]=E[Y]=0;E[X^2]=E[Y^2]=1$. Then, following the proposed procedure of Section 2.3, we will consider a sequence of cross-constraint sets for determining the hierarchy of lower MI bounds.

The most obvious cross moment to be considered is the $XY$ expectation, equal to the Gaussian correlation $c_g\equiv cor(G_{\widehat{X}},G_{\widehat{Y}})=cor(X,Y)$ between the ‘Gaussianized’ variables $(X,Y)$. The difference between $c_g$and the linear correlation $c=\text{cor}(\widehat{X},\widehat{Y})$ is easily expressed as:

$$c_g-c=E\left[\left(G_{\widehat{X}}-\widehat{X}\right)\left(G_{\widehat{Y}}-\widehat{Y}\right)\right]$$

(9)

The signal of the factor $G_{\widehat{X}}-\widehat{X}$ in Equation (9) roughly depends on the skewness $sk(\widehat{X})=E[{\widehat{X}}^3]$ and excess of kurtosis $kur(\widehat{X})=E[{\widehat{X}}^4]-3$ through the rule of thumb, stating that $\mathrm{sgn}(G_{\widehat{X}}-\widehat{X})$ is approximated by $-\mathrm{sgn}(sk(\widehat{X}))$ and $-\mathrm{sgn}(kur(\widehat{X}))\hspace{0.17em}\mathrm{sgn}(\widehat{X})$, respectively for a skewed $\widehat{X}$ PDF and a symmetric $\widehat{X}$ PDF (idem for $\widehat{Y}$). Therefore, $c_g$ can result in an enhancement of correlation $c$ or in the opposite effect, as shown in [21] for the RV pair of meteorological variables ($\widehat{X}$ = North Atlantic Oscillation Index, $\widehat{Y}$ = monthly precipitation). The Gaussian correlation is a concordance measure like the rank correlation and Kendall τ, being thus invariant for a monotonically growing smooth homeomorphism of both $\widehat{X}$and $\widehat{Y}$. Those measures are expressed as functionals of the bivariate copula-function $c[u={\displaystyle {\mathrm{\int }}_{-\infty }^{X_1}{\rho }_{X_1}(x)dx},v={\displaystyle {\mathrm{\int }}_{-\infty }^{X_2}{\rho }_{X_2}(y)dy}]={\rho }_{X_1{X}_2}(X_1,X_2)/\left({\rho }_{X_1}(X_1){\rho }_{X_2}(X_2)\right)$, which is uniquely dependent on the cumulated marginal probabilities and equal to the density ratio, independently from the specific forms of marginal PDFs [24]. In particular, the Gaussian correlation is given by:

$$c_g={\displaystyle {\mathrm{\int }}_0^1{\displaystyle {\mathrm{\int }}_0^{1}c[u,v]{\mathrm{\Phi }}^{-1}(u)}{\mathrm{\Phi }}^{-1}(v)du\hspace{0.17em}dv}$$

(10)

3.2. Gaussian and Non-Gaussian MI

The purpose of this sub-section is to express which part of MI comes from joint non-Gaussianity. If the ‘Gaussianized’ variables $(X,Y)$ are jointly non-Gaussian, then the original standardized variables $(\widehat{X},\widehat{Y})$, obtained from $(X,Y)$ by invertible smooth monotonic transformations, are jointly non‑Gaussian as well. However, the converse is not true. In fact, any nonlinear transformation $(\widehat{X},\widehat{Y})$ of jointly Gaussian RVs $(X,Y)$ leads to non-Gaussian $(\widehat{X},\widehat{Y})$ with the joint non-Gaussianity arising in a trivial way. Therefore, the ‘genuine’ joint non-Gaussianity can only be diagnosed in the space of Gaussianized variables $(X,Y)$. When the marginal PDFs of $(\widehat{X},\widehat{Y})$ are non-Gaussian then $c_g\ne c$ in general. In particular, if the correlation $c$ is null as it occurs when $(\widehat{X},\widehat{Y})$ are principal components (PCs), then $c_g$ can be non-null, leading to statistically dependent PCs since they are non-linearly correlated.

The MI between Gaussianized variables $(X,Y)$ or $(\widehat{X},\widehat{Y})$ is expressed as $I(\widehat{X},\widehat{Y})=I(X,Y)=2H_g-H_{{\rho }_{XY}}$, where $H_g\equiv 1/2\log (2\pi e)$ is the Shannon entropy of the standard Gaussians $X$, $Y$ and $H_{{\rho }_{XY}}$ is the $(X,Y)$joint entropy, associated to the joint PDF ${\rho }_{XY}$. Given the above equality, the MI is decomposed into two positive terms: $I(\widehat{X},\widehat{Y})=I{(\widehat{X},\widehat{Y})}_g+I{(\widehat{X},\widehat{Y})}_{ng}$. The first term is the Gaussian MI [21] given by $I{(\widehat{X},\widehat{Y})}_g=I{(X,Y)}_g=-1/2\log (1-c_g^2)\equiv I_g(c_g)\ge 0$, as a function of the Gaussian correlation $c_g$ (see its graphic in Figure 1), and the second one is the non-Gaussian MI $I{(\widehat{X},\widehat{Y})}_{ng}=I{(X,Y)}_{ng}$, which is due to joint non-Gaussianity and nonlinear statistical relationships among variables.

The MI $I(\widehat{X},\widehat{Y})$ is related to the negentropy $J(\widehat{X},\widehat{Y})$, i.e., to the KL divergence between the PDF and the Gaussian PDF with the same moments of order one and two. That is shown by:

• Theorem 2: Given ${({\widehat{X}}_r,{\widehat{Y}}_r)}^T=A{(\widehat{X},\widehat{Y})}^T$, a pair of rotated standardized variables (A being an invertible 2 × 2 matrix), one has the following result with proof in Appendix 2:

  $$\begin{array}{l}J(\widehat{X},\widehat{Y})=J(\widehat{X})+J(\widehat{Y})+I(\widehat{X},\widehat{Y})-I_g(cor(\widehat{X},\widehat{Y}))=\\ =J({\widehat{X}}_r,{\widehat{Y}}_r)=J({\widehat{X}}_r)+J({\widehat{Y}}_r)+I({\widehat{X}}_r,{\widehat{Y}}_r)-I_g(cor({\widehat{X}}_r,{\widehat{Y}}_r))\end{array}$$

  (11)

A simple consequence is that in the space of uncorrelated variables (i.e., $I_g(cor(\widehat{X},\widehat{Y}))=0$), the joint negentropy is the sum of marginal negentropies with the MI, thus showing that there are intrinsic and joint sources of non-Gaussianity. One interesting corollary is derived from that.

• Corollary 1: For standard Gaussian variables $(X,Y)$ and standardized rotated ones${(X_r,Y_r)}^T=A{(X,Y)}^T$, we have

  $$I_{ng}(X,Y)=J(X,Y)=J(X_r,Y_r)=J(X_r)+J(Y_r)+I(X_r,Y_r)-I_g(cor(X_r,Y_r))$$

  (12)

For the proof it suffices to consider Gaussian variables ${(\widehat{X},\widehat{Y})}^T={(X,Y)}^T$ in (11). Their self-negentropy vanishes by definition and the correlation term is the Gaussian MI.

Negentropy has the property of being invariant for any orthogonal or oblique rotations of the Gaussianized variables $(X,Y)$. However, this invariance does not extend to $I{(X,Y)}_{ng}$. From (12), in particular when $(X_r,Y_r)$ are uncorrelated (e.g., standardized principal components of $(X,Y)$ or uncorrelated standardized linear regression residuals), the negentropy equals the KL divergence between the joint PDF and that of an isotropic Gaussian with the same total variance. That KL divergence is the compactness (level of concentration around a lower-dimensional manifold), as defined in [25]. This measure is invariant under orthogonal rotations. The last term of (12) vanishes due the fact that the null correlation allows one to decompose of $I{(X,Y)}_{ng}$ into positive contributions: the self-negentropies of uncorrelated variables $(X_r,Y_r)$ and their MI $I(X_r,Y_r)=I{(X_r,Y_r)}_g+I{(X_r,Y_r)}_{ng}$. These variables can be ‘Gaussianized’ and rotated, leading to further decomposition of $I{(X_r,Y_r)}_{ng}$ until the possible “emptying”/depletion of the initial joint non-Gaussianity into Gaussian MIs and univariate negentropies. The PDF of the new rotated variables will be closer to an isotropic spherical Gaussian PDF. Since it is algorithmically easier to compute univariate rather than multivariate entropies, the above method can be used for an efficient estimation of MIs.

The search for rotated variables maximizing the sum of individual negentropies $J(X_r)+J(Y_r)$ in (12) with minimization of $I(X_r,Y_r)$ or their statistical dependency is the goal of Independent Component Analysis (ICA) [4].

A natural generalization of the MI decomposition is possible when $(X,Y)$ is obtained from a generic ME-anamorphosis by decomposing the MI into a term associated to correlation, under the constraint that marginals are set to given ME-PDFs (the equivalent to I_g), and to a term not explained by correlation (the equivalent to I_ng). There is however no guarantee that this decomposition is unique as in the case of non-Gaussians, since there is no natural bivariate extension of univariate prescribed PDFs with a given correlation [30].

By looking again at Equation (12), we notice that when original variables are correlated $(c=cor(X_r,Y_r)\ne 0)$, the sum of marginal negentropies is not necessarily a lower bound of the joint negentropy $J(\widehat{X},\widehat{Y})$, because in some cases $I(\widehat{X},\widehat{Y})-I_g(c)$ can be negative. This means that $I_g(c)$ is not generally a proper lower bound of the MI. An example of that is given by the following discrete distribution with support on four points: $(\widehat{X},\widehat{Y})=(1,1),(1,-1),(-1,1),(-1,-1)$, with mass probabilities $P_{\widehat{X}\widehat{Y}}$, respectively of (1 + c)/4, (1 + c)/4, (1 − c)/4 and (1 − c)/4. This 4-point distribution has $\widehat{X}$ and $\widehat{Y}$ zero means, unit variances and Pearson correlation $c$. The PDF is made of four Dirac-Deltas. In this case, the mutual information $I(\widehat{X},\widehat{Y})$ is easily computed from the 4-point discrete mean of $\log (P_{\widehat{X}\widehat{Y}}/(P_{\widehat{X}}{P}_{\widehat{Y}}))$ while the marginal mass probabilities are $P_{\widehat{X}}(1)=P_{\widehat{Y}}(1)=P_{\widehat{X}}(-1)=P_{\widehat{Y}}(-1)=1/2$. After some lines of algebra the MI becomes:

$$\hspace{0.17em}I(X,Y)=\log \left[2{\left(\frac{1+\left|c\right|}{2}\right)}^{\left(\frac{1+\left|c\right|}{2}\right)}{\left(\frac{1-\left|c\right|}{2}\right)}^{\left(\frac{1-\left|c\right|}{2}\right)}\right]\equiv I_{d4}(c)\hspace{0.17em}$$

(13)

The MI is bounded for any value of c and $I_{d4}(c)\to \log (2)<\infty$ as $(|c|\hspace{0.17em}\to 1)$. The graphic of $I_{d4}(c)\hspace{0.17em}$is included in Figure 1, showing that $I_g(c)>I_{d4}(c)$ for about $\left|c\right|>0.45$. This behavior is also reproduced by finite discontinuous PDFs (e.g., replacing the Dirac-Deltas by cylinders of probability) as well as continuous PDFs. We test it by approximating the discrete PDF by the weighted superposition of 4 spherical bivariate Gaussian PDFs, all of which have a sharp isotropic standard deviation σ = 0.001 and are centered at the 4 referred centroids. Figure 1 depicts successively growing MI lower bounds for each value of c, using the maximum entropy method (labels, $I_g$, $I_g+I_{ng,4}$, $I_g+I_{ng,6}$, $I_g+I_{ng,8}$, see Section 3.3 for details), showing the convergence to the asymptotic value log(2).

Figure 1. Semi-logarithmic graphs of $I_g(c)$ (black thick line), $I_4(c)$, (grey thick line) and of the successive growing estimates of $I_4(c)$: $I_g$, $I_g+I_{ng,4}$, $I_g+I_{ng,6}$ and $I_g+I_{ng,8}$ (grey thin lines). See text for details.

Figure 1. Semi-logarithmic graphs of $I_g(c)$ (black thick line), $I_4(c)$, (grey thick line) and of the successive growing estimates of $I_4(c)$: $I_g$, $I_g+I_{ng,4}$, $I_g+I_{ng,6}$ and $I_g+I_{ng,8}$ (grey thin lines). See text for details.

3.3. The Sequence of Non-Gaussian MI Lower Bounds from Cross-Constraints

In order to build a monotonically increasing sequence of lower bounds for $I(X,Y)$ (cf. Section 2.3), we have considered a sequence of encapsulated sets of functions whose moments will constrain the ME-PDFs. Those functions consist of single (univariate) and combined (bivariate) monomials of the standard Gaussians X and Y. The numerical implementation of joint ME-PDFs constrained by polynomials in dimensions d = 2, 3 and 4 was studied by Abramov [17,18,19], with particular emphasis on the efficiency and convergence of iterative algorithms. Here, we use the algorithm proposed in [21] and explained in the Appendix 1. Let us define the information moment set $T_p$ as the set of bivariate monomials with total order up to p:

$$T_p\equiv \left(X^i{Y}^j:\hspace{0.17em}\hspace{0.17em}1\le i+j\le p,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(i,j)\in N_0^2\right)\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}p\in N$$

(14)

This set is decomposed into marginal (independent) and cross monomials as $T_p\equiv \left(T_{p\hspace{0.17em}\hspace{0.17em}ind},T_{p\hspace{0.17em}\hspace{0.17em}cr}\right)$ with the corresponding vector of expectations ${\mathrm{\theta }}_p\equiv ({\mathrm{\theta }}_{p\hspace{0.17em}\hspace{0.17em}ind},{\mathrm{\theta }}_{p\hspace{0.17em}\hspace{0.17em}cr})=(E[T_{p\hspace{0.17em}\hspace{0.17em}ind}],E[T_{p\hspace{0.17em}\hspace{0.17em}cr}])$ referring respectively to 2p and p(p − 1)/2 independent and cross moments. The components of ${\mathrm{\theta }}_p$ are of the form$m_{i,j}\equiv E[X^i{Y}^j],1\le i+j\le p$, where the independent part is set to the moment values of the standard Gaussian: $E_g[X^p]=0$ for an odd p and $E_g[X^p]=(p-1)!!=(p-1)(p-3)\mathrm{...}$ for an even p (idem for Y), where $E_g$denotes expectation with respect to the standard Gaussian. The cross moments are bounded by a simple application of the Schwartz inequality as ${(m_{i,j})}^2\le E_g[X^{2i}]E_g[Y^{2j}]$. The monomials of $T_p$ are linearly independent and space-generating by linear combinations, thus forming a basis in the sense of vector spaces.

For an even integer p and constraint set $T_p$ there exists an integrable ME-PDF with support in $ℝ$, of the form ${\rho }^{*}(X,Y)=C\exp (P_p(X,Y))$, where $P_p(X,Y)$ is a p^th order polynomial given by a linear combination of monomials in $T_p$ with weights given by Lagrange multipliers and C being a normalizing constant. That happens because $|(X^i{Y}^j)|\le O(X^p+Y^p),\forall X,Y\in ℝ,1\le i+j\le p,$thus allowing for $P_p(X,Y)\to -\infty$ as $\left|X\right|,\left|Y\right|\to \infty$, leading to convergence of ME-integrals [18]. For odd p, those integrals diverge because the dominant monomials $X^p,Y^p$ change sign from positive to negative real and therefore the ME-PDF is not well defined.

In order to build a sequence of MI lower bounds and use the procedure of Section 2.3, we have considered the sequence of information moment sets of even order $(T_2,{\mathrm{\theta }}_2),(T_4,{\mathrm{\theta }}_4),(T_6,{\mathrm{\theta }}_6),\mathrm{...}\hspace{0.17em}$ with any pair of consecutive sets satisfying the premises of Theorem 1, i.e., all independent moment sets are ME-congruent. This will lead to the corresponding monotonically growing sequence of lower bounds of MI, denoted as $I{(X,Y)}_{g,2}\le I{(X,Y)}_{g,4}\le I{(X,Y)}_{g,6}\le \mathrm{...}$, where the subscript g means that variables are marginally standard Gaussian. The first term of the sequence is the Gaussian MI $I_{g,2}(X,Y)\equiv I{(X,Y)}_g$, dependent upon the Gaussian correlation. The difference between the subsequent terms and the first one leads to the non-Gaussian MI of order p defined as $I{\left(X,Y\right)}_{ng,p}\equiv I{\left(X,Y\right)}_{g,p}-I{\left(X,Y\right)}_{g,2}\le I{\left(X,Y\right)}_{ng}$, which increases with p and converges to $I{\left(X,Y\right)}_{ng}$ as $p\to \infty$, under quite general conditions [15].

In the same manner as stated in (12), the lower bound $I{\left(X,Y\right)}_{ng,p}$ of $I{\left(X,Y\right)}_{ng}$ is also a lower bound for the joint negentropy, which is invariant for any affine transformations ${(X_r,Y_r)}^T=A{(X,Y)}^T$, i.e.,

$$I{\left(X,Y\right)}_{ng,p}=H{(X,Y)}_2-H{(X,Y)}_p=H{(X_r,Y_r)}_2-H{(X_r,Y_r)}_p$$

(15)

where $H{(X,Y)}_p$ is the bivariate ME associated to $(T_p,{\mathrm{\theta }}_p)$ and $H{(X,Y)}_2=2H_g-I{(X,Y)}_g$. The successive differences $I{\left(X,Y\right)}_{ng,p+2}-I{\left(X,Y\right)}_{ng,p}$ are non-negative, with the extra MI being explained by moments of orders p + 1 and p + 2, while not explained by lower-order moments.

There is no analytical closed formula for the dependence of non-Gaussian MI on cross moments. However, under the scenario of low joint non-Gaussianity (small KL divergence to the joint Gaussian), the ME-PDF can be approximated by the Edgeworth expansion [31], based on orthogonal Hermite polynomials and _ng approximated as a polynomial of joint bivariate cumulants: $k^{[i,j]}\equiv E[X_r^i{Y}_r^j]-E_g[X_r^i]E_g[Y_r^j]$, i + j > 2, of any pair of uncorrelated standardized variables ${(X_r,Y_r)}^T=A{(X,Y)}^T$. Cross-cumulants are nonlinear correlations measuring joint non-Gaussianity [32], vanishing when $\left(X,Y\right)$ are jointly Gaussian. For example $I{\left(X,Y\right)}_{ng,p=4}$ is approximated as in [13] by the sum of squares:

$$\begin{array}{l}I{(X,Y)}_{ng(Ed,p=4)}\equiv \frac{1}{12}\left[3{\left(k^{[2,1]}\right)}^2+3{\left(k^{[1,2]}\right)}^2+{\left(k^{[3,0]}\right)}^2+{\left(k^{[0,3]}\right)}^2\right]\\ +\frac{1}{48}\left[{\left(k^{[4,0]}\right)}^2+4{\left(k^{[3,1]}\right)}^2+6{\left(k^{[2,2]}\right)}^2+4{\left(k^{[1,3]}\right)}^2+{\left(k^{[0,4]}\right)}^2\right]+\\ \frac{1}{72}[5{\left(k^{[3,0]}\right)}^4+30{\left(k^{[3,0]}{k}^{[2,1]}\right)}^2+\frac{2}{3}{\left(9k^{[3,0]}{k}^{[1,2]}+6k^{[2,1]}{k}^{[2,1]}\right)}^2+\frac{2}{3}{\left(9k^{[0,3]}{k}^{[2,1]}+6k^{[1,2]}{k}^{[1,2]}\right)}^2+\\ \frac{1}{2}{\left(2k^{[3,0]}{k}^{[0,3]}+18k^{[2,1]}{k}^{[1,2]}\right)}^2+30{\left(k^{[1,2]}{k}^{[0,3]}\right)}^2+5{\left(k^{[0,3]}\right)}^4]=I{(X,Y)}_{ng}+O\left({n_{eq}}^{-3/2}\right)\end{array}$$

(16)

where $(X,Y)$ is assumed to be the arithmetic average of an equivalent number n_eq of independent and identically distributed (iid) bivariate RVs. Therefore, from the multidimensional Central Limit Theorem [33], the larger n_eq is, the closer the distribution is to joint Gaussianity, and the smaller the absolute value of cumulants become.

3.4. Non-Gaussian MI across the Polytope of Cross Moments

The $(X,Y)$ cross moments in the expectation vector ${\mathrm{\theta }}_p$ (p even) are not completely free. Rather, they satisfy to Schwarz-like inequalities defining a compact set $D_p$ within which cross-moments lie. That set portrays all the possible non-Gaussian ME-PDFs with p-order independent moments equal to those of the standard Gaussian. Under these conditions $I{\left(X,Y\right)}_{ng}=I{\left(X,Y\right)}_{ng,p}$. In order to have a better feel on how _ng behaves, we have numerically evaluated $I_{ng,p=4}$ along the allowed set of cross-moments.

In order to determine that set, let us begin by invoking some generalities about polynomials. Any bivariate polynomial $P_p(x,y):{ℝ}^2\to \hspace{0.17em}ℝ\hspace{0.17em}$ of total order p is expressed as a linear combination of linearly independent monomials from the basis $T_p^1\equiv T_p\cup \left\{1\right\}$, obtained from $T_p$ including unity. Then:

$$P_p(x,y)={\displaystyle {\sum }_i{a}_i{T}_{i,p}(x,y)}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{T}_{i,p}\in T_p^1$$

(17)

If the condition $P_p(x,y)\ge 0,\hspace{0.17em}\forall (x,y)\in {ℝ}^2$ holds i.e., $P_p(x,y)$ is positive semi-definite (PSD), then its expectation is non-negative: $E[P_p(x,y)]={\displaystyle {\sum }_i{a}_i{\mathrm{\theta }}_{i,p}\ge 0}$, where $\hspace{0.17em}{\mathrm{\theta }}_{i,p}\equiv E[T_{i,p}]$, thus imposing a constraint on the components of ${\mathrm{\theta }}_p$. A sufficient condition for the positiveness is that $P_p(x,y)$ is a sum of squares (SOS) or, without loss of generality, the square of a certain polynomial $Q_{p/2}(x,y)=b^T{T}_{p/2}^1$ of total order p/2, where $b$ is a column vector of coefficients multiplying monomials of $T_{p/2}^1$. Then, $P_p(x,y)$ is written as a quadratic form $P_p(x,y)={\left(Q_{p/2}(x,y)\right)}^2=b^T\left(T_{p/2}^1{T_{p/2}^1}^T\right)b\ge 0$. By taking the expectation operator we have $b^{T}E[\left(T_{p/2}^1{T_{p/2}^1}^T\right)]b\ge 0,\hspace{0.17em}\hspace{0.17em}\forall b$, which implies the positiveness of the matrix of moments $E[\left(T_{p/2}^1{T_{p/2}^1}^T\right)]\equiv C_p$, which is given in terms of components of ${\mathrm{\theta }}_p$.

When p = 4 and d = 2, the case of bivariate quartics, any PSD polynomial is a SOS [34] and vice versa. However, for p ≥ 6 there are PSD-non-SOS polynomials (e.g., those coming from the inequality between arithmetic and geometric means [35]). Therefore, a necessary and sufficient condition among fourth-order moments is that $E[\left(T_2^1{T_2^1}^T\right)]\equiv C_4$ be a PSD matrix. Let us study the conditions for that.

By ordering the set $T_2^1\equiv {\left(1,x,y,x^2,y^2,xy\right)}^T$, one has the 6 × 6 matrix $C_4$, written in the simplified form in terms of moments:

$$C_4=\left[\begin{array}{l}(1,0,0,1,1,c_g)\\ (0,1,c_g,0,m_{1,2},m_{2,1})\\ (0,c_{g}1,m_{2,1},0,m_{1,2})\\ (1,0,m_{2,1}3,m_{2,2},m_{3,1})\\ (1,m_{1,2}0,m_{2,2}3,m_{1,3})\\ (c_g,m_{2,1},m_{1,2},m_{3,1},m_{1,3},m_{2,2})\end{array}\right]$$

(18)

A necessary and sufficient condition for the positiveness of $C_4$ is given by the application of the Sylvester criterion, stating that the determinants d1, d2, d3, d4, d5 and d6 of the 6 upper sub-matrices of $C_4$ are positive. From these, only those of orders 4, 5 and 6 lead to nontrivial relationships, given with help of Mathematica® [36] as:

$$d_4\equiv 2-2c_g^2-{m_{2,1}}^2 \ge 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_{4,dual}\equiv 2-2c_g^2-{m_{1,2}}^2 \ge 0$$

(19)

$$d_5\equiv -2 \left(m_{2,1}^2+m_{2,1}^2\right)+m_{1,2}^2 m_{2,1}^2-2 m_{1,2}{m}_{2,1} \left(m_{2,2}-1\right) c_g+\left(m_{2,2}-3\right) \left(m_{2,2}+1\right) \left(-1+c_g^2\right)\ge 0$$

(20)

$$\begin{array}{l}{d}_6\equiv \left(9-m_{2,2}^2\right)c_g^4+2\left[\left(m_{1,3}+m_{3,1}\right)\left(m_{2,2}-3\right)+m_{1,2} m_{2,1} m_{2,2}\right]c_g^3+\\ \left[\begin{array}{l}\left(m_{2,1}^2+m_{1,2}^2\right)\left(9-2m_{2,2}\right)-m_{1,2} m_{2,1}\left(2\left(m_{1,3}+m_{3,1}\right)+m_{1,2} m_{2,1}\right)+\\ 2\left(m_{3,1}^2+m_{1,3}^2+m_{1,3} m_{3,1}\left(1-m_{2,2}\right)\right)+\left(m_{2,2}-3\right)\left(3+m_{2,2}\left(2+m_{2,2}\right)\right)\end{array}\right]c_g^2+\\ 2\left[\begin{array}{l}\left(3-2m_{2,1}^2-2m_{1,2}^2-m_{2,2}\right)\left(m_{1,3}+m_{3,1}\right)+m_{2,2}\left(m_{1,2}^2{m}_{3,1}+m_{2,1}^2{m}_{1,3}\right)+\\ m_{1,2} m_{2,1}\left(m_{1,3} m_{3,1}+m_{2,1}^2+m_{1,2}^2-3+m_{2,2}\left(5-2m_{2,2}\right)\right)\end{array}\right]c_g+\\ \left[\begin{array}{l}-m_{2,2}^3+\left(m_{2,1}^2+m_{1,2}^2+2\right)m_{2,2}^2+\left(3-4\left(m_{2,1}^2+m_{1,2}^2\right)+m_{1,2} m_{2,1}\left(3m_{1,2} m_{2,1}-2\left(m_{1,3}+m_{3,1}\right)\right)\right)m_{2,2}+\\ m_{1,2} m_{2,1}\left(-2\left(m_{1,2}^2{m}_{3,1}+m_{2,1}^2{m}_{1,3}+m_{1,2} m_{2,1}\right)+6\left(m_{1,3}+m_{3,1}\right)\right)+2\left(m_{2,1}^4+m_{1,2}^4\right)+\\ -3\left(m_{2,1}^2+m_{1,2}^2\right)-2\left(m_{1,3}^2+m_{3,1}^2\right)-2m_{1,3} m_{3,1}+{\left(m_{1,3} m_{2,1}\right)}^2+{\left(m_{3,1} m_{1,2}\right)}^2\end{array}\right]\ge 0\end{array}$$

(21)

In Equation (22), the inequality for $d_4$ has a dual relationship (d_4,dual), its sign being reversed by swapping the two indices in m_i,j, whereas $d_5$ and $d_6$ are symmetric with respect to indicial swap. The term d₆ of Equation (21) is a fourth-order polynomial of c_g. The inequalities for$d_4$,$d_{4\hspace{0.17em}dual}$, $d_5$ and $d_6$ hold inside a compact domain denoted D₄, with the shape of what resembles a rounded polytope in the space of cross moments ($m_{1,2}, m_{2,1}, m_{1,3}, m_{3,1}, m_{2,2}$), for each value of the Gaussian correlation cg. The case of Gaussianity lies within the interior of D₄, corresponding to $m_{1,2}=m_{2,1}=0; m_{1,3}=m_{3,1}=3c_g\text{; }{m}_{2,2}=2c_g^2+1$, thus defining the hereafter called one-dimensional ‘Gaussian manifold’.

In order to illustrate how non-Gaussianity depends on moments of third and fourth order, we have computed the non-Gaussian MI of order 4 ($I_{ng,p=4}$), along a set of 2-dimensional cross-sections of D₄ crossing the Gaussian manifold and extending up the boundary of D₄. For Gaussian moments, $I_{ng,4}$ vanishes, being approximated by the Edgeworth expansion (16) near the Gaussian manifold.

In order to get a picture of $I_{ng,p=4}$, we have chosen six particular cross-sections of D₄ by varying two moments and setting the remaining to their ‘Gaussian values’. The six pairs of varying parameters are: A (c_g, m_2,1), B (c_g, m_3,1), C (c_g, m_2,2), D (m_2,1, m_3,1) at c_g=0, E (m_2,1, m_2,2) at c_g=0 and F (m_3,1, m_2,2) at c_g = 0, with the contours of the corresponding I_ng field shown in Figure 2a–f. The fields are retrieved from a discrete mesh of 100 × 100 in moment space. The Gaussian state lies at: (c_g, 0), (c_g, 3c_g), (c_g, 2c_g² + 1), (0, 0), (0, 1) and (0, 1), respectively for cases A up to F. The moment domains are obtained by solving inequalities for d₄, d₅ and d₆ and applying the restrictions imposed by the crossing of the Gaussian manifold (e.g., m_1,2 = 0, m_1,3 = m_3,1 = 3c_g, m₂,₂ = 2c_g² + 1 for case A). We obtain the following restrictions for cases A–F:

$$\text{Case} \mathrm{A}:\hspace{0.17em}{m}_{2,1}^2\le \frac{1}{2}\left[3+3c_g^2\left(1-2c_g^2\right)-\left(1-c_g^2\right){\left(1+52c_g^2+28c_g^4\right)}^{1/2}\right]$$

(22)

$$\text{Case} \mathrm{B}:\hspace{0.17em}{c}_g+2c_g^3-\sqrt{2\left(1-c_g^2-c_g^4+c_g^6\right)}\le m_{3,1}\le c_g+2c_g^3+\sqrt{2\left(1-c_g^2-c_g^4+c_g^6\right)}$$

(23)

$$\text{Case} \mathrm{C}:\hspace{0.17em}\frac{1}{2}\left[-1+c_g^2+{\left(1+34c_g^2+c_g^4\right)}^{1/2}\right]\le m_{2,2}\le 3$$

(24)

$$\text{Case} \mathrm{D}:\hspace{0.17em}{m}_{3.1}^2\le 2-3m_{2,1}^2+4m_{2,1}^4\hspace{0.17em};\hspace{0.17em}\left|m_{21}\right|\le 1\hspace{0.17em};\hspace{0.17em}\left|m_{21}\right|\le \sqrt{2}$$

(25)

$$\text{Case} \mathrm{E}:\hspace{0.17em}{m}_{2.1}^2\le m_{2,2}\le 1+{\left(4-2m_{2,1}^2\right)}^{1/2}\hspace{0.17em};\hspace{0.17em}{m}_{2,1}^2\le \sqrt{3}$$

(26)

$$\text{Case} \mathrm{F}: \hspace{0.17em}{m}_{2,2}^3-2m_{2,2}^2-3m_{2,2}+2m_{3,1}^2\le 0$$

(27)

The analytical boundary of allowed domains is emphasized with thick solid lines in all Figure 2a–f. There are some common aspects among the figures. As expected, I_ng vanishes at the Gaussian states, the Gaussian manifold, marked with G in Figure 2. I_ng grows monotonically towards the boundary of the moment domains D₄. There, $\det (C_4)=0$, meaning that $C_4$ is singular and there is a vector $(b\ne 0)\in Ker(C_4)$. This holds if one gets the deterministic relationship $Q_{p/2}(x,y)=b^T{T}_2^1=0$, leading to a Dirac-Delta-like ME-PDF along a one-dimensional curve. This in turn leads to I_ng = ∞, except possibly in a set of singular points of D₄ on which I_ng is not well defined. In practice, infinity is not reached due to stopping criteria for convergence of the iterative method used for obtaining the ME‑PDF.

Figure 2. Field of the non-Gaussian MI $I_{ng,4}$ along 6 bivariate cross-sections of the set of allowed moments. Two varying moments are featured in each cross-section: A (c_g, m_2,1), B (c_g, m3_,1), C (c_g, m_2,2), D (m_2,1, m_3,1) at c_g=0, E (m_2,1, m_2,2) at c_g = 0 and F (m_3,1, m_2,2) at c_g = 0. The letter G indicates points or curves where Gaussianity holds.

Figure 2. Field of the non-Gaussian MI $I_{ng,4}$ along 6 bivariate cross-sections of the set of allowed moments. Two varying moments are featured in each cross-section: A (c_g, m_2,1), B (c_g, m3_,1), C (c_g, m_2,2), D (m_2,1, m_3,1) at c_g=0, E (m_2,1, m_2,2) at c_g = 0 and F (m_3,1, m_2,2) at c_g = 0. The letter G indicates points or curves where Gaussianity holds.

At states where |c_g| = 1, I_g = ∞ and I_ng has a second-kind singularity discontinuity where the contours merge together without a well-defined limit for I_ng. In the neighborhood of the Gaussian state with c_g = 0 in Figure 2d–f, I_ng is approximated by the quadratic form (16) as is confirmed by the elliptic shape of I_ng contours. The value of I_ng can surpass I_g, thus emphasizing the fact that in some cases much of the MI may come from nonlinear (X,Y) correlations.

The joint entropy is invariant for a mirror symmetry in one or both variables: $X\to -X$ or $Y\to -Y$, because the absolute value of the determinant of that transformation equals 1. As a consequence, the dependence of the Gaussian and non-Gaussian MI on moments also reflects these intrinsic mirror symmetries. For instance, in Figure 2d, the symmetry $X\to -X$ leads to the dependency relations I_ng(m_2,1, m_3,1) = I_ng(m_2,1, −m_3,1), where arguments are the varying moments, while symmetry $Y\to -Y$ leads to I_ng(m_2,1, m_3,1) = I_ng(−m_2,1, −m_3,1). The symmetries corresponding to the remaining figures are: I_ng(c_g, m_2,1) = I_ng(−c_g, m_2,1) = I_ng(c_g, −m_2,1); I_ng(c_g, m_3,1) = I_ng(−c_g, −m_3,1); I_ng(c_g, m_2,2) = I_ng(−c_g, m₂₂); I_ng(m_2,1, m_2,2) = I_ng(−m_2,1, m_2,2); I_ng(m_3,1, m_2,2) = I_ng(−m_3,1, m_2,2), respectively for Figure 2a–c, e and f.

Near the boundary of $D_4$, the ME problem is very ill-conditioned. The Hessian matrix of the corresponding ME functional is the inverse $M_4^{-1}$, where $M_4\equiv E[T_4{T}_4^T]-E[T_4]E[T_4^T]$ is the covariance matrix of $T_4$, made by 8 independent plus 6 cross moments. The condition number CN (ratio between the largest and smallest eigenvalue) of $M_4$ (and of $M_4^{-1}$) tends to ∞ at the boundary of $D_4$ where $M_4$ is singular, due to the above deterministic relationship. Therefore, the closer that boundary is, the closer the ME-PDF is to a deterministic relationship, the more ill conditioned the ME problem is and the slower the numerical convergence of the optimization algorithm becomes.

4. The Effect of Noise and Nonlinearity on Non-Gaussian MI

The aim of this section is an exploratory analysis of the possible sources of non-Gaussianity in a bivariate statistical relationship. Towards that aim, we explore the qualitative behavior of I_g and I_ng between a standardized signal $\widehat{X}$(with null mean and unit variance) and an $\widehat{X}$-dependent standardized response variable $\widehat{Y}$ contaminated by noise. For this purpose, a full range of signal-to-noise variance ratios (snr) shall be considered, from pure signal to pure noise. The statistics are evaluated from one‑million-long synthetic $(\widehat{X},\widehat{Y})$ iid independent realizations produced by a numeric Gaussian random generator. Many interpretations are possible for the output variable: (i) $\widehat{Y}$ taken as the observable outcome emerging from a noisy transmission channel fed by $\widehat{X}$ (ii) given by the direct or indirect observation affected by measurement and representativeness errors corresponding to a certain value $\widehat{X}$ of the model state vector [37] (iii) the outcome from a stochastic or deterministic dynamical system [38].

In order to estimate $I(\widehat{X},\widehat{Y})$, the working variables $(\widehat{X},\widehat{Y})$ are transformed by anamorphosis into standard Gaussian variables $(X,Y)$.

We consider, without loss of generality, $X=\widehat{X}$. The variable $Y$ is given by Gaussian anamorphosis $Y=G_{\widehat{Y}}(\widehat{Y})~N(0,1)\hspace{0.17em}$ as in Equation (8), with:

$$\widehat{Y}=s\hspace{0.17em}F(X)+\hspace{0.17em}{(1-s^2)}^{1/2}\hspace{0.17em}n(X,W)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}s\in [0,1]$$

(28)

where $F(X)$ is a purely deterministic transfer function and $n(X,W)$ is a scalar noise uncorrelated with $F(X)$, depending in general on $X$ (e.g., multiplicative noise) and from a vector $W$ of independent Gaussians contaminating the signal. Both $F(X)$ and $n(X,W)$ have unit variance with $n(X,0)=0$. The signal-to-noise variance ratio is $snr=s^2/(1-s^2)$.

Then, the Gaussian MI I_g is computed for each value of s∈[0,1] and compared among several scenarios of $F(X)$ and $n(X,W)$. A similar comparison is done for the non-Gaussian MI, approximated here by I_ng,p=8. Six case studies have been considered (A, B, C, D, E and F); their signal and noise terms are summarized in Table 1, along with the colors with which they are represented in Figure 3 further below.

Table 1. Types of signal and noise in Equation (35) and corresponding colors used in Figure 3.

Case study	$F(X)$	$n(X,W)$	Color
A—Gaussian noise (reference)	$X$	$W$	Black
B—Additive non-Gaussian independent noise	$X$	$W^3/\sqrt{15}$	Red
C—Multiplicative noise	$W\hspace{0.17em}X$	Blue
D—Smooth nonlinear homeomorphism	$X^3/\sqrt{15}$	$W$	Magenta
E—Smooth non-injective transformation	$(X^3-X)/\sqrt{10}$	Green
F—Combined non-Gaussianity	$(X^3-X)/\sqrt{10}$	$X{W}_1/\sqrt{2}+W_2^3/\sqrt{30}$	Cyan

Table 1. Types of signal and noise in Equation (35) and corresponding colors used in Figure 3.

Case study	$F(X)$	$n(X,W)$	Color
A—Gaussian noise (reference)	$X$	$W$	Black
B—Additive non-Gaussian independent noise	$X$	$W^3/\sqrt{15}$	Red
C—Multiplicative noise	$W\hspace{0.17em}X$	Blue
D—Smooth nonlinear homeomorphism	$X^3/\sqrt{15}$	$W$	Magenta
E—Smooth non-injective transformation	$(X^3-X)/\sqrt{10}$	Green
F—Combined non-Gaussianity	$(X^3-X)/\sqrt{10}$	$X{W}_1/\sqrt{2}+W_2^3/\sqrt{30}$	Cyan

In Table 1, $W,W_1,W_2$ are independent standard Gaussian noises. We begin with a reference case, A, which refers to Gaussian noise. Case B refers to a symmetric leptokurtic (i.e., with kurtosis larger than that of the Gaussian) non-Gaussian noise. In case C the multiplicative noise depends linearly on the signal X. In case D, the signal is a nonlinear cubic homeomorphism of the real domain. For case E, the signal is nonlinear and not injective in the interval [−1, 1], thus introducing ambiguity in the relationship $(X,Y)$. Finally, in case F all the factors—non-Gaussian noise, multiplicative noise and signal ambiguity—are pooled together.

Figure 3. Graphs depicting the total MI (a), Gaussian MI (b) and non-Gaussian MI (c) of order 8 for 6 cases (A–F) of different signal-noise combinations with the signal weight s in abscissas varying from 0 up to 1. See text and Table 1 for details about the cases and their color code.

Figure 3. Graphs depicting the total MI (a), Gaussian MI (b) and non-Gaussian MI (c) of order 8 for 6 cases (A–F) of different signal-noise combinations with the signal weight s in abscissas varying from 0 up to 1. See text and Table 1 for details about the cases and their color code.

Figure 3a,b show the graphics of the total MI, estimated by I_g + I_ng,8 and of the Gaussian MI (I_g) for the six cases (A to F). The graphic of non-Gaussian MI as approximated by I_ng,8 is depicted in Figure 3c for five cases (B to F). In Figure 4, we show a ‘stamp-format’ collection of the contouring of ME-PDFs of polynomial order p = 8 for all cases (A to F) and extreme and intermediate cases of the snr: s = 0.1, s = 0.5 and s = 0.9. This illustrates how the snr and the nature of both the transfer function and noises influence the PDFs.

For the Gaussian noise case (A), the non-Gaussian MI is theoretically set to zero since the joint distribution of $(X,Y)$ is Gaussian. In all scenarios, both I_g and the total MI I_g + I_ng grow, as expected, with the snr. This is in accordance to the Bruijn’s equality stating the positiveness of the MI derivative with respect to snr and established in the literature of signal processing for certain types of noise [39,40]. On the contrary, the monotonic behavior as a function of snr is not a universal characteristic of the non-Gaussian MI.

By observing Figure 3a–c, the following qualitative results are worth mentioning. We begin by comparing the total MI in three cases (A, B and C), which share the same linear signal but feature noise of different kinds (Figure 3a). Both the red (B) and blue (C) lines lie above the black line (A) for each given s, thus indicating that the total MI is lowest when the noise is Gaussian. This means that the Gaussian noise is the most signal degrading of noises with the same variance [41]. The extra MI found in the B and C cases come, respectively, from the Gaussian MI (see case B in Figure 3c) and from the non-Gaussian MI (see case C in Figure 3a), as it is also apparent by looking at ME-PDFs for cases B, C (s = 0.1) (Figure 4).

We consider now the cases A, D, E, all of which have a Gaussian noise. Their differences lie in the signals, with the one in A being linear and the ones in D and E being nonlinear. By comparing these cases it is seen that I_g is highest for the linear signal, the black curve lying above the magenta (D) and green (E) curves for each s in Figure 3b. This indicates that the Gaussian MI, measuring the degree of signal linearity, is lower when the signal introduces nonlinearity (cases D and E) than when no nonlinearity is present (case A).

Figure 4. Collection of stamp-format ME-PDFs for cases A–F (see text for details) and signal weight s = 0.1 (a), s = 0.5 (b) and s = 0.9 (c) over the [−3, 3]² support set.

Figure 4. Collection of stamp-format ME-PDFs for cases A–F (see text for details) and signal weight s = 0.1 (a), s = 0.5 (b) and s = 0.9 (c) over the [−3, 3]² support set.

It is worth noting that, while the signals in A and D are injective, the one in E is not, thus introducing ambiguity. This will imply loss of information in E, which is visible in the total MI depicted for each s in Figure 3a. In fact, there the green curve (case E) lies lower than the black (A) and magenta (D) curves for every s. The effect of nonlinearity is quite evident in ME-PDFs, in particular for high s value (Figure 4, cases D, E, s = 0.9).

We focus now on the non-Gaussian MI, depicted in Figure 3c for each s. The curve representing the case B, with a linear signal and a state-independent noise, indicates that the non-Gaussian MI is null for both s = 0 and s = 1. The first zero of non-Gaussian MI (at s = 0) is justified by the noise being state‑independent, whereas the second zero (at s = 1) is due to the signal being linear, which means that all the MI resides in the Gaussian MI. The non-Gaussian MI is thus positive and maximum at intermediate values of s.

By looking at case C (multiplicative noise), it is seen that the non-Gaussian MI remains roughly unchanged for every s < 1. This holds even at s = 0 (pure noise), since the noise is state-dependent and thus some information is already present. At s = 1 the non-Gaussian MI is null due to the signal being linear (as in case B).

By observing the cases with Gaussian noise and nonlinear signals (D and E) in Figure 3c, it can be seen that their non-Gaussian MI grows with s (and thus with the relative weight of the signal), due to their signals being nonlinear. This gradual behavior is also reflected in the ME-PDFs (Figure 4, cases D, E along the s values).

Finally, we consider the case in which the signal is nonlinear and the noise comprises a multiplicative and a non-Gaussian additive component (case F). As compared with E (which differs from it in that the noise is Gaussian), it can been that non-Gaussian MI is always larger in F independently of s. This is due to the fact that in F there is information even at s = 0, due to the state-dependence of its noise. For all values of s, the ME-PDF exhibits quite a large deviation from Gaussianity.

7. Discussion and Conclusions

We have addressed the problem of finding the minimum mutual information (MinMI), or the least noncommittal MI between d = 2 random variables, consistent with a set of marginal and joint expectations. The MinMI is a proper MI lower bound when marginals are set to ME-PDFs through appropriate nonlinear single anamorphoses. Moreover, the MinMI increases as long as one increases number of independent cross-constraints of the bivariate ME problem. Considering a sequence of moments, we have obtained a hierarchy of lower MI bounds approximating the total MI value. The method can easily be generalized for d > 2 variables with the necessary adaptations.

One straightforward application of that principle follows from the MI estimation from ‘Gaussianized’ variables with real support, where the marginals are rendered standard Gaussian N(0,1) by Gaussian anamorphosis. This allows for the MI decomposition into two positive contributions: a Gaussian term I_g, which depends uniquely on the Gaussian correlation c_g (Pearson correlation in the space of ‘Gaussianized’ variables), and a non-Gaussian term I_ng depending on nonlinear correlations. This term is equal to the joint negentropy, which is invariant for any oblique or orthogonal rotation of the ‘Gaussianized’ variables and is related to the ‘compactness’ measure or the closeness of the PDF towards a low manifold deterministic relationship. The Gaussian MI is also a ‘concordance’ measure, invariant for any monotonically growing homeomorphisms of marginals and consequently expressed as a functional of the copula-density function, which is exclusively dependent on marginal cumulated probabilities. In certain extreme cases, very far from Gaussianity, the Pearson correlation among non‑Gaussian variables is not a proper measure of the mutual information. An example of that situation is given.

Cross moments under marginal standard Gaussians are bounded by Schwarz-like inequalities defining compact sets, the shape of which resemble a rounded polytope where cross moments live. The allowed moment values portray all possible joint PDFs with Gaussian marginals. Inside that set lies the so called one-dimensional Gaussian manifold, parametrized by c_g, where joint Gaussinity holds. There, I_ng vanishes, growing towards infinity as far as the boundary is approached, where variables satisfy a deterministic relationship and the ME problem is ill conditioned. This behavior is illustrated in cross-sections of the polytope of cross moments of total order p = 4.

In order to systematize the possible sources of Gaussian and non-Gaussian MI, we have computed it in the context of nonlinear noisy channels. The MI has been computed between a Gaussian input and a panoply of (linear and/or nonlinear) outputs contaminated by different kinds of noise for a full range of the signal-to-noise variance ratio. Sources of non-Gaussian MI include: (a) the nonlinearity of the signal transfer function, (b) multiplicative noise and (c) non-Gaussian additive noise. This paper is followed by a companion one [27] on the estimation of non-Gaussian MI from finite samples with practical applications.