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

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}$$

• This is bilinear interpolation. – hardmath Jan 7 '19 at 22:36
• @hardmath Yes, but unfortunately WP also only shows the cases of an unit square and an axes-aligned rectangle. – flawr Jan 7 '19 at 23:07
• I take your points! A more general solution can be found if the four points are an affine image of corners of the unit square. For some four points the determinant vanishes and no matrix inverse exists. I'll add a note. – hardmath Jan 7 '19 at 23:34
• @hardmath The affine image of a rectangle is a parallelogram; thus we still haven't an answer for the general case where the 4 points can be any set of four points. – Jean Marie Jan 8 '19 at 5:55

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*}

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 :

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

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

• First of all, thank you for your extensive answer. I see how this simplified result can easily be extended to the case where the four points form an arbitrary rectangle that is aligned with the axes, by modifying the four constraints accordingly. For $(x,y) = (u,v),(u',v),(u',v'),(u,v')$ we'd get $z = p\times (1-x/u)(1-y/v)+ q \times x/u'(1-y/v) + \ldots$, which we could expand in the same way as you did to get the coefficients. But I fail to see how this provides a solution (or inverse of the matrix) for the general case. Could you elaborate on that? – flawr Jan 7 '19 at 23:05
• I agree on the fact (as I just did in my comment to @hardmath) that bilinear interpolation doesn't cover the general case (when I gave my answer, I recognize that I had completely overlooked this aspect). I just computed the (formal) inverse of the general matrix you give : it is very complicated... I think to two ways to treat the general case : 1) the usage of projective geometry that can provide a correspondence between any quadrilateral and a rectangle or 2) the usage of B-splines that can provide a way analogous to bilinear interpolation, providing a way to specify a local influence – Jean Marie Jan 8 '19 at 6:15
• Similar questions are treated in the very interesting (rather old) book "Group-Theoretical Methods in Image Understanding" by Ken-ichi Kanatani (Springer, 1990) – Jean Marie Jan 8 '19 at 6:27
• I have added an edit to my answer. – Jean Marie Jan 8 '19 at 7:08