The left hand side is:
$$\sum_{\mathbf r \in \Lambda} e^{-i\mathbf k \cdot \mathbf r}=\sum_{n_1=-\infty}^{\infty}\sum_{n_2=-\infty}^{\infty}\sum_{n_3=-\infty}^{\infty}e^{-ik_1a_1n_1}e^{-ik_2a_2n_2}e^{-ik_3a_3n_3}$$
$$=\prod_{i=1}^3 \sum_{n_i=-\infty}^{\infty}e^{-ik_ia_in_i}=\prod_{i=1}^3 \sum_{n_i=-\infty}^{\infty}e^{-i2\pi m_in_i}$$
where I've used the fact that $k_i \equiv \frac{2\pi}{a_i}m_i$ for some $m_i\in \mathbb Z$. Now the inside sum is of the form $f(x)\equiv\sum_{n=-\infty}^\infty e^{-i2\pi nx}$. I claim that in fact:
$$f(x)=\sum_{l=-\infty}^\infty \delta(x-l) \qquad (*)$$
How do we show this? Well, first of all, $f(x)$ is periodic with period $1$, because:
$$f(x+1)=\sum_{n=-\infty}^\infty e^{-i2\pi n(x+1)}=\sum_{n=-\infty}^\infty 1\times e^{-i2\pi nx}=f(x)$$
Now, let's restrict ourselves to one period, for example $x\in(-\frac12,\frac12]$. The integral of $f$ in this region is:
$$\int_{-1/2}^{1/2}dx f(x) = \sum_{n=-\infty}^\infty\int_{-1/2}^{1/2}dx e^{-i2\pi nx}=\sum_{n=-\infty}^\infty \frac{\sin(\pi n)}{\pi n}=1$$
The last sum is one since $\sin(n\pi)\equiv 0$, so the only term contributing is the $n=0$ term which looks like a $\frac{\sin x}{x} \to 1$ kind of limit. Also note that in this period, $f(x)=0$ for any nonzero $x$. This is because corresponding to each term $e^{-i2\pi nx}$ in the sum, there is also a term which cancels it out (evaluate the geometric sum explicitly if this doesn't convince you!).
Putting all this together, our function $f$ satisfies all of the defining properties of the Dirac distribution in each period (the integral is one, and it's zero everywhere other than the origin). All we now have to do is to put infinitely many of these Dirac deltas together (with a spacing of the period 1), to ensure periodicity; meaning that we have derived equation $(*)$.
Using this result, the desired sum is:
$$\sum_{\mathbf r \in \Lambda} e^{-i\mathbf k \cdot \mathbf r}=\prod_{i=1}^3 \sum_{l_i=-\infty}^{\infty}\delta(m_i-l_i)$$
Finally, we want to express the result in terms of the momentum $k_i = \frac{2\pi}{a_i}m_i$. Using the fact that $\delta(ax)=\delta(x)/|a|$, we have:
$$\sum_{\mathbf r \in \Lambda} e^{-i\mathbf k \cdot \mathbf r}=\prod_{i=1}^3 \sum_{l_i=-\infty}^{\infty}\delta(\frac{a_i}{2\pi}k_i-l_i)$$
$$=\prod_{i=1}^3 \sum_{l_i=-\infty}^{\infty}(\frac{2\pi}{a_i})\delta(k_i-\frac{2\pi}{a_i}l_i)$$
$$=\frac{(2\pi)^3}{a_1a_2a_3} \sum_{l_1=-\infty}^{\infty}\sum_{l_2=-\infty}^{\infty}\sum_{l_3=-\infty}^{\infty}\delta(k_1-\frac{2\pi}{a_1}l_1)\delta(k_2-\frac{2\pi}{a_2}l_2)\delta(k_3-\frac{2\pi}{a_3}l_3)$$
Finally, realizing that $\frac{(2\pi)^3}{a_1a_2a_3} \equiv V^*$, and that the $\frac{2\pi}{a_i}l_i$ terms are simply components of momenta situated on the reciprocal lattice, we get the final result:
$$\sum_{\mathbf r \in \Lambda} e^{-i\mathbf k \cdot \mathbf r}=V^* \sum_{\mathbf q \in \Lambda^*}\delta(\mathbf k - \mathbf q)$$
Note that the $V^*$ is essential in making the units consistent!
