This answer is a small modification of the answer of Denis Serre. I added for reader's convenience: (1) the result is slightly more general; (2) the answer contains much more details; (3) I am using a convolution by mollification approximation instead of inf-convolution.

Since convex functions satisfy
$$
\langle \nabla f(x)-\nabla f(y),x-y\rangle\geq 0,
$$
it suffices to prove the following more general result.

**Theorem.**
Let $f\in C^1(\mathbb{R}^n)$ and let $L>0$.Then the following conditions are equivalent:
\begin{equation}
(1)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ |\langle\nabla f(x)-\nabla f(y),x-y\rangle|\leq L|x-y|^2
\quad
\text{for all $x,y\in\mathbb{R}^n$.}
\end{equation}
\begin{equation}
(2)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ |\nabla f(x)-\nabla f(y)|\leq L|x-y|
\quad
\text{for all $x,y\in\mathbb{R}^n$.}
\end{equation}

**Proof.**
While the implication (2) to (1) is obvious the other is not so we will prove the implication from (2) to (1) now.
Assume first that $f\in C^\infty(\mathbb{R}^n)$. For $|u|=1$, (1) yields,
$$
\left|\left\langle\frac{\nabla f(x+tu)-\nabla f(x)}{t},u\right\rangle\right|\leq L,
$$
so passing to the limit as $t\to 0$ gives
$$
|\langle D^2f(x)u,u\rangle|\leq L.
$$
Since $D^2 f(x)$ is a symmetric matrix, the spectral theorem implies that the operator norm of the matrix $D^2f(x)$ satisfies
$$
\Vert D^2f(x)\Vert = \sup_{|u|=1}|\langle D^2f(x)u,u\rangle|\leq L.
$$
This estimate however, easily implies the result
\begin{equation}
\begin{split}
& |\nabla f(x)-\nabla f(y)|=
\left|\int_0^1\frac{d}{dt}\nabla f(y+t(x-y))\, dt\right|\\
&\leq |x-y|\int_0^1\Vert D^2f(y+t(x-y))\Vert\, dt\leq L|x-y|.
\end{split}
\end{equation}
This completes the proof when$f\in C^\infty$. Assume now that $f\in C^1$ and let $f_\epsilon=f*\varphi_\epsilon$ be a standard approximation by convolution. Recall that $f_\epsilon\in C^\infty$ and $\nabla f_\epsilon=(\nabla f)*\varphi_\epsilon$. We have
\begin{equation}
\begin{split}
&
|\langle \nabla f_\epsilon(x)-\nabla f_\epsilon(y),x-y\rangle|=
\Big|\Big\langle\int_{\mathbb{R}^n} (\nabla f(x-z)-\nabla f(y-z))\varphi_\epsilon(z)\, dz,x-y\Big\rangle\Big|\\
&\leq
\int_{\mathbb{R}^n}
\big|\big\langle \nabla f(x-z)-\nabla f(y-z)),(x-z)-(y-z)\big\rangle\big|\, \varphi_\epsilon(z)\, dz
\leq L|x-y|^2,
\end{split}
\end{equation}
where the last inequality is a consequence of (1) and $\int_{\mathbb{R}^n}\varphi_\epsilon=1$.
Since $f_\epsilon\in C^\infty$, the first part of the proof yields
$$
|\nabla f_\epsilon(x)-\nabla f_\epsilon(y)|\leq L|x-y|
$$
and the result follows upon passing to the limit as $\epsilon\to 0$.