Given the surface $$z=f(x,y)=\begin{bmatrix}1-x&x\end{bmatrix}\begin{bmatrix}b&c\\a&d\end{bmatrix}\begin{bmatrix}1-y\\y\end{bmatrix}$$ defined on $x,y\in[0,1]$, (where $a, b, c$ and $d$ are parameters that equal $f(1, 0), f(0, 0), f(0, 1)$ and $f(1,1)$ respectively), how can I find analytically the point $\mathbf x_c=\left(x_c,y_c,f(x_c,y_c)\right)$ closest to some given point $\mathbf x_b=(x_b,y_b,z_b)$?
Here's what I've tried so far:
If you define a function $$D^2(x,y)=(x_b-x)^2+(y_b-y)^2+(z_b-f(x,y))^2$$ then $\mathbf x_c$ will be the value that minimizes this function -- that is, $$\frac{dD^2}{dx}(x_c,y_c)=0\\ \frac{dD^2}{dy}(x_c,y_c)=0$$
This expands to:
$$0=-2(x_b-x_c)-2(z_b-f)\frac{df}{dx}\\ 0=-2(y_b-y_c)-2(z_b-f)\frac{df}{dy}$$
The problem is, when you actually to solve for $x_c$ and $y_c$, the $f$ and $\frac{df}{dk}$ terms expand into such a mess of constants that it becomes nearly impossible to solve by hand. Wolfram|Alpha also refuses to get anywhere, although I don't know if that's because of the number of terms or because it's actually impossible to solve.
I thought I might be able to solve it using the fact that $\mathbf x_b-\mathbf x_c$ will always be normal to $f$, but that ends up as the same equation:
$$\begin{bmatrix}-\frac{df}{dx}\\ -\frac{df}{dy} \\ 1\end{bmatrix}\text{is a vector normal to }f\\ \mathbf x_b - \mathbf x_c=(\text{some scalar }s)\begin{bmatrix}-\frac{df}{dx}\\ -\frac{df}{dy} \\ 1\end{bmatrix}\\ x_b-x_c=-s\frac{df}{dx}\\ y_b-y_c=-s\frac{df}{dy}\\ z_b-f(x_c,y_c)=s\\ \Rightarrow x_b-x_c=-(z_b-f)\frac{df}{dx}\\ \Rightarrow y_b-y_c=-(z_b-f)\frac{df}{dy}$$
Does anyone know of a way to simplify the algebra, or a proof that it cannot be solved exactly? (I've already found some approximations that will work fine for what I need, so I'm not interested in iterative methods; I'm just curious at this point)