Optimal way to compute pairwise mutual information using numpy

Question

For an m x n matrix, what's the optimal (fastest) way to compute the mutual information for all pairs of columns (n x n)? By mutual information, I mean: I(X, Y) = H(X) + H(Y) - H(X,Y) where H(X) refers to the Shannon entropy of X. Currently I'm using `np.histogram2d` and `np.histogram` to calculate the joint (X,Y) and individual (X or Y) counts. For a given matrix `A` (e.g. a 250000 X 1000 matrix of floats), I am doing a nested `for` loop, ``` n = A.shape[1] for ix = arange(n) for jx = arange(ix+1,n): matMI[ix,jx]= calc_MI(A[:,ix],A[:,jx]) ``` Surely there must be better/faster ways to do this? As an aside, I've also looked for mapping functions on columns (column-wise or row-wise operations) on arrays, but haven't found a good general answer yet. Here is my full implementation, following the conventions in the Wiki page: ``` import numpy as np def calc_MI(X,Y,bins): c_XY = np.histogram2d(X,Y,bins)[0] c_X = np.histogram(X,bins)[0] c_Y = np.histogram(Y,bins)[0] H_X = shan_entropy(c_X)…

Accepted Answer

I can't suggest a faster calculation for the outer loop over the n*(n-1)/2 vectors, but your implementation of `calc_MI(x, y, bins)` can be simplified if you can use scipy version 0.13 or scikit-learn.

In scipy 0.13, the `lambda_` argument was added to `scipy.stats.chi2_contingency` This argument controls the statistic that is computed by the function. If you use `lambda_="log-likelihood"` (or `lambda_=0`), the log-likelihood ratio is returned. This is also often called the G or G2 statistic. Other than a factor of 2*n (where n is the total number of samples in the contingency table), this is the mutual information. So you could implement `calc_MI` as:

```
from scipy.stats import chi2_contingency

def calc_MI(x, y, bins):
    c_xy = np.histogram2d(x, y, bins)[0]
    g, p, dof, expected = chi2_contingency(c_xy, lambda_="log-likelihood")
    mi = 0.5 * g / c_xy.sum()
    return mi
```

The only difference between this and your implementation is that this implementation uses the natural logarithm instead of the base-2 logarithm (so it is expressing the information in "nats" instead of "bits"). If you really prefer bits, just divide `mi` by log(2).

If you have (or can install) `sklearn` (i.e. scikit-learn), you can use `sklearn.metrics.mutual_info_score`, and implement `calc_MI` as:

```
from sklearn.metrics import mutual_info_score

def calc_MI(x, y, bins):
    c_xy = np.histogram2d(x, y, bins)[0]
    mi = mutual_info_score(None, None, contingency=c_xy)
    return mi
```

Optimal way to compute pairwise mutual information using numpy

Solution

Problem