Show that $\nabla f \in C^{1,1}$ if the difference between $f$ and its linear interpolant has a quadratic bound Let $f : \mathbb{R}^d \to \mathbb{R}$ be a differentiable function such that for every $x,y$ and every $\lambda \in [0,1]$: $$|\lambda f(x) + (1-\lambda) f(y) -f(\lambda x + (1-\lambda)y)| \leqslant |x-y|^2 $$ holds.
How can one show that there is a constant $C>0$ such that we have : $$|\nabla f (x) - \nabla f (y)| \leqslant  C |x-y|,$$ for every $x,y$.
I did not find a way to solve this. 
 A: I claim this holds with $C=32$, which can probably be improved.  Translating $f$ and adding an affine function to it does not affect the condition or conclusion. So we may normalize things so that $f(0)=0$ and $\nabla  f(0)=0$, and our goal is to prove $|\nabla f(x)|\le 32|x|$. 
With $\lambda=1/2$, the stated condition becomes
$$\left|\frac{f(x) + f(y)}{2} - f\left(\frac{x+y}{2}\right)\right| \le |x-y|^2,\quad x,y \in\mathbb R^d\tag{1}$$
For the function $g(x)=4|x|^2$ we have 
$$\frac{g(x) + g(y)}{2} - g\left(\frac{x+y}{2}\right) = 2|x|^2+2|y|^2 -|x+y|^2 = |x-y|^2\tag{2}$$
Hence, both functions $h_+=g+f$ and $h_-=g-f$ are midpoint-convex:
$$\frac{h_{\pm}(x) + h_{\pm}(y)}{2} - h_{\pm}\left(\frac{x+y}{2}\right)\ge 0\tag{3}$$
Since Midpoint-Convex and Continuous Implies Convex, both $h_+$ and $h_-$ are convex. In the following, $h$ stands for either $h_+$ or $h^-$. 
By construction, $\nabla h_{\pm}(0)=0$, which in view of their convexity implies $h_{\pm}(x)\ge 0$ for all $x$. Hence $|f(x)|\le 4|x|^2$. This in turn implies $$h_{\pm}(x)\le 8|x|^2\tag{4}$$ 
In the following, $h$ stands for either $h_+$ or $h_-$. 
Fix $x_0$ and let $v = \nabla h(x_0)$. Let $x_1 = x_0+\frac{|x_0|}{|v|}v$.  By convexity, 
$$
h(x_1) \ge h(x_0)+\nabla h(x_0)\cdot (x_1-x_0) \ge |v| |x_0|  \tag{5}
$$
On the other hand, (4) implies
$$
h(x_1) \le 8|x_1|^2 \le 8|2x_0|^2 = 32|x_0|^2 \tag{6}
$$
Combine (5) and (6) to get $|v|\le 32|x_0|$. We thus have proved 
$$
|\nabla h_{\pm}(x)|\le 32|x|,\quad \forall x\in\mathbb R^d \tag{7}
$$
Since $f=(h_+ - h_-)/2$, it follows that $|\nabla f(x)|\le 32|x|$, as claimed. 
Remark
We don't need $f$ to be differentiable everywhere: the inequality $|\nabla f(x)-\nabla f(y)|\le 32|x-y|$ holds provided $f$ is differentiable at the points $x$ and $y$.
A: The desired inequality holds with $C=16$ (which again may not be sharp). Choose a unit vector $\nu$ so that 
$$ |\nabla f(x)-\nabla f(y)|=(\nabla f(x)-\nabla f(y))\cdot\nu.$$
For $n$ large, interpolate a grid between $x$ and $y$, setting  $$h=\frac{|x-y|}n , \qquad z_j = y+jh\frac{x-y}{|x-y|}\quad\mbox{ for  $j=0\ldots n$.}$$
Then the quantity 
$$ \Delta(h) = \frac{f(x+\nu h)-f(x)}h-\frac{f(y+\nu h)-f(y)}h = \sum_{j=1}^n \Delta_j$$
where 
$$
\Delta_j = \frac{f(z_j+\nu h)-f(z_j)}h-\frac{f(z_{j-1}+\nu h)-f(z_{j-1})}h.
$$
Because the midpoint $m_j$ between $z_j+\nu h$ and $z_{j-1}$ is also the midpoint between $z_{j-1}+\nu h$ and $z_j$, 
$$
\frac12 \Delta_j = \frac{f(z_j+\nu h)+f(z_{j-1})-2f(m_j)}{2h}
-\frac{f(z_j)+f(z_{j-1}+\nu h)-2f(m_j)}{2h}.
$$
Applying the assumption with $\lambda=\frac12$ in each case, we find
that because $|z_j-z_{j-1}|=h$,
$$
|\Delta_j|\le \frac2h \left(|\nu h+z_j-z_{j-1}|^2+|\nu h+z_{j-1}-z_j|^2\right) \le \frac4h(2h)^2 = 16h,
$$
for each $j$. By consequence, $\Delta(h)\le 16nh=16|x-y|$,
for each $n$.  Consequently 
$$|\nabla f(x)-\nabla f(y)| = \lim_{n\to\infty}\Delta(h) \le 16|x-y|.$$
