Where does the relation $\nabla^2(1/r)=-4\pi\delta^3({\bf r})$ between Laplacian and Dirac delta function come from? It is often quoted in physics textbooks for finding the electric potential using Green's function that 
$$\nabla ^2 \left(\frac{1}{r}\right)=-4\pi\delta^3({\bf r}),$$ 
or more generally 
$$\nabla ^2 \left(\frac{1}{|| \vec x - \vec x'||}\right)=-4\pi\delta^3(\vec x - \vec x'),$$
where $\delta^3$ is the 3-dimensional Dirac delta distribution.
However I don't understand how/where this comes from. Would anyone mind explaining?

We can use the simplest method to display the results, as shown below : -
$$ \nabla ^2 \left(\frac{1}{r}\right) = \nabla \cdot \nabla \left( \frac 1 r \right) = \nabla \cdot \frac {-1 \mathbf {e_r}} {r^2} $$ 
Suppose there is a sphere centered on the origin, then the total flux on the surface of the sphere is : -
$$ \text {Total flux} = 4 \pi r^2 \frac {-1} {r^2} = -4 \pi $$
Suppose the volume of the sphere be $ \mathbf {v(r)}$, so by the definition, the divergence is : -
$$ \lim_{\text {volume} \to zero} \frac {\text {Total Flux}} {\text {Volume}} = \lim_{\text {v(r)} \to 0} \left(\frac {-4 \pi} {v(r)}\right)  $$
So obviously,
$$ \lim_{\text {r} \to 0} \left[ \nabla ^2 \left(\frac{1}{r}\right) \right]= \lim_{\text {r} \to 0} \left[ \nabla \cdot \nabla \left( \frac 1 r \right) \right]= \lim_{\text {r} \to 0} \left(\frac {-4 \pi} {v(r)} \right) = \text {infinite} $$
$$ \lim_{\text r\to 0} \int \nabla ^2 \left( \frac 1 r \right) dv(r) = \lim_{\text r\to 0}\int \frac {-4 \pi} {v(r)} dv(r) = -4\pi$$
Since the laplacian is zero everywhere except $r \to 0$, it is true and real that :-
$$ \nabla ^2 \left(\frac{1}{r}\right)=-4\pi\delta^3({\bf r}) $$
