Keep in mind that the left-side $$\int_{-\infty}^{\infty} f(x) \delta(x - a) dx$$ is not an inner product of two functions, since the Dirac delta function isn't a function. It's properly understood in the context of distribution theory, or regarded as a generalized function. That being said, you're right: The utility isn't in computing the value at $a$, but rather in knowing that the integral can be simplified.
Another nice reason to study the delta function is that it can be regarded (again, in the sense of distributions) as the derivative of the Heaviside function.
As an example of this that comes up in differential equations, suppose you have a mass hanging on a spring; the displacement $y$ from the equilibrium can then be modeled by
$$y'' + y = 0$$
which has oscillatory solutions. More generally, if we have a driving force $f(t)$ applied to the system, we have
$$y'' + y = f(t)$$
Now suppose we want to deliver a sharp impulse to the system to get it going initially; the delta function can be thought of as the limit of very short, but strong impulses acting on the system (and hence, as an instantaneous shock). This can be studied via the "ODE"
$$y'' + y = \delta(t - 1)$$
subject to initial conditions $y(0) = y'(0) = 0$ (so our block is not moving, and not displaced at time $0$; then we hit it $1$ second later). A standard way to study such equations is via the Laplace transform; recall that
$$\mathcal{L}\Big(f(t)\Big)(s) = \int_0^{\infty} e^{-st} f(t) dt$$
Now notice that
$$\mathcal{L}\Big(\delta(t-1)\Big)(s) = \int_0^{\infty} e^{-st} \delta(t - 1) dt = e^{-s}$$
The left-hand side of the equation can be transformed easily, and then standard techniques (e.g. a table of Laplace transforms) can give the final solution quite easily from here.
The key point is that the niceness with which $\delta$ interacts with integrals allows us to extend some techniques to difficult situations.
$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
$$
\mbox{In spherical coordinates,}\quad
\delta\pars{\vec{r}}={\delta\pars{r}\delta\pars{\cos\pars{\theta}}\delta\pars{\phi} \over r^{2}}\quad
\mbox{such that}
$$
\begin{align}
\color{#66f}{\large\int_{{\mathbb R}^{3}}\delta\pars{\vec{r}}\,\dd^{3}\vec{r}}
&=\int_{0^{-}}^{\infty}\dd r\,r^{2}\int_{0}^{\pi}\dd\theta\,\sin\pars{\theta}
\int_{0}^{2\pi}\dd\phi\,{\delta\pars{r}\delta\pars{\cos\pars{\theta}}\delta\pars{\phi} \over r^{2}}
\\[3mm]&=\underbrace{\bracks{\int_{0^{-}}^{\infty}\delta\pars{r}\,\dd r}}
_{\ds{=\ 1}}\
\underbrace{\bracks{%
\int_{0}^{\pi}\delta\pars{\cos\pars{\theta}}\sin\pars{\theta}\,\dd\theta}}
_{\ds{=\ 1}}\
\underbrace{\bracks{\int_{0}^{2\pi}\delta\pars{\phi}\,\dd\phi}}_{\ds{=\ 1}}\
\\[3mm]&=\ \color{#66f}{\Large 1}
\end{align}
$$\mbox{Note that}\quad
\int_{{\mathbb R}^{3}}\delta\pars{\vec{r} - \vec{r}_{0}}\,\dd^{3}\vec{r}
=\int_{{\mathbb R}^{3}}\delta\pars{\vec{r}}\,\dd^{3}\vec{r}
$$
Best Answer
Coming from a measure-theoretic point of view, one can write the following:
$\int_{-\infty}^{+\infty}\delta(x-a)dx = \int_{\mathbb R} 1 d \mu_a(x)$, where $\mu_a$ denotes the point measure with mass $1$ at $a$ and $0$ everywhere else.
Now, if we want to integrate over an interval $I$, we integrate the characteristic function $\chi_I(x)$ over $\mathbb R$:
$$\int_I \delta(x-a)dx= \int_{\mathbb R}\chi_I(x)d\mu_a(x).$$
And here comes a important distinction which we dont need to make when using the Riemann integral: If $a$ is a border point of the interval $I$, do we consider the open interval $I_1 =(-\infty,a)$, or the interval $I_2=(-\infty,a]$ ? Because the mass of $\mu_a$ lies at $a$, this is an important distinction:
$$\int_{(-\infty,a)} \delta(x-a)dx= \int_{\mathbb R}\chi_{(-\infty,a)}(x)d\mu_a(x) = 0 $$
while
$$\int_{(-\infty,a]} \delta(x-a)dx= \int_{\mathbb R}\chi_{(-\infty,a]}(x)d\mu_a(x) = 1.$$
So we have to be very careful with what we mean by $\int_{-\infty}^a\delta(x-a)dx.$
Summing up: The fact that it's 'infinitely thin' translates to the measure having a point mass, and this means that the inclusion or the exclusion of single points will change the value of the integral.