Let $f(x,y):\mathbb R^2\to\mathbb R$ be convex, and let $(x_0,y_0)$ and $(x_1,y_1)$ be distinct points. Define $g(x,y)=a x y+b x+cy+d$ to be the unique multilinear function that agrees with $f(x,y)$ in $(x_0,y_0)$ and $(x_1,y_1)$ and such that $\nabla g$ agrees with $\nabla f$ in those points as well. (This defines 6 equationis in $a,b,c,d$. We only care about the cases where this system is solvable.)
Can we say that $f(x,y) \ge g(x,y)$ for all $(x,y)\in\mathbb R^2$?
If $f$ was a one-argument convex function $\mathbb R\to \mathbb R$, we would know that if $f$ agreed with a linear function in a point, and also agreed with its derivative, then $f(x)\ge g(x)$ for all $x\in\mathbb R$.
I'm wondering if the same thing holds in two dimensions? Or if not, are there any further requirements we can put on $f$ such that this global condition holds?
In particular, I'm interested in whether this holds for the function $f(x,y) = h(x)h(y)$ where $h(x)=((1+x)^s+(1-x)^s)^{1/s}$ for some $s\in[1,2]$.
Update: Ok, a counterexample is $f(x,y)=|x+y|$, which agrees with $g(x,y)=x y+1$ at $(x,y)=(-1,-1)$ and $(1,1)$, but clearly $xy+1>|x+y|$ for $x$ and $y$ large enough.
The question should then really be if products of convex functions (that is $f(x,y)=h_1(x)h_2(y)$) have the property mentioned above.
Update 2: So another counterexample is $f(x,y)=|x|y$ which of course agrees with $xy$ on all of $x\ge 0$, but $|x|y\le xy$ for $x\le0,y\le 0$.
It still appears to hold for my example function, but maybe we also need to assume that $h_1=h_2$ or that the functions are non-negative.