This is a very interesting story which you usually don't get in an undergraduate course on electromagnetism.
For simplicity, consider just two charges $q_1$ at $\mathbf{r}_1$ and $q_2$ at $\mathbf{r}_2$. Their fields are
$$\mathbf{E}_1=q_1\frac{\mathbf{r}-\mathbf{r}_1}{|\mathbf{r}-\mathbf{r}_1|^3}$$
and
$$\mathbf{E}_2=q_2\frac{\mathbf{r}-\mathbf{r}_2}{|\mathbf{r}-\mathbf{r}_2|^3}$$
where I'll use Gaussian units instead of SI units since factors of $\frac{1}{4\pi\epsilon_0}$ are just noise.
The integral for the electrostatic potential energy stored in the combined field $\mathbf{E}_1+\mathbf{E}_2$ is
$$U=\frac{1}{8\pi}\int(\mathbf{E}_1+\mathbf{E}_2)^2\,d^3\mathbf{r}.$$
Expanding the square of the field, there are three terms:
$$U_1=\frac{q_1^2}{8\pi}\int\frac{d^3\mathbf{r}}{|\mathbf{r}-\mathbf{r}_1|^4},$$
$$U_2=\frac{q_2^2}{8\pi}\int\frac{d^3\mathbf{r}}{|\mathbf{r}-\mathbf{r}_2|^4},$$
and
$$U_{12}=\frac{q_1q_2}{4\pi}\int\frac{\mathbf{r}-\mathbf{r}_1}{|\mathbf{r}-\mathbf{r}_1|^3}\cdot\frac{\mathbf{r}-\mathbf{r}_2}{|\mathbf{r}-\mathbf{r}_2|^3}\,d^3\mathbf{r}.$$
The first two are divergent "self-energy" integrals representing the infinite electrostatic energy of each point charge by itself. They are infinite positive constants that don't depend on the separation of the two charges. They can be ignored; think of them as being a contribution to the mass-energy of the charges, infinitely "renormalizing" their masses, in a classical version of the renormalization you encounter in quantum field theory.
The third, the "interaction-of-$\mathbf{E}_1$-with-$\mathbf{E}_2$" integral, is a finite integral. With some effort it can be performed analytically. (I will not give the details; it's a nontrivial calculation and this answer is already long.) The result turns out to be
$$U_{12}=\frac{q_1q_2}{|\mathbf{r}_1-\mathbf{r}_2|},$$
which is simply the usual formula for the electrostatic potential energy of two point charges. This integral can be positive or negative depending on the signs of $q_1$ and $q_2$. It isn’t positive-definite because the integrand isn't the square of anything; it’s $\mathbf{E}_1\cdot\mathbf{E}_2$ and dot products can be positive or negative. In fact, in some regions the integrand is positive and in others it’s negative, as you can see by drawing a diagram with some field arrows for the two fields.
So now you know that the usual formula for the electrostatic potential energy of two point charges is the interaction energy of their two individual fields, not the total field energy (which is infinite).