# Meaning of InterpolationOrder -> All for multidimensional interpolation

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


• Dunno, but the return value of if["InterpolationOrder"] that I get is {9223372036854775806, 9223372036854775806}. Oo Commented Apr 16, 2019 at 8:49
• @HenrikSchumacher Oops ... It seems I tried this with M12.0 (it's available in the cloud). Commented Apr 16, 2019 at 8:54
• Anyways, very good questions. I am also curious what works there in the background. Commented Apr 16, 2019 at 9:00
• That sounds as if they were using straight-forward global interpolation by a polynomial of degree up to n. Then you have Binomial[n, 2] basis functions. In that case, this should become nasty for higher point counts due to Runge's phenomenon and ill-conditioned linear systems (for solving for the coefficients). So I presume, that they will switch to another method when the point count becomes larger... Commented Apr 16, 2019 at 9:06
• In this context, I just want to point to the "Obtuse" package for doing multidimensional interpolation on scattered data points. You can download it (for free) from familydahl.se/mathematica. The package implements five different methods: the Delaunay interpolation method, the Voronoi interpolation method, the Shepard interpolation method, the RBF interpolation method and the ObtuseAngle interpolation method. For the ObtuseAngle method the InterpolationOrder can be set to 0, 1, 2, and 3, but not to All. The number of dimensions is not limited. Commented Jun 12, 2019 at 19:54

This is code that has been written many moons ago... here is 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]];
SystemPrivateEvaluateListPolynomial[vals, {0.5, 0.5}]

0.268157