Invert a $4 \times 4$ matrix with a given structure?

Question

When tryig to fit $f(x,y) = a+bx+cy+dxy$ to the values of four points, we will have to invert following matrix. (Let us assume that $x_i,y_i$ are chosen suitably such that it is regular.) We could obviously solve this with any method (LU etc), but does the special structure of this matrix maybe allow for a compact representation of its inverse?

I thought there might be a way because this matrix looks similar to a Vandermonde matrix (but it is obviously not the same), which do sometimes have inverses that are particularly easy to compute.

$$M=\begin{bmatrix} 1 & x_1 & y_1 & x_1 y_1\\ 1 & x_2 & y_2 & x_2 y_2\\ 1 & x_3 & y_3 & x_3 y_3\\ 1 & x_4 & y_4 & x_4 y_4\end{bmatrix}$$

Answer 1

Given four distinct points $(x_1,y_1),(x_2,y_2),(x_3,y_3),(x_4,y_4)$, we can exhibit the cases for which the matrix $M$ is invertible, and hence for which all the bilinear interpolation problems are uniquely solvable, by evaluating the determinant of $M$ as a quartic polynomial in the $x_i,y_i$.

It will simplify our work, and help to provide a geometric interpretation of the result, if we begin with a change of variable translating the first of these points to the origin:

$$ u = x - x_1 $$ $$ v = y - y_1 $$

Supposing that with the same function values at points $(u_i,v_i) = (x_i-x_1,y_i-y_1)$ as in the corresponding original points, a bilinear interpolant exists:

$$ g(u,v) = \hat a + \hat b u + \hat c v + \hat d uv $$

Then the solution to the original problem will be:

$$ \begin{align*} f(x,y) &= \hat a + \hat b (x-x_1) + \hat c (y-y_1) + \hat d (x-x_1)(y-y_1) \\ &= (\hat a - \hat b x_1 - \hat c y_1 + \hat d x_1 y_1) + (\hat b - \hat dy_1)x \\ &\; \; \; + (\hat c - \hat dx_1)y + \hat d xy \end{align*} $$

It follows that this solution would uniquely determine the desired bilinear interpolant $f(x,y)$ at the original points:

$$ \begin{pmatrix} a \\ b \\ c \\ d \end{pmatrix} = \begin{pmatrix} 1 & -x_1 & -y_1 & x_1y_1 \\ 0 & 1 & 0 & -y_1 \\ 0 & 0 & 1 & -x_1 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \hat a \\ \hat b \\ \hat c \\ \hat d \end{pmatrix} $$

In that sense we can assume without loss of generality that $(x_1,y_1) = (0,0)$, and we proceed to find the determinant of $M$ in this slightly special case:

$$ \begin{align*} \det M &= \left| \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 1 & x_2 & y_2 & x_2 y_2 \\ 1 & x_3 & y_3 & x_3 y_3 \\ 1 & x_4 & y_4 & x_4 y_4 \end{array} \right| \\ &= \left| \begin{array}{ccc} x_2 & y_2 & x_2 y_2 \\ x_3 & y_3 & x_3 y_3 \\ x_4 & y_4 & x_4 y_4 \end{array} \right| \\ &= x_2x_4y_3(y_4-y_2) + x_3x_4y_2(y_3-y_4) + x_2x_3y_4(y_2-y_3) \end{align*} $$

With this expression we can identify some geometric conditions on the four points for which the matrix has zero determinant and is thus not invertible. One such case is for all four points to be collinear. For example, $(0,0),(1,1),(2,2),(3,3)$. Proof: If all four points are collinear, then the first two columns of the $3\times 3$ determinant above are linearly dependent.

Another case is for two points to have the same $x$-coordinate while the other two points share the same $y$-coordinate. For example, $(0,0),(0,1)$ and $(2,1),(2,2)$. Proof: Suppose WLOG that the two points with equal $x$-coordinates are $x_1=x_2=0$ as previously argued. Then $y_3=y_4$ and all the summands in our expression for $|M|$ vanish.

So we can give a formula for the $4\times 4$ inverse exactly when determinant of $M$ is not zero. Putting the expression for the determinant back in its more general form, not assuming $(x_1,y_1)=(0,0)$, makes it look a bit more similar to the Vandermonde matrix determinant:

$$ \begin{align*} \det M &= (x_2-x_1)(x_4-x_1)(y_3-y_1)(y_4-y_2) \\ &\; \; \; + (x_3-x_1)(x_4-x_1)(y_2-y_1)(y_3-y_4) \\ &\; \; \; + (x_2-x_1)(x_3-x_1)(y_4-y_1)(y_2-y_3) \end{align*} $$

Answer 2

Obtaining coefficients $a, b, c, d$ of the surface (quadric) with equation

$$z=a+bx+cy+dxy$$

passing through 4 points doesn't necessitate to solve a system of equations or the inversion of a linear system. It can be obtained directly. Here is how (Meanwhile, we will deepen our understanding because, in this way, one will get the inverse matrix of the initial matrix) :

Let us take the simplified case where the four points are

$$(x,y) = (0,0), \ \ (1,0), \ \ (0,1), \ \ (1,1)$$

with prescribed values :

$$f(0,0)=p, \ f(1,0)=q, \ f(0,1)=r, \ f(1,1)=s. \ \ \ \ (1)$$

Let us consider the following expression : in it, each one of the 4 terms has a specialized "task" ensuring that the four constraints (1) are fulfilled :

$$ \begin{equation} z=\left\{\begin{array}{llll} p&\times& \ (1-x)(1-y) \ +... \ & \text{ ``activated'' iff} \ x=0 \ \& \ y=0\\ q&\times& \ x(1-y) \ +... \ & \text{ ``activated'' iff} \ x=1 \ \& \ y=0\\ r&\times& \ (1-x)y \ +... \ & \text{ ``activated'' iff} \ x=0 \ \& \ y=0\\ s&\times& \ xy; \ & \text{ ``activated'' iff} \ x=1 \ \& \ y=1 \end{array}\right. \end{equation} $$

Expanding all, we get :

$$z=p(1-x-y+xy)+q(x-xy)+r(y-xy)+sxy \ \ \iff$$

$$z=\underbrace{p}_{= \ a}+\underbrace{(-p+q)}_{= \ b}x+\underbrace{(-p+r)}_{= \ c}y+\underbrace{(p-q-r+s)}_{= \ d}xy$$

or, with a matrix formulation :

$$\begin{pmatrix}a\\b\\c\\d\end{pmatrix}= \begin{pmatrix} \ \ 1& \ \ 0& \ \ 0&0\\ -1& \ \ 1& \ \ 0&0\\ -1& \ \ 0&\ \ 1&0\\ \ \ 1&-1&-1&1 \end{pmatrix} \begin{pmatrix}p\\q\\r\\s\end{pmatrix}$$

which is not surprizing, being naturally the inverse of system

$$\begin{pmatrix}p\\q\\r\\s\end{pmatrix}= \begin{pmatrix} 1&x_1&y_1&x_1y_1\\ 1&x_2&y_2&x_2y_2\\ 1&x_3&y_3&x_3y_3\\ 1&x_4&y_4&x_4y_4\\ \end{pmatrix} \begin{pmatrix}a\\b\\c\\d\end{pmatrix} \ \ \iff \ \ $$

(with $(x_1,y_1)=(0,0), \ \ (x_2,y_2)=(1,0), \ \ (x_3,y_3)=(0,1), (x_4,y_4)=(1,1))$

$$\begin{pmatrix}p\\q\\r\\s\end{pmatrix}= \begin{pmatrix} 1&0&0&0\\ 1&1&0&0\\ 1&0&1&0\\ 1&1&1&1 \end{pmatrix} \begin{pmatrix}a\\b\\c\\d\end{pmatrix}.$$

---

Edit :

1) This method (called "bilinear interpolation" as noted by @hardmath) as the OP has make the remark provides an answer only in the case of points constituting a rectangle with sides parallel to the axes.

2) If we assume that the coordinates axes constitute an orthonormal basis, the determinant of the matrix given in the question can receive a geometrical interpretation when expanded along its fourth column :

$$-x_1y_1\underbrace{\begin{vmatrix} 1&x_2&y_2\\ 1&x_3&y_3\\ 1&x_4&y_4 \end{vmatrix}}_{2 \ \times \ \text{area of} \ P_2P_3P_4} \ + \ x_2y_2\underbrace{\begin{vmatrix} 1&x_1&y_1\\ 1&x_3&y_3\\ 1&x_4&y_4 \end{vmatrix}}_{2 \ \times \ \text{area of} \ P_1P_3P_4} \ + \ \cdots$$

thus a weighted combination of (oriented) areas of the different triangles one can make with points $P_1, P_2, P_3, P_4$ the weights being the product of coordinates of these points.

As one can freely choose the coordinate axes, one can take axes with $P_1$ the origin, and axes with basis $\overrightarrow{P_1P_2}, \overrightarrow{P_1P_3}$ (under the condition that $P_1,P_2,P_3$ aren't aligned ; this is of course no longer an orthonormal basis in general!). Thus, the determinant reduces to :

$$x_4y_4\begin{vmatrix} 1&x_1&y_1\\ 1&x_2&y_2\\ 1&x_3&y_3 \end{vmatrix}$$

which is nonzero unless point $P_4$ is situated on one of the coordinate axes (i.e. aligned with $P_1$ and $P_2$ or aligned with $P_1$ and $P_3$).