How to compute joint entropy of high-dimensional data? Normally, I compute the (empirical) joint entropy of some data, using the following code:
import numpy as np

def entropy(x):
    counts = np.histogramdd(x)[0]
    dist = counts / np.sum(counts)
    logs = np.log2(np.where(dist > 0, dist, 1))
    return -np.sum(dist * logs)

x = np.random.rand(1000, 5)
h = entropy(x)

This works perfectly fine, as long as the number of features does not get too large (histogramdd can maximally handle 32 dimensions, i.e. features). Assuming that I would like to compute the joint entropy $H(X_1, X_2, \ldots, X_{728})$ of the MNIST dataset, is it possible to compute this? If yes, how can this be done?
 A: Generally, estimating the entropy in high-dimensions is going to be intractable. What you can do instead is estimate an upper bound on the entropy. 
Note that entropy can be written as an expectation:
$$H(X_1, \ldots, X_n) = -\mathbb E_p \log p(x)$$
Here, $\mathbb E_p$ is an expectation over the distribution $p(x)$. 
Imagine that you fit some other generative model, $q(x)$, that you can calculate exactly. This won't be exactly the same as $p(x)$ but it can help you get a upper bound on the entropy of $p(x)$.
The KL divergence can be written as:
$$D(p(x)\| q(x)) = \mathbb E_p \log p(x) - \mathbb E_p \log q(x)$$
Using Jensen's inequality, we can see that KL divergence is always non-negative, and therefore, $H(X) = -\mathbb E_p \log p(x) \leq - \mathbb E_p \log q(x)$. The quantity on the right is what people sometimes call the negative log-likelihood of the data (drawn from $p(x)$) under the model, $q(x)$.  Because $D(p(x)\| p(x)) = 0$ and $D(p(x)\| q(x)) \geq 0$, this implies that no model, $q$, can give a better score for negative log likelihood than the true distribution, $p$. 
The negative log likelihood is often reported in papers as a measure of how well you have modeled the data, here's one example (see Table 1) that links to others.
Intuitively, why can't we exactly calculate the entropy, or provide nearly tight lower bounds? Entropy measures the optimal compression for the data. If you know the true entropy, you are saying that the data can be compressed this much and not a bit more. That's difficult in high-dimensions because there could always be some hidden structure that could help you compress a little more but that you might not observe with a small number of samples. 
EDIT: I forgot one really important component from your question. You are estimating entropy by binning your data. This is definitely going to fail in high dimensions. Suppose you have 2 bins for each dimension (maybe greater or less than 0.5). Then in $d=784$ dimensions, the total number of bins is $2^{784}$. With only $60,000$ samples, almost every bin will be empty. You don't have enough to samples empirically estimate the frequency of each bin. That's why papers like the one I linked use more sophisticated strategies for modeling $q(x)$ that have a small number of parameters that can be estimated more reliably.
