An intuitive explanation is that we want $q$ to be large whenever either $p$ or $|f|$ is large. Otherwise, our estimate of $E_p[f]$ might have a lot of error, since we're "missing out" on sampling the most influential regions of the real number line. A proof is below:
So we want to minimize
$$\begin{align}
\text{Var}_q\left[ \frac{p(x)f(x)}{q(x)} \right]
&= E_q\left[ \left( \frac{p(x)f(x)}{q(x)} \right)^2 \right] - E_q\left[\frac{p(x)f(x)}{q(x)} \right]^2 \\
\end{align}$$
The second term is constant with respect to $q$. In fact, it comes out to exactly $E_p[f(x)]^2$, so we can drop it from the optimization, which leaves with
$$\begin{align}
E_q\left[ \left( \frac{p(x)f(x)}{q(x)} \right)^2 \right] &= \int \frac{p(x)^2 f(x)^2}{q(x)} dx
\end{align}$$
We also have the constraint that $\int q(x) dx = 1$. Since this is a constrained optimization problem, we can write the lagrangian:
$$\begin{align}
L(q, \lambda) &= \int \frac{p(x)^2 f(x)^2}{q(x)} dx + \lambda \left( \int q(x) dx - 1 \right)
\end{align}$$
We want the functional derivative with respect to $q$:
$$\begin{align}
\frac{\partial L(q, \lambda)}{\partial q(x)} &= \lim_{\epsilon \rightarrow 0} \frac{\partial}{\partial \epsilon}
\left[ \int \frac{p(x)^2 f(x)^2}{q(x) + \epsilon \eta(x)} dx + \lambda \left( \int (q(x)+\epsilon \eta(x)) dx - 1 \right) \right]
\end{align}$$ for any arbitrary $\eta$.
This comes out to
$$\begin{align}
&\lim_{\epsilon \rightarrow 0}
\left[ \int -\frac{p(x)^2 f(x)^2}{(q(x) + \epsilon \eta(x))^2} \eta(x)\ dx + \lambda \int \eta(x) dx \right] \\
&= \int -\frac{p(x)^2 f(x)^2}{q(x)^2} \eta(x)\ dx + \lambda \int \eta(x) dx \\
&= \int \eta(x)\left( \lambda -\frac{p(x)^2 f(x)^2}{q(x)^2}\right) dx \\
\end{align}$$
Since we want the derivative to be $0$ for all $\eta(x)$, then we must have
$$\begin{align}
0 &=\lambda -\frac{p(x)^2 f(x)^2}{q(x)^2} \\
\lambda &= \frac{p(x)^2 f(x)^2}{q(x)^2} \\
q(x)^2 &= \frac{p(x)^2 f(x)^2}{\lambda} \\
q(x) &= \frac{p(x)|f(x)|}{\sqrt{\lambda}}
\end{align}$$
And it must be the case that $\sqrt{\lambda} = Z$.
