This is a question coming from physics, but its nature is purely mathematical. Given some continuous distribution of charge $\rho$ (take it compactly supported, or "nice enough" depending on the problem you are treating), we define the electrostatic energy as:
$$E[\rho]=\frac{1}{8\pi}\int d^3x\int d^3y \frac{\rho(\vec{x})\rho(\vec{y})}{|\vec{x}-\vec{y}|}$$
Is it positive definite? How to show it?
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.
Here is another nice angle:
In our take, we can allow the charge distribution to be a signed (Borel-)measure $\mu=\mu_+-\mu_-$ ($\mu_+$ and $\mu_-$ ordinary Borel measures). The integrals $$I_+:=\int_{\mathbb{R}^3 \times \mathbb{R}^3}d\mu_+(x)d\mu_+(y)\,|x-y|^{-1},\qquad I_-:=\int_{\mathbb{R}^3 \times \mathbb{R}^3}d\mu_-(x)d\mu_-(y)\,|x-y|^{-1}$$ are well-defined. Let us first assume that $I_+<+\infty>I_-$. We will borrow a bit from the identity $$|x-y|^{-1}=\pi^{-3} \int_{\mathbb{R}^3}d\lambda(z)\,|x-z|^{-2}|z-y|^{-2}.$$ For our purpose we do not need to know that the exact proportionality constant in that identity is $\pi^{-3}$ and then the proof of the identity becomes exceedingly simple after observing that the RHS is homogeneous of degree $-1$ and rotational invariant. The following is then a consequence of Fubini's Theorem (which, recall, requires only the negative or positive part of an integrand to be integrable for its validity): $$\int_{\mathbb{R}^3 \times \mathbb{R}^3}d\mu(x)d\mu(y)\,|x-y|^{-1}\\ =\pi^{-3}\int_{\mathbb{R}^3}d\lambda(z)\left(\left(\int_{\mathbb{R}^3}d\mu(x)\,|x-z|^{-2}\right)\left(\int_{\mathbb{R}^3}d\mu(y)\,|y-z|^{-2}\right)\right)\\ =\pi^{-3}\int_{\mathbb{R}^3}d\lambda(z)\left(\int_{\mathbb{R}^3}d\mu(x)\,|x-z|^{-2}\right)^2\geq 0.$$
