Solution of $\int p_\theta(z) \log q(z) dz$ of Gaussian case

Following is from the original paper of concept of VAE(variational autoencoder) by Kingma,Welling 2014

B. Solution of $$D_{KL}(p_\phi(z)||q_\theta(z))$$ of Gaussian case

The variational lower bound (the objective to be maximized) contains a KL term that can often be integrated analytically. Here we give the solution when both the prior $$q_{\theta}(z) = N (0, I)$$ and the posterior approximation $$p_\phi(z|x^{ (i)})$$ are Gaussian. Let $$J$$ be the dimensionality of $$z$$. Let $$\mu$$ and $$\sigma$$ denote the variational mean and standard deviation evaluated at datapoint $$i$$, and let $$\mu_j$$ and $$\sigma_j$$ simply denote the $$j$$-th element of these vectors. Then: $$\int p_\theta(z) \log q(z) dz\\ = \int N (z; \mu,\sigma^2 ) \log N (z; 0, I) dz\\ = − {J\over 2} \log(2\pi) − {1\over 2} \sum_{j=1}^{J} (\mu_j{^2} + \sigma_j^2 )$$

At the equation above can't understand how the second equality calculated. Any hint to understand those eqaulity?

In this case, I find it easier to start from a simple case and then build up in complexity. The simplest case is to consider $$J=1$$.
\begin{align} -\int_{-\infty}^{\infty} p(z) \log(q(z)) dz &= \frac{1}{2}\log(2\pi\sigma_2^2) - \int p(z) \left(-\frac{\left(z - \mu_2\right)^2}{2 \sigma_2^2}\right)dz \\ &= \frac{1}{2}\log(2\pi\sigma_2^2) + \frac{\mathbb{E}_{z\sim p}[z^2] - 2 \mathbb{E}_{z\sim p}[z]\mu_2 +\mu_2^2} {2\sigma_2^2} \\ &= \frac{1}{2}\log(2\pi\sigma_2^2) + \frac{\sigma_1^2 + \mu_1^2-2\mu_1\mu_2+\mu_2^2}{2\sigma_2^2} \\ &= \frac{1}{2}\log(2\pi\sigma_2^2) + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2\sigma_2^2}\\ \int_{-\infty}^{\infty} p(z) \log(q(z)) dz = V &= -\frac{1}{2}\log(2\pi\sigma_2^2)-\frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2\sigma_2^2} \end{align} The key is recognizing that we can expand the quadratic in the first line; this gives us a sum of several integrals. Then we apply the law of the unconscious statistician, and we use the fact that $$\text{Var}(z)=\mathbb{E}[z^2]-\mathbb{E}[z]^2$$. The rest is just rearranging.
In this special case, we know that $$q$$ is a standard normal, so \begin{align} -\frac{1}{2}\log(2\pi\sigma_2^2)-\frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2\sigma_2^2} &= -\frac{1}{2} \log(2\pi) -\frac{1}{2}\left(\sigma_1^2 + \mu_1^2\right) \end{align}
How can we generalize this to $$J>1$$? Consider the case of a diagonal covariance matrix. In this case, the $$z_j$$ are independent. So the solution arises from the sum of $$V_j$$. If you're not convinced, then you'll need to crank through the matrix arithmetic, using multivariate normal $$p$$ and $$q$$. It's not particularly hard, it's just tedious. Here's some threads to get started: