Calculus method
The Taylor series of a function of $d$ variables is as follows:
$$f(\mathbf{x} + \mathbf{y}) = \sum_{n_1=0}^{\infty}\ldots\sum_{n_d=0}^{\infty}
\frac{(y_1-x_1)^{n_1} \ldots (y_d - x_d)^{n_d}}{n_1!\ldots n_d!}
\left. \left(
\partial_1^{n_1}\ldots \partial_d^{n_d} \right) f \right|_{\mathbf{x}}.$$
where $\partial_i^{n_i} f$ means "the $n_i^{\text{th}}$ order partial derivative of $f$ with respect to the $i^{\text{th}}$ coordinate".
Consider only the terms where $\sum_{i=1}^{\infty} n_i = 1$.
These are the terms for which there is only one derivative being taken, a.k.a. the "linear terms".
Keeping the $0^{\text{th}}$ order term and the linear terms gives
$$f(\mathbf{x} + \mathbf{y}) = f(\mathbf{x})+\sum_{i=1}^{d}(y_i - x_i) \left. (\partial_i f \right)|_{\mathbf{x}} . \qquad (*)$$
Define the displacement $\epsilon = \mathbf{y} - \mathbf{x}$. Then
$$\epsilon \cdot \nabla = \sum_{i=1}^d (y_i - x_i) \partial_i$$
by definition of what $\cdot$ means. Therefore, keeping only up to linear terms, $(*)$ becomes
$$f(\mathbf{y} - \mathbf{x}) = f(\mathbf{x}) + \left. (\epsilon \cdot \nabla)\right|_{\mathbf{x}}f . $$
So you see, your formula is correct, just in a particular notation you may not be used to.
Algebraic method
For the sake of typing less, I'm going to explain this in one dimension, but nothing in the actual calculation is restricted to one dimension.
The state vector can be thought of as a linear combination of eigenvectors of the position operator
$$|\Psi\rangle = \int_x \Psi(x) |x\rangle \qquad (1)$$
where $\hat{x}|x\rangle = x|x\rangle$.
Bringing in a bra $\langle y |$ from the left and using $\langle y | x \rangle = \delta(x-y)$ in (1) gives
$$\langle y | \Psi \rangle = \Psi(y) . \qquad (2)$$
Combining (1) and (2) gives
$$|\Psi\rangle = \int_x |x\rangle \langle x| \Psi \rangle \qquad (3)$$
which shows that $\int_x |x\rangle \langle x | = \text{Identity}$.
Let's really understand what all this means:
Eq. (1) just says that you can write a vector as a linear combination of basis vectors.
In this case, $\Psi(x)$ are the coefficients in the linear combination.
In other words, the wave function is just the coefficients of a linear combination expansion written in the position basis.
Eq. (2) makes this explicit by showing that the inner product of a position basis vector $|y\rangle$ with $|\Psi\rangle$ is precisely $\Psi(y)$.
Eq. (3) just expresses a neat way of expressing the identity operator in a way which we will find very useful.
Note that this works with any basis.
For example,
$$\int_p |p\rangle \langle p | = \text{identity} .$$
We can use this to express a position eigenvector in terms of momentum eigenvectors,
$$|x\rangle = \int_p |p\rangle \langle p | x \rangle = \int_p e^{-i x p/\hbar}|p\rangle ,$$
where we've used the fact that $\langle x | p \rangle = e^{i p x/\hbar}$ [1].
Now let's get back to the question at hand. Define $|\Psi'\rangle \equiv e^{i \epsilon \hat{p}/\hbar}|\Psi\rangle$.
Let us evaluate $\Psi'(y) \equiv \langle y|\Psi'\rangle$:
\begin{align}
\Psi'(y) &= \langle y | \Psi' \rangle \\
&= \langle y| e^{i \epsilon \hat{p}/\hbar} |\Psi \rangle \\
&= \langle y| e^{i \epsilon \hat{p}/\hbar} \int_x \Psi(x) |x\rangle \\
&= \int_x \Psi(x) \langle y | e^{i \epsilon \hat{p}/\hbar} |x \rangle .
\end{align}
To proceed we need to compute $e^{i \epsilon \hat{p}/\hbar} |x \rangle$.
We can do this using our expression for the identity in terms of momentum states,
$$
\begin{align}
e^{i \epsilon \hat{p}/\hbar} |x \rangle
&= e^{i \epsilon \hat{p}/\hbar} \left( \int_p |p\rangle \langle p | \right) | x \rangle \\
&= \int_p e^{i \epsilon \hat{p} / \hbar} |p\rangle \langle p | x \rangle \\
&= \int_p e^{i \epsilon p / \hbar} |p\rangle \langle p | x \rangle \\
&= \int_p e^{i \epsilon p /\hbar} |p\rangle e^{-i p x/\hbar}\\
&= \int_p e^{-i(x - \epsilon) p / \hbar} |p\rangle \\
&= |x - \epsilon \rangle
\end{align}$$
This is actually the thing you should remember about this question:
$$e^{i\epsilon\hat{p}/\hbar}|x\rangle = |x - \epsilon \rangle .$$
This is an excellent equation because it tells you something about how an the $e^{i \epsilon \hat{p} /\hbar}$ operator changes a position eigenvector without having to write it out in a particular basis.
Now that we have that, we can go back to what we were trying to compute,
$$
\begin{align}
\Psi'(x) &= \int_x \Psi(x) \langle y | e^{i \epsilon \hat{p}/\hbar} |x \rangle \\
&= \int_x \Psi(x) \langle y | x - \epsilon \rangle \\
&= \int_x \Psi(x) \delta(y - (x - \epsilon)) \\
&= \Psi(y + \epsilon) .
\end{align}
$$
This is the equation you wanted to prove.
If I messed up any minus signs I hope someone will edit :)
[1] I'm not proving that, and if you want to know why that's true it would make a good question on this site.
Best Answer
No, it's not an approximation. There are a number of ways in which one might define the momentum operator, and one of those ways is as follows. In the following I set $\hbar = 1$ for convenience.
If you have a family of unitary operators $U_\lambda$ indexed by a continuous variable $\lambda\in \mathbb R$ which has the properties that $U_0 = \mathbb I$ (the identity operator) and that $U_\lambda \circ U_\rho = U_{\lambda + \rho}$, then it turns out that you can write $U_\lambda = \exp[-i \lambda A]$ where $A$ is a self-adjoint operator which is called the infinitesimal generator of the family of operators $U_\lambda$. In order to extract $A$ from $U_\lambda$, we use a difference quotient$^1$: $$A\psi(x) := i\lim_{\lambda\rightarrow 0} \frac{\big(U_\lambda \psi\big)(x)- \psi(x)}{\lambda}$$
In this particular case, we let $\big(U_\lambda\psi\big)(x) = \psi(x-\lambda)$, which means that $U_\lambda$ are the spatial translation operators$^2$ which shift the wavefunction by $\lambda$; we then find that $$A\psi(x) := i\lim_{\lambda\rightarrow 0} \frac{\psi(x-\lambda) - \psi(x)}{\lambda} = -i \psi'(x) \tag{$\star$}$$ and therefore that $A = -i \frac{d}{dx}$.
This approach is formally equivalent to writing $U_\epsilon \approx 1 -i \epsilon A$ and expanding $U_\epsilon \psi$ to first order in $\epsilon$. This yields $$ \big(U_\epsilon \psi\big)(x) = \psi(x-\epsilon)$$ $$\implies \big([1 - i\epsilon A]\psi\big)(x) + \mathcal O(\epsilon^2)=\psi(x) - \epsilon \psi'(x) + \mathcal O(\epsilon^2) $$ $$\implies \big(A\psi\big)(x) = -i \psi'(x)\tag{$\star\star$}$$ just as before.
The main point here is the relationship between families of unitary operators $U_\lambda$ and their corresponding self-adjoint infinitesimal generators $A$. Given a family $U_\lambda$, $A$ is technically defined via a difference quotient as in $(\star)$, but can also be obtained by expanding $U_\lambda$ to first order, as in $(\star\star)$. In the latter case, one might say that $A$ is defined to be the operator in the first-order term when one expands $U_\epsilon$ about $\epsilon=0$; though the first order expansion is an approximation to $U_\epsilon$, the identification of $A$ as the coefficent of the first-order term is not.
$^1$It's important to note that for a generic $\psi$, this difference quotient does not necessarily converge to a well-defined limit; in such cases, the domain of $A$ is defined to be those $\psi$ in the Hilbert space such that the difference quotient does have a well-defined limit.
$^2$In Hamiltonian mechanics, the momentum observable is the infintesimal generator of spatial translations (see e.g. here). Therefore, we define the momentum operator in quantum mechanics to be the infinitesimal generator of spatial translations by way of analogy.