In the lecture by D.Blei: https://www.cs.princeton.edu/courses/archive/fall11/cos597C/lectures/variational-inference-i.pdf
Variational inference is explained and he shows how to derive the optimal variational distribution. However, he does not show the mathematical details on how he did to go from equation (22) to equation (24).
Could someone help me understand in detail how that happened or point me to a good resource? 
From (22) to (23), it involves functional derivatives. The same question can be found here.
From (23) to (24), it just set (23) equals to 0 and get the corresponding $q_j(z_j)$. The constant in the exponential function is left out, and the constrain that $\sum_k q^*(z_k)=1$ is omitted, so the proportional sign is used. To get the complete $q^*(z_k)$, the normalizing term should be divided: $q^*(z_k) = \frac{q^*(z_k)}{\sum_k q^*(z_k)}$
Instead of using the gradients/derivative, the update equation can be derived using the KL divengence.
Given the equation (21) in the question, the ELBO can be written as a function of $q_j$ without expanding the integration.
$ ELBO(q_j) = E_j[E_{-j}[\text{log}p(z_j,z_{-j},x)]] - E_j[\text{log}q_j(z_j)] + \text{constant} $
Where the first term is derived by writting the first term in (21) using the iterated expectation. The second term is derived by decomposing the second term in (21) as $q_j$ are independent with each other.
Actually, it's the same as the equation (22) in the question, just in another form.
Then, it can be found that $E_j[E_{-j}[\text{log}p(z_j,z_{-j},x)]] - E_j[\text{log}q_j(z_j)]$ is the negative KL divengence between $q_j$ and $\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}$.
$ELBO(q_j)=-\text{KL}(q_j(z_j)||\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}) + \text{constant}$
So, in order to maximize the ELBO, we have to minimize the KL divengence between $q_j$ and $\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}$. The minimal happens when $q_j$ is equal to the $\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}$.
So, we can update the $q_j$ with:
$q^*_j\propto{\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\}}$
The full distribution of $q^*_j$ can be get by normalizing according to the proportional relationship.
$q^*_j = \frac{ \text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\} }{ \int_{z_j}\text{exp}\{E_{-j}[\text{log}p(z_j,z_{-j},x)]\} }$
But in reality, in order to update the variational parameter, the full distribution is usually not necessary to calculate if the conjugate exponential family distribution is used.
The core idea is explained above. For more detail and extension about how the exponential family is used in VI, refer the following link:
