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

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: 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.

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 }}$.

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:

Its difference with respect to $I\left(X,Y\right)$ is given by: 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:

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.

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).

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: 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:

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:

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.

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:

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).

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:

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., 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: 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.

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:

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:

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: