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!
 A: Let 
$$
h(x)=\left\{\begin{array}{lll}\mathrm{e}^{-1/x^2} & \text{if} & x>0,\\
0 & \text{if} & x\le 0.\end{array}\right.
$$
Then $h\in C^\infty(\mathbb R)$. Then set
$$
j(\boldsymbol{x})=c\,h\big(1-\|\boldsymbol{x}\|^2\big),
$$
where $\boldsymbol{x}\in\mathbb R^n$, and $c>0$, so that $\int_{\mathbb R^n}j(\boldsymbol{x})\,d\boldsymbol{x}=1$. Clearly, $j\ge 0$, $j\in C^\infty(\mathbb R^n)$
and $\,\mathrm{supp}\,j\subset B(0,1)$ - the unit ball.
Next define $j_e(\boldsymbol{x})=\varepsilon^{-n}j(\varepsilon^{-1}\boldsymbol{x})$. Then $\int_{\mathbb R^n}j_\varepsilon(\boldsymbol{x})\,d\boldsymbol{x}=1$ and let the function
$$
f_\varepsilon=f*j_\varepsilon,
$$
i.e.,
$$
f_\varepsilon(\boldsymbol{x})=\int_{\mathbb R^n} f(\boldsymbol{y})\,j_\varepsilon(\boldsymbol{x}-\boldsymbol{y})\,d\boldsymbol{y}=\int_{\mathbb R^n} f(\boldsymbol{x}-\boldsymbol{y})\,j_\varepsilon(\boldsymbol{y})\,d\boldsymbol{y}=
\frac{1}{\varepsilon^n}\int_{B(0,\varepsilon)} f(\boldsymbol{x}-\boldsymbol{y})\,j(\boldsymbol{y}/\varepsilon)\,d\boldsymbol{y}=\frac{1}{\varepsilon^n}\\=\int_{B(0,1)} f(\boldsymbol{x}-\varepsilon\boldsymbol{y})\,j(\boldsymbol{y})\,d\boldsymbol{y}.
$$
Clearly $f_\varepsilon\in C^\infty(\mathbb R^n)$. Next
$$
f_\varepsilon(\boldsymbol{x}_1)-f_\varepsilon(\boldsymbol{x}_2)=
\int_{\mathbb R^n} \big(f(\boldsymbol{x}_1-\boldsymbol{y})-f(\boldsymbol{x}_2-\boldsymbol{y})\big)\,j_\varepsilon(\boldsymbol{y})\,d\boldsymbol{y}
$$
and hence
$$
\lvert\,f_\varepsilon(\boldsymbol{x}_1)-f_\varepsilon(\boldsymbol{x}_2)\rvert\le
\int_{\mathbb R^n} \lvert\, f(\boldsymbol{x}_1-\boldsymbol{y})-f(\boldsymbol{x}_2-\boldsymbol{y})\rvert\,j_\varepsilon(\boldsymbol{y})\,d\boldsymbol{y}\le\|\boldsymbol{x}_1-\boldsymbol{x}_2\|\int_{\mathbb R^n}j_\varepsilon(\boldsymbol{y})\,d\boldsymbol{y}=\|\boldsymbol{x}_1-\boldsymbol{x}_2\|.
$$
Finally
$$
\lvert\,f_\varepsilon(\boldsymbol{x})-f(\boldsymbol{x})\rvert\le
\left|\int_{B(0,1)} \big(f(\boldsymbol{x}-\varepsilon\boldsymbol{y})-f(\boldsymbol{x})\big)\,j(\boldsymbol{y})\,d\boldsymbol{y}.\,\right|\le \cdots\le \varepsilon.
$$
A: 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.
