To add some clarity to the original paper you cited in the previous post and make notations consistent, I'll write the true density as $p(x)$ and the approximating density as $q_\theta(x)$ parameterized by $\theta$. Since you haven't specified the true density (the one being approximated), I'll try to (re-)explain what the proof means in each step.
The goal is to show that as long as the approximating density $q_\theta(x)$ belongs to an exponential family, minimizing the Kullback-Leibler (KL) divergence $\mathrm{KL}(p\| q_\theta)$ only requires matching the sufficient statistics. First, look at the definition of the KL divergence:
\begin{align}
\mathrm{KL}(p\| q_\theta) &= \int\log \frac{p(x)}{q_\theta(x)}\, p(x)\, dx \\
&= \mathrm{E}_{p(x)}\left(\log \frac{p(x)}{q_\theta(x)} \right) \\
&= \mathrm{E}_{p(x)}(\log p(x)) - \mathrm{E}_{p(x)}(\log q_\theta(x)).
\end{align}
Since we need to minimize this (as a function of the parameters $\theta$), we will make this point clear by rewriting it as $f(\theta) = \mathrm{KL}(p\| p_\theta)$, and use the first-order condition $\nabla_\theta f(\theta) = 0$. We see that $\mathrm{E}_{p(x)}(\log p(x))$ in the KL divergence disappears upon differentiation because it's not a function of $\theta$. Therefore,
$$
\nabla_\theta f(\theta) = -\nabla_\theta \mathrm{E}_{p(x)}(\log q_\theta(x)).
$$
Now, let's go back to the definition of exponential family: $q_\theta(x) = h(x) \exp\{\theta^\top \phi(x) - A(\theta) \} $ where $A(\theta)$ is the log-normalizing constant. Taking the logarithm of this density yields
$$
\begin{align}
\log q_\theta(x) &= \log h(x) + \theta^\top \phi(x) - A(\theta)\\
\mathrm{E}_{p(x)}(\log q_\theta(x)) &= \mathrm{E}_{p(x)}(\log h(x)) + \theta^\top \mathrm{E}_{p(x)}(\phi(x)) - A(\theta)\\
\nabla_\theta \mathrm{E}_{p(x)}(\log q_\theta(x)) &= \mathrm{E}_{p(x)}(\phi(x)) - \nabla_\theta A(\theta).
\end{align}
$$
But as I've derived in this answer, $\nabla_\theta A(\theta) = \mathrm{E}_{q_\theta(x)}(\phi(x))$. Therefore, the first-order condition $\nabla_\theta f(\theta) = 0$ gives us $\mathrm{E}_{p(x)}(\phi(x)) = \mathrm{E}_{q_\theta(x)}(\phi(x))$.
To verify that solving the first-order condition for $\theta$ gives us the minimizer, we must compute the Hessian matrix and check if it's positive-definite. From $\nabla_\theta \mathrm{E}_{p(x)}(\log q_\theta(x)) = \mathrm{E}_{p(x)}(\phi(x)) - \nabla_\theta A(\theta)$, it's easily observed that $\nabla_\theta^2 f(\theta) = \nabla_\theta^2 A(\theta)$, which is the covariance matrix of the sufficient statistics. I will not prove this because this is a standard result in mathematical statistics, but refer to lecture notes like this if you're interested. Again, $[\nabla_\theta^2 A(\theta)]_{ij} = \mathrm{Cov}(\phi_i(x),\phi_j(x))$. Covariance matrices are by definition positive-definite, and therefore the first-order condition when solved for $\theta$ indeed produces the minimizer.
Now, since you've asked how this plays out for normal-gamma, we've already established that the sufficient statistics are $\phi_1(x) = \log T$, $\phi_2(x) = T$, $\phi_3(x)=TX$, and $\phi_4(x)=TX^2$. To obtain the covariance matrix in full, you should compute 10 second-order derivatives $\frac{d^2}{d\theta_i d\theta_j} A(\theta)$ for $i,j=1,\ldots,4$ for
$$
A(\theta)= -\left(\theta_1+\dfrac{1}{2} \right)\log\left(\frac{\theta_3^2}{4\theta_4} - \theta_2 \right) -\frac{1}{2}\log(-2\theta_4) +\log\Gamma\left(\theta_1+\frac{1}{2}\right) + \frac{1}{2}\log(2\pi),
$$
for which I seriously recommend that you consider using tools like Wolfram Alpha.
Best Answer
First, note that normal-gamma distribution is an exponential family and has the following form: $$ f(x;\eta) = h(x)\exp\{\color{blue}{\eta^\top t(x)} \color{red}{-A(\eta)} \}, $$ where $\eta$ is the natural parameter, $t(x)$ is the sufficient statistic, $h(x)$ is the underlying measure, and $A(\eta)$ is the log-normalizing constant. What's interesting about the log-normalizing constant $A(\eta)$ is that its derivatives provide the moments of the sufficient statistic, i.e., $$ \begin{align} \frac{d}{d\eta}A(\eta) &= \dfrac{d}{d\eta}\left(\log \int \exp\{\eta^\top t(x)\} h(x)\,dx \right)\\ &= \frac{\int t(x)\exp\{\eta^\top t(x)\}h(x)\,dx}{\int \exp\{\eta^\top t(x)\}h(x)\,dx}\\ &= \int t(x) \exp\{\eta^\top t(x)-A(\eta)\}h(x)\,dx\\ &= E(t(X)). \end{align} $$
Let's rewrite the density in terms of its natural parameter and sufficient statistic. Observe that $$ \begin{align} f(X,T;\mu,\lambda,\alpha,\beta) &= \exp\Bigg\{\color{blue}{\left(\alpha-\frac{1}{2} \right)\log T -\left(\beta + \frac{\lambda \mu^2}{2} \right)T + \lambda \mu XT - \frac{\lambda}{2} TX^2} \\ &\quad \quad\color{red}{+\alpha \log \beta + \frac{1}{2}\log \lambda - \log \Gamma(\alpha) - \frac{1}{2}\log (2\pi)}\Bigg\}. \end{align} $$ This part is merely matching the shape of the general form of an exponential family and a specific example. I colored the parts so each component is clearly separated from one another. Now, it's pretty clear from the blue part that $$ \eta_1 = \alpha - \frac{1}{2},\;\eta_2 = -\beta - \frac{\lambda\mu^2}{2},\;\eta_3=\lambda\mu,\;\eta_4=-\frac{\lambda}{2}, $$ and $$ T_1 = \ln T,\;T_2=T,\; T_3=TX,\;T_4=TX^2, $$ so that $\eta^\top t(x) = \eta_1T_1 + \eta_2T_2 + \eta_3T_3 + \eta_4T_4$.
Next step is to reexpress $A(\eta)$ in terms of $\eta$. Note that $\alpha= \eta_1+\frac{1}{2}$, $\beta=-\eta_2 + \frac{\eta_3^2}{4\eta_4}$, $\mu = -\frac{\eta_3}{2\eta_4}$, and $\lambda = -2\eta_4$. Then, $$ \begin{align} A(\eta)&= -\left(\eta_1+\dfrac{1}{2} \right)\log\left(\frac{\eta_3^2}{4\eta_4} - \eta_2 \right) -\frac{1}{2}\log(-2\eta_4) +\log\Gamma\left(\eta_1+\frac{1}{2}\right) + \frac{1}{2}\log(2\pi). \end{align} $$
From now on, it's good ol' differentiation: $$ \begin{align} \frac{d}{d\eta_1}A(\eta) &= -\log\left(\frac{\eta_3^2}{4\eta_4} - \eta_2 \right)+ \psi\left( \eta_1+\frac{1}{2}\right) = \psi(\alpha)-\log(\beta),\\ \frac{d}{d\eta_2}A(\eta) &= \frac{\eta_1+\frac{1}{2}}{\frac{\eta_2^2}{4\eta_4}-\eta_2} = \frac{\alpha}{\beta},\\ \frac{d}{d\eta_3}A(\eta) &= -\left(\eta_1+\frac{1}{2} \right)\frac{\eta_3/(2\eta_4)}{\frac{\eta_3^2}{4\eta_4} - \eta_2} = \frac{\alpha \mu}{\beta},\\ \frac{d}{d\eta_4}A(\eta) &= -\left(\eta_1+\frac{1}{2} \right)\frac{-\frac{\eta_3^2}{4\eta_4^2}}{\frac{\eta_3^2}{4\eta_4} - \eta_2} - \frac{1}{2\eta_4} = \frac{\alpha\mu^2}{\beta} + \frac{1}{\lambda}. \end{align} $$