Ttyjfyk

What specific method does Interpolation use for unstructured multi-dimensional data when we set InterpolationOrder -> All? Documentation links are welcome.

Example 2D data:

data = RandomReal[1, 20, 3];

When the data points are not on a grid, the only allowed settings for InterpolationOrder are 1 and All, according to the error message issued when trying something else.

With 1, it is clear how it works: a Delaunay triangulation is computed and linear interpolation is done over each triangle.

But how does All work, and what determines the actual order that is chosen?

if = Interpolation[data, InterpolationOrder -> All];

if["InterpolationOrder"]
(* 5 *)

Show[
 Plot3D[if[x, y], x, 0, 1, y, 0, 1],
 Graphics3D[PointSize[Large], Point[data]]
]

This is code that has been written many moons ago... first an example:

d = 0.4138352728412389, 0.02365673668161028, 0.5509946389658635, 
 0.7254061374370833, 0.14521595926324116, 
 0.6528630823305817, 0.48768962246740544, 
 0.22066264105073286, 0.8309710560928056, 
 0.3496966364384875, 0.4553589220242207, 
 0.9383446951847001, 0.2126873262146789, 
 0.017512080396716145, 0.967248982535015, 
 0.6211273372083488, 0.3548669163916416, 
 0.737108322193581, 0.6919974835480842, 0.9322403408098401;
f = 0.9953617542392983, 0.14070666511222818, 
0.285662339441511, 0.7988192898854105, 0.3592646208757597, 
0.565455746009103, 0.22110814761432618, 0.2735048548887764, 
0.08792348530403005, 0.4202942851818514;
data = Join[d, f, 2];
if = Interpolation[data, InterpolationOrder -> All];
if[0.5, 0.5]

0.268157

And here is roughly what it does:

dt = Transpose[d];
temp = Join[ConstantArray[1., Length[d]], dt, dt[[1]]^2, 
 dt[[1]]*dt[[2]], dt[[2]]^2, dt[[1]]^3, 
 dt[[1]]^2*dt[[2]], dt[[1]]*dt[[2]]^2, dt[[2]]^3];
p = Transpose[temp];
ls = LinearSolve[p];
vals = ls[Flatten[f]];
System`Private`EvaluateListPolynomial[vals, 0.5, 0.5]

0.268157