interpolation - How does ListInterpolation work?

ListInterpolation takes an array of data and interpolates between the entries. Let us interpolate between the elements of the identity matrix. For reference, this is what the identity matrix looks like:

MatrixPlot[IdentityMatrix[5]]

Now let's interpolate, and request first order interpolation:

if = ListInterpolation[IdentityMatrix[5], {{0, 1}, {0, 1}}, InterpolationOrder -> 1];
Plot[if[x, x], {x, 0, 1}]

The result is a piecewise function as expected, but not a piecewise linear function! Why? What does first order interpolation mean here?

---

The kind of linear interpolation I am familiar with is what ListDensityPlot does:

ListDensityPlot[IdentityMatrix[5], DataRange -> {{0, 1}, {0, 1}}, Mesh -> All]

Construct a Delaunay triangulation of the data and interpolate on each triangle.

This is not what ListInterpolation does (or even Interpolation when the data is strictly on a square grid). What does ListInterpolation do, then?

---

Inspired by this question:

Answer

I'll just treat the simplest case of interpolating across four noncoplanar points. The extension to general grids is straightforward (one just joins multiple pieces appropriately).

Consider the following bivariate interpolating polynomial:

ip[x_, y_] = InterpolatingPolynomial[{{{0, 0}, 1}, {{0, 1}, -1/2},
                                      {{2, 0}, 1/2}, {{2, 1}, 0}}, {x, y}];

and the following InterpolatingFunction[]:

if = ListInterpolation[{{1, -1/2}, {1/2, 0}}, {{0, 2}, {0, 1}}, InterpolationOrder -> 1];

They are the same:

Plot3D[ip[x, y] - if[x, y] // Chop, {x, 0, 2}, {y, 0, 1}]

Let's have a look at what ip[] looks like:

Expand[ip[x, y]]

   1 - x/4 - (3 y)/2 + (x y)/2

Huh. That last term certainly isn't linear. What gives?

In a bilinear interpolation such as this one, although the procedure is something along the lines of "linearly interpolate for each grid line in one direction, and then linearly interpolate those results", the result is inevitably quadratic. This should come as no surprise to people familiar with the hyperbolic paraboloid, which is one of the simplest examples of a ruled surface, or a surface formed by sweeping out a line in space. In fact, z == 1 - x/4 - (3 y)/2 + (x y)/2 is indeed a hyperbolic paraboloid.

{Plot3D[ip[x, y], {x, 0, 2}, {y, 0, 1}, MeshFunctions -> {#1 - #2 &}], 
 Plot[ip[x, x], {x, 0, 1}]} // GraphicsRow