Why variational inference and not maximum likelihood?

When using variational inference scheme, we assume latent variable $$\mathbf z$$, model $$p(\mathbf x, \mathbf z)$$, and maximize $$\log p(\mathbf x)$$.

Introducing variational distribution $$q(\mathbf z)$$, the loss (negative ELBO) can be then written as: $$-\mathop{\mathbb{E}}[\log p(\mathbf x)] + \text{KL}[q(\mathbf z) || p(\mathbf z)]$$. In amortized inference we model $$q(\mathbf z | \mathbf x)$$, which together with the reparametrization trick gives the so-called variational autoencoder model (VAE).

My question is, why not start from the fact that we want to maximize $$\log p(\mathbf x)$$ and write $$\log p(\mathbf x) = \int \log p(\mathbf x | \mathbf z) p(\mathbf z) d\mathbf z = \mathop{\mathbb{E}}_{p(\mathbf z)} [\log p(\mathbf x | \mathbf z)]$$ Now we can just optimize this with Monte Carlo by sampling from $$p(\mathbf z)$$.

Is the issue that we don't really get the true posterior $$p(\mathbf z | \mathbf x)$$? If so, why do we really care about the posterior? Distribution $$q(\mathbf z)$$ is anyways an approximation to the true posterior (often not really good), and ELBO is just the lower bound on the function we are trying to maximize.

Or is it the case that we need the second term always as a regularization? In that case, why do we want to approximate posterior $$p(\mathbf z | \mathbf x)$$, and not put regularization on $$p(\mathbf z)$$. As I see it, everything is intractable so there is no much use of strict Bayesian behavior when we approximate on multiple levels.

An example is GMM model. We can just set parameters (means, variances...) as learnable and learn with gradient based optimization from data.

The equation you use is not correct:

$$\log p(\mathbf x) \color{red}{\neq} \int \log p(\mathbf x | \mathbf z) p(\mathbf z) d\mathbf z = \mathop{\mathbb{E}}_{p(\mathbf z)} [\log p(\mathbf x | \mathbf z)].$$

$$\log p(\mathbf x) = \log \int p(\mathbf x | \mathbf z) p(\mathbf z) d\mathbf z = \log \mathop{\mathbb{E}}_{p(\mathbf z)} [p(\mathbf x | \mathbf z)].$$

You could now sample from within the log, but then you would not get an unbiased estimate.

A further option is to use Jensen's identity to get a lower bound and optimise that. And that's exactly what variational inference does.