The authors proposed to use a neural network to find an encoder function $q_\psi(x)$ to approximate the true posterior $p_{\theta^}(z|x)$, and a decoder function to approximate the true generative model $p_{\theta^}(x|z)$ (data distribution)
For any probability density function of the observed data $x$, we can maximize the corresponding log-likelihood with our input data. The log-likelihood of our input data can be decomposed into:
$\text{log}p(x)=\text{log}p(x)\int_zq(z|x)dz$
$=\int_zq(z|x)\text{log}p(x)dz$
$=\int_zq(z|x)\text{log}(\frac{p(x, z)}{p(z|x)})dz$
$=\int_zq(z|x)\text{log}(\frac{p(x, z)}{q(z|x)}\frac{q(z|x)}{p(z|x)})dz$
$=\int_zq(z|x)\text{log}(\frac{p(x, z)}{q(z|x)})dz+\int_zq(z|x)\text{log}(\frac{q(z|x)}{p(z|x)})dz$
, where $q(z|x)$ is the posterior distribution of $z$ given the input data $x$.
Let the first term on the right-hand side be:
$L_b=\int_zq(z|x)\text{log}(\frac{p(x, z)}{q(z|x)})dz$
, and the second term on the right hand side would be:
$D_{KL}(q(z|x)\;||\;p(z|x)) = \int_zq(z|x)\text{log}\frac{q(z|x)}{p(z|x)}dz$
Since Kullback–Leibler divergence (KL divergence) is always $\geq 0$ (proof):
$\text{log}p(x)=L_b+D_{KL}(q(z|x)\;||\;p(z|x))\geq L_b+0=L_b$
To maximize the log-likelihood of $p(x)$, we try to maximize the lower bound $L_b$, this is also known as the evidence lower bound (ELBO) of log-likelihood of $p(x)$.