This is a late response but I hope it may help:
Proof (this proof is basically a summary of the explanations from the author):
To go for the proof, first we can follow the procedure given by the author which allows us to have a more convenient expression:
$$\begin{aligned}
\theta_{ML} &=\arg \max_\theta p_{model}(\mathbb{X};\theta)\\
&= \arg\max_\theta\prod_{i=1}^m p_{model}(x^{(i)};\theta) \\
&= \arg\max_\theta\sum_{i=1}^m \log(p_{model}(x^{(i)};\theta)) \\
&= \arg\max_\theta \frac{1}{m}\sum_{i=1}^m \log(p_{model}(x^{(i)};\theta))\\
&= \arg\max_\theta \sum_{i=1}^m \frac{1}{m} \log(p_{model}(x^{(i)};\theta))
\end{aligned}$$
Note this last expression is the expectation of this $\log(p_{model}(x;\theta))$ function with respect to the empirical distribution defined by the training data ($\hat{p}_{data}$) which puts a probability of $1/m$ on each of the $m$ points $x^{(1)},x^{(2)},...,x^{(m)}$. So we can equivalently write this last expression as:
$$ \theta_{ML} = \arg\max_{\theta}\mathbb{E}_{x\sim \hat{p}_{data}} \log(p_{model}(x;\theta))$$
Hence, obtaining the value of $\theta$ which satisfies this expresion will maximize the likelihood of $p_{model}(x,\theta)$ being the statistical model that best fits our set of data samples $\mathbb{X}$.
But we can also get the parameter $\theta$ that maximizes the likelihood using the KL divergence of the probability distributions $\hat{p}_{data}$ (empirical distribution of $\mathbb{X}$) and $p_{model}$ (our statistical model that we are using to fit $\mathbb{X}$):
$$ D_{\text{KL}}(\hat{p}_{data} \parallel p_{model}) = \mathbb{E}_{x\sim \hat{p}_{data}} [
\log(\hat{p}_{data}(x)) - \log(p_{model}(x;\theta))]$$
This is because of the explanation given by the author, which is that $\mathbb{E}_{x\sim \hat{p}_{data}}
\log(\hat{p}_{data}(x))$ does not depend on $\theta$ (it only depends on the data generating process), so it can be trated as a constant. Hence we can adress the same problem of finding the value of $\theta$ that maximizes the likelihood by minimizing this KL divergence, because this is the same as minimizing:
$$ \mathbb{E}_{x\sim \hat{p}_{data}} [- \log(p_{model}(x;\theta))]$$
Which is just the negative form of the simplified expression for $\theta_{ML}$ that we have written earlier from the perspective of the likelihood.
So, to sum up, we can also calculate the parameter $\theta_{ML}$ by:
$$
\theta_{ML} = \arg\min_{\theta} D_{\text{KL}}(\hat{p}_{data} \parallel p_{model})= \arg\min_{\theta} \mathbb{E}_{x\sim \hat{p}_{data}} [- \log(p_{model}(x;\theta))]
$$
Intuition:
With this said, I believe we can think of using the KL divergence for maximizing the likelihood as a way of making the predicted distribution $(p_{model}(\mathbb{X},\theta))$ as close as possible to the empirical distribution.
Thereby with our predicted distribution and by sampling it, we would be able to obtain a set of samples similar to the initial ones ($\mathbb{X}$). So this may mean that we have correctly calculate the true distribution of the data $\mathbb{X}$.
