How to approximate a globally Lipschitz function by differentiable functions with bounded derivatives? for some positive integer $d \geq 1$
I have a globally Lipschitz continuous function $f \colon \mathbb{R}^d \to \mathbb{R}$ with Lipschitz constant $1$ and would like to approximate it by a sequence $f_\varepsilon$ with the following properties:


*

*For all $\varepsilon > 0$, the function $f_\varepsilon$ is $k$ times partially differentiable and all its partial derivatives up to order $k$ are bounded. Here, $k$ is some positive integer ($k=3$ suffices for my purposes but arbitrary $k$ would be nicer);

*For all $\varepsilon > 0$, it holds that $\| f_\varepsilon \|_{\text{Lip}} \leq \| f \|_{\text{Lip}}$, where $\|\cdot\|_\text{Lip}$ denotes the usual Lipschitz norm;

*For $\varepsilon \to 0$, it holds that $\| f_\varepsilon - f \|_{\infty} \to 0$.


For $d=1$ and $k=2$ I have found the following example in the literature which is stated without proof (it should continue to work for arbitrary $d$ if one replaces the one dimensional Gaussian distribution by a $d$-dimensional one) but don't understand why the second derivative exists (in the classical sense) and is bounded: $$f_\varepsilon(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} g(x + \varepsilon y) e^{- y^2/2} \, \mathrm{d}y.$$ 
Can anybody give me a hint or provide an example of such a sequence?
Thanks in advance!

Rewrite $f_{\epsilon}(x)$ by substituting $y=y'-\frac{1}{\epsilon}x$:
$$
\begin{align}
    f_{\epsilon}(x) & = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}g(\epsilon y')
             e^{-(y'-\frac{1}{\epsilon}x)^{2}/2}\,dy' \\
        & = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}g(\epsilon y')e^{-(\epsilon y'-x)^{2}/2\epsilon^{2}}\,dy'
\end{align}
$$
Now you can see that the derivatives respect to $x$ fall onto the Gaussian, instead onto the function $g$. You can further substitute $y''=\epsilon y$ to make the integral a little more transparent.
