Two seemingly equivalent formulas for bilinear interpolation yield different results

Question

Suppose that I have a function $f:[0,1]^2 \to \mathbb{R}$ and I would like to approximate it from its values in the corners of the domain by the well-known formula $$ f(s,t) \approx (1-s,s) \begin{pmatrix} f(0,0) & f(0,1) \\ f(1,0) & f(1,1) \end{pmatrix} \begin{pmatrix} 1-t\\ t \end{pmatrix} \enspace, $$ which can be found on Wikipedia.

I tried two ways of writing this in Mathematica:

bilinearExpl[f_, s_, t_] :=
        (1 - s) (1 - t) f[0, 0] + (1 - s) t f[0, 1] + (1 - t) s f[1, 0] + s t f[1, 1]
bilinearMatr[f_, s_, t_] := {1 - s, s}.{{f[0, 0], f[0, 1]}, {f[1, 0], f[1, 1]}}.{1 - t, t}

If I try them out, they produce results that are quite much the same:

bilinearExpl[f, s, t]

gives

whereas

bilinearMatr[f, s, t]

returns

Taking a function

scalarFunction[x_, y_] := x^2 y

as an example, they also produce the same approximation of scalarFunction.

{Plot3D[scalarFunction[s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}],
 Plot3D[bilinearExpl[scalarFunction, s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}],
 Plot3D[bilinearMatr[scalarFunction, s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}]}

So far so good. The same formula should work for a function $f:[0,1]^2 \to \mathbb{R}^2$ and both methods should give the same results. But they don't! Take, for instance,

vectorFunction[x_, y_] := {x, y}

Then

{ParametricPlot[vectorFunction[s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}],
 ParametricPlot[bilinearExpl[vectorFunction, s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}],
 ParametricPlot[bilinearMatr[vectorFunction, s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}]}

gives

and the two right plots are certainly not the same.

Question: why are the results of the two (seemingly equivalent) formulas different?

I suppose there is some dimension mismatch of the arrays.

Answer 1

FWIW, Mathematica knows about bilinear interpolation:

InterpolatingPolynomial[{{{0, 0}, f[0, 0]}, {{1, 0}, f[1, 0]},
                         {{0, 1}, f[0, 1]}, {{1, 1}, f[1, 1]}}, {s, t}]

For some reason, InterpolatingPolynomial[] chokes when f is vector-valued; I'm not sure why.

In any event: you seem to know about Mathematica's symbolic capabilities, so you should use it to help you debug when things are not going the way you wanted them. Here's what happens when you try to bilinearly interpolate a vector-valued function with your code:

With[{f = {#1, #2} &}, (* bilinearExpl *)
     (1 - t) ((1 - s) f[0, 0] + s f[1, 0]) + t ((1 - s) f[0, 1] + s f[1, 1])]
   {s (1 - t) + s t, t}

This is the correct result; effectively, the interpolation formula's scalar form is mapped across your vector function's components.

With[{f = {#1, #2} &}, (* bilinearMatr *)
     {1 - s, s}.{{f[0, 0], f[0, 1]}, {f[1, 0], f[1, 1]}}.{1 - t, t}]
   {s (1 - t), s (1 - t) + t}

See the difference?

The cure in this case is to reorder your dot products:

With[{f = {#1, #2} &},
     {1 - t, t}.({1 - s, s}.{{f[0, 0], f[0, 1]}, {f[1, 0], f[1, 1]}})]
   {s (1 - t) + s t, t}

(Note the use of parentheses so that the matrix-vector product is done first.)

The reordered dot product formula still makes sense even if the vector-valued function has more than two components:

With[{f = {#1 - #2, #1 + #2, #1 #2} &},
     {1 - t, t}.({1 - s, s}.{{f[0, 0], f[0, 1]}, {f[1, 0], f[1, 1]}})]
   {s (1 - t) + (-1 + s) t, s (1 - t) + (1 + s) t, s t}

while your version will choke:

With[{f = {#1 - #2, #1 + #2, #1 #2} &},
     {1 - s, s}.{{f[0, 0], f[0, 1]}, {f[1, 0], f[1, 1]}}.{1 - t, t}]

Dot::dotsh: Tensors {{s, s, 0}, {-1 + s, 1 + s, s}} and {1 - t, t} have
incompatible shapes.
   {{s, s, 0}, {-1 + s, 1 + s, s}}.{1 - t, t}

Answer 2

I think this might help. Dot documentation says:

Now consider this, noting the operational grouping via parentheses:

rules = {a -> {a1, a2}, b -> {b1, b2}, c -> {c1, c2}, d -> {d1, d2}};

{r, s}.({{a, b}, {c, d}} /. rules).{t, u}

({r, s}.{{a, b}, {c, d}}.{t, u}) /. rules

{(a1 r + c1 s) t + (a2 r + c2 s) u, (b1 r + d1 s) t + (b2 r + d2 s) u}

{(a1 r + c1 s) t + (b1 r + d1 s) u, (a2 r + c2 s) t + (b2 r + d2 s) u}

One possible approach:

ClearAll[bilinearMatr]

Block[{f, s, t},
  bilinearMatr[f_, s_, t_] =
    {1 - s, s}.{{f[0, 0], f[0, 1]}, {f[1, 0], f[1, 1]}}.{1 - t, t}
];

{ParametricPlot[vectorFunction[s, t], {s, 0, 1}, {t, 0, 1}, Mesh -> {11, 11}], 
 ParametricPlot[bilinearExpl[vectorFunction, s, t], {s, 0, 1}, {t, 0, 1}, 
  Mesh -> {11, 11}], 
 ParametricPlot[bilinearMatr[vectorFunction, s, t], {s, 0, 1}, {t, 0, 1}, 
  Mesh -> {11, 11}]}

For automation of this Block/Set work-around see:

• How to make a function like Set, but with a Block construct for the pattern names