Is Variational Inference to make an inference of conditional probability from joint probability? While studying variational inference, I am told that it is to make an inference from conditional probability from joint probability, directly. Somehow, it is kind of hard to understand why it is so, and I am also told that the reason why we do this is that sometimes it is very hard to do integrals. What can be the case for it?
 A: In Bayesian inference, you postulate a statistical model over some observable quantities $\mathbf{x}=(x_1,\dots,x_n)$ and some unobservable quantities $\mathbf{z}=(z_1,\dots,z_m)$ (often model parameters, but possibly also latent variables). The model usually takes the form of a likelihood function $p(\mathbf{x}|\mathbf{z})$ and a prior distribution over the latent variables $p(\mathbf{z})$. The goal of the inference is to use data over the observable quantities $\mathbf{x}$, to infer the posterior distribution $p(\mathbf{z}|\mathbf{x})$ of the latent variables, given observed data (Bayes' theorem):
$$p(\mathbf{z}|\mathbf{x})=\frac{p(\mathbf{x}|\mathbf{z})p(\mathbf{z})}{p(\mathbf{x})} $$
The denominator is the marginal density $p(\mathbf{x})$ of the observables, and it's often called the evidence. Computing this marginal density requires computing a multidimensional integral in the $m$ unobservable variables $z_1,\dots,z_m$ :
$$p(\mathbf{x})=\int{p(\mathbf{x}|\mathbf{z})p(\mathbf{z})d\mathbf{z}}$$
The computation of a multidimensional integral can be a real hassle: for example, if you use a tensored Gaussian quadrature rule, the number of quadrature points grows exponentially with $m$, which precludes this approach for high-dimensional approaches. There are of course more efficient multidimensional quadrature rules, but even these have shortcomings and limitations. Sometimes you don't even have an explicit expression for the integrand, because you don't have an explicit formula for the likelihood. If you are at least able to evaluate it for arbitrary values of the parameters, then you can either use quadrature rules (again, if the problem has low or moderate dimensionality), or you can use methods which estimate the posterior without having to compute the evidence, such as MCMC. However, these approaches, even if less computationally intensive than numerical quadrature in high dimensions, are still quite intensive. Also, judging convergence of an MCMC sampler is not simple.
An alternative, deterministic approach exists, which doesn't compute $p(\mathbf{z}|\mathbf{x})$ "exactly", but it gives an analytical approximation to it, which is often surprisingly good. This is Variational Inference or Variational Bayes. The concept is very simple: we consider a family $\mathcal{D}$ of distributions $q(\mathbf{z})$ which are simpler than the posterior, but which could be good approximations to it. The goal of Variational Inference is to find $q^* \in \mathcal{D}$ such that $q^*(\mathbf{z})$ is the best approximation to $p(\mathbf{z}|\mathbf{x})$, among all distributions in $\mathcal{D}$. "Closeness" is usually measured by the Kullback–Leibler divergence*:
$$q^*(\mathbf{z}) = \arg\!\min_{q\in\mathcal{D}} \text{KL}(q||p)$$
Note that we talk about "closeness" and not distance, because KL-divergence is not a distance, or metric (it's not symmetric). Now, we can prove that
$$\text{KL}(q||p)=\mathbb{E}[\log{q(\mathbf{z})}]-\mathbb{E}[\log p(\mathbf{z}|\mathbf{x})]$$
where expectations are with respect to $q(\mathbf{z})$. Computing $$\mathbb{E}[\log p(\mathbf{z}|\mathbf{x})]$$ would still require to compute the evidence and thus it would be useless. For this reason, we use the definition of conditional distribution and we get
$$\text{KL}(q||p)=\mathbb{E}[\log{q(\mathbf{z})}]-\mathbb{E}[\log p(\mathbf{z},\mathbf{x})]+\mathbb{E}[\log p(\mathbf{x})]=\mathbb{E}[\log{q(\mathbf{z})}]-\mathbb{E}[\log p(\mathbf{z},\mathbf{x})]+\log p(\mathbf{x})$$
where in the last passage we used the fact that $p(\mathbf{x})$ is simply a constant. Since $p(\mathbf{x})$ is a constant, minimizing $\text{KL}(q||p)$ correponds to maximizing the quantity
$$\text{ELBO}(q)=\mathbb{E}[\log p(\mathbf{z},\mathbf{x})]-\mathbb{E}[\log{q(\mathbf{z})}]$$
$\text{ELBO}(q)$ is called the evidence lower bound for a very good reason: since $\text{KL}(q||p)\geq 0 \ \forall p,q$ and since $\log p(\mathbf{x})= \text{KL}(q||p) + \text{ELBO}(q)$, of course $\log p(\mathbf{x})\geq\text{ELBO}(q)$. Note that $\text{ELBO}(q)$ is a functional, i.e., is a function which maps a function $q(\mathbf{z})$ to a real number. We use Variational Calculus to maximize or minimize functionals, and that's why this approach to Bayesian inference is called Variational Bayes.
Thus, finally we reduced our problem to maximizing $\text{ELBO}(q)$ over all distributions belonging to $\mathcal{D}$, and in doing so we got rid of the complicated integral which would be needed to compute the evidence, which answers your latter point (why Variational Inference allows us not to compute "hard" integrals). Concerning your former point, I think it's imprecise: we don't approximate a conditional distribution with a joint one (which wouldn't make much sense), because the  approximation to $p(\mathbf{z}|\mathbf{x})$ is $q^*(\mathbf{z})$ which is a distribution over the $\mathbf{z}$ (like $p(\mathbf{x}|\mathbf{z})$), not a joint distribution over $\mathbf{z}$ and $\mathbf{x}$. However, it's true that in $\text{ELBO}(q)$ (the functional we want to maximize), the joint distribution $p(\mathbf{z},\mathbf{x})$ appears instead than the posterior $p(\mathbf{z}|\mathbf{x})$ (which is a conditional distribution). Thus maybe this could be the sense in which we "went" from a conditional to a joint distribution...(?)
Finally, very good, open access references on this topic are:
Variational Inference: A Review for Statisticians
A Tutorial on Variational Bayesian Inference
The Wikipedia page on Variational Bayes

* in an earlier version of this question I claimed that when we use KL divergence to measure "closeness", then We obtain the so-called mean-field Variational Bayes approach. This is wrong: mean-field VB does indeed choose KL as a measure of closeness, but it also postulates that the latent variables are independent, so that
$$q(\mathbf{z}) =\prod_{i=1}^m q(z_i)$$
This simplifies inference at the cost of accuracy. Forcing $q(\mathbf{z})$ to belong to a full-factorized family of distributions will usually lead to a worse approximation to $p(\mathbf{z}|\mathbf{x})$, unless of course the latent are actually mutually independent conditionally on observed data.
