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.
On
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|.$$
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$.