Is electrostatic energy positive definite?

One of the integrations can be regarded as a convolution that yields the potential:

$$ E[\rho]=\frac1{8\pi}\int\mathrm d^3x\rho(\vec x)V(\vec x) $$

with

$$ V(\vec x)=\int\mathrm d^3y\frac{\rho(\vec y)}{|\vec x-\vec y|}\;. $$

This is $V=\rho*\dfrac1r$, where the asterisk denotes convolution. A convolution in real space corresponds to multiplication in Fourier space, so $\mathcal F(V)=\mathcal F(\rho)\mathcal F(1/r)$, where $\mathcal F(\cdot)$ denotes the Fourier transform. Since the Fourier transform is unitary, we have, with the Fourier transform $\mathcal F(1/r)=4\pi/k^2$ of the Coulomb potential,

$$ \begin{align} E[\rho] &= \frac1{8\pi}\int\mathrm d^3x\rho(\vec x)V(\vec x) \\ &= \frac1{8\pi}\int\mathrm d^3k\mathcal F(\rho)(\vec k)\mathcal F(V)(\vec k) \\ &= \frac1{8\pi}\int\mathrm d^3k\mathcal F(\rho)(\vec k)\mathcal F(\rho)(\vec k)\mathcal F (1/r)(\vec k) \\ &= \frac12\int\mathrm d^3k\left(\mathcal F(\rho)(\vec k)\right)^2\frac1{k^2}\;, \end{align} $$

which is manifestly positive definite.

Joriki's method in the answer above allows one to say something about other sorts of energy $E[\rho]=\int dxdy \left(\rho(x)\rho(y)f(|x-y|)\right)$ as well. In the case where $f(r) = 1/r$ (and in three dimensions) there is however a more physical proof which can be summarized as \begin{equation} \int \frac{\rho(x)\rho(y)}{4\pi \epsilon_0 |x-y|}=\int \rho V = \epsilon_0\int \left(\nabla \cdot \mathbf{E}\right)V \underbrace{=}_{P.I.} -\epsilon_0\int \mathbf{E}\cdot\nabla V = \epsilon_0\int \mathbf{E}^2 \geq 0. \end{equation} There are of course some minor technical assumptions to be made for these manipulations to make sense.