Solution (source code)

= Solution

<Normal-normal conjugacy> gives
$$
p(\theta\mid y)=N(\widetilde\theta,\sigma_\theta^2),
\qquad
\sigma_\theta^2=\frac{\sigma^2\tau^2}{\sigma^2+\tau^2},
\qquad
\widetilde\theta=\frac{\tau^2}{\sigma^2+\tau^2}y.
$$
The three quantities in the proposed identity are
$$
\log Z=-\frac12\left[
\log(2\pi(\sigma^2+\tau^2))+
\frac{y^2}{\sigma^2+\tau^2}
\right],
$$
$$
\mathbb E_{\theta\mid y}\log L(\theta)
=-\frac12\log(2\pi\sigma^2)
-\frac{(y-\widetilde\theta)^2+\sigma_\theta^2}{2\sigma^2},
$$
and, by the <Kullback-Leibler divergence between normal distributions>,
$$
D_{\mathrm{KL}}(p(\theta\mid y)\Vert\pi)
=\frac12\left[
\log\frac{\tau^2}{\sigma_\theta^2}
+\frac{\sigma_\theta^2+\widetilde\theta^2}{\tau^2}-1
\right].
$$
Substitution and simplification show that the latter two expressions differ by exactly $\log Z$, so the equality holds.