New Directions Of Interpolation

| 0 Comments

We have spent a few months looking at how we might interpolate between sets of points \(\left(x_i, y_i\right)\), where the \(x_i\) are known as nodes and the \(y_i\) as values, to approximate values of \(y\) for values of \(x\) between the nodes, either by connecting them with straight lines[1] or with cubic curves[2][3].
Last time[4], in preparation for interpolating between multidimensional vector nodes, we implemented the ak.grid type to store ticks on a set of axes and map their intersections to ak.vector objects to represent such nodes arranged at the corners of hyperdimensional rectangular cuboids.
With this in place we're ready to take a look at one of the simplest multidimensional interpolation schemes; multilinear interpolation.

One Dimension, Two Dimension, Three Dimension, More

Figure 1: Rectangular Interpolation
The easiest way to describe this algorithm is to start with the special case of two dimensional nodes, known as bilinear interpolation, in which the nodes sit at the corners of a set of rectangles that divide up the area bounded by the four nodes that sit on the first and last tick marks on both axes.
To interpolate within one of those rectangles we can first choose a pair of opposing edges and draw a line perpendicular to them through the point of interest, then use linear interpolation along both of those edges to approximate the values at their points of intersection with that line and again between them to approximate the value at the point itself, as illustrated by figure 1 which marks the point of interest in red, the nodes in green and the points of intersection in black.

In three dimensions the nodes will sit at the corners of rectangular cuboids or, to put it more plainly, boxes. To interpolate within these we can pick a pair of opposing faces, draw a line perpendicular to them both through the point whose value we're approximating, use bilinear interpolation to approximate the values at the points of intersection of that line with each face and linearly interpolation between those points to approximate the point's value.

Finally, and I wouldn't be at all surprised if you've already guessed where we're going with this, in \(n\) dimensions the nodes sit at the corners of \(n\) dimensional rectangular cuboids which we can interpolate within by picking an opposing pair of \(n-1\) dimensional hyperfaces, drawing a line perpendicular to them through the point, use \(n-1\) dimensional multilinear interpolation to approximate the values at the points of intersection and finish the calculation with linear interpolation between them.

Any Which Way You Choose

One thing that might seem slightly dubious about this algorithm is that the order in which opposing hyperfaces are chosen isn't specified, but it turns out that it doesn't actually matter.

To understand why, first note that multilinear interpolation within an \(n\) dimensional rectangular cuboid is equivalent to a function
\[ f_n\left(x_0,\dots,x_{n-1}\right) = \sum_{i=0}^{2^n-1} a_i \prod_{j=0}^{n-1} x_j^{i_j} \]
for some constants \(a_i\), where \(\sum\) is the summation sign, \(\prod\) is the product sign and \(i_j\) is the \(j^\mathrm{th}\) bit in the binary representation of \(i\) so that
\[ i = \sum_{j=0}^{n-1} 2^{i_j} \]
as proven in derivation 1.

Derivation 1: The Multilinear Interpolation Function
Firstly, since the structure of the formula is symmetric in \(x_i\), so that
\[ \sum_{i=0}^{2^n-1} a_i \prod_{j=0}^{n-1} x_j^{i_j} = \sum_{i=0}^{2^n-1} a^\prime_i \prod_{j=0}^{n-1} x_{p(j)}^{i_j} \]
for any reordering \(p(j)\) of the indices \(j\), if it's correct for any ordering of the axes it must be correct for all of them.

Next, for some ordering of the axes of an \(n+1\) dimensional multilinear interpolation choose the opposing hyperfaces lying on the first \(n\) of them and draw a line through the point \(\left(x_0,\dots,x_{n-1},x_n\right)\) to the points of intersection \(\left(x_0,\dots,x_{n-1},x_{0,n}\right)\) and \(\left(x_0,\dots,x_{n-1},x_{1,n}\right)\). If the formula is correct then their \(n\) dimensional multilinear interpolated values will be
\[ \begin{align*} f_{0,n}\left(x_0,\dots,x_{n-1}\right) &= \sum_{i=0}^{2^n-1} a_{0,i} \prod_{j=0}^{n-1} x_j^{i_j}\\ f_{1,n}\left(x_0,\dots,x_{n-1}\right) &= \sum_{i=0}^{2^n-1} a_{1,i} \prod_{j=0}^{n-1} x_j^{i_j} \end{align*} \]
for some constants \(a_{0,i}\) and \(a_{1,i}\). Linearly interpolating between these yields
\[ \begin{align*} f_{n+1}(x_0,\dots,x_n) &= f_{0,n}\left(x_0,\dots,x_{n-1}\right) + \frac{x_n-x_{0,n}}{x_{1,n}-x_{0,n}} \left(f_{1,n}\left(x_0,\dots,x_{n-1}\right)-f_{0,n}\left(x_0,\dots,x_{n-1}\right)\right)\\ &= \left(\sum_{i=0}^{2^n-1} a_{0,i} \prod_{j=0}^{n-1} x_j^{i_j}\right) + \frac{x_n-x_{0,n}}{x_{1,n}-x_{0,n}} \left(\left(\sum_{i=0}^{2^n-1} a_{1,i} \prod_{j=0}^{n-1} x_j^{i_j}\right) - \left(\sum_{i=0}^{2^n-1} a_{0,i} \prod_{j=0}^{n-1} x_j^{i_j}\right)\right)\\ &= \left(\sum_{i=0}^{2^n-1} a_{0,i} \prod_{j=0}^{n-1} x_j^{i_j}\right) + \frac{x_n-x_{0,n}}{x_{1,n}-x_{0,n}} \sum_{i=0}^{2^n-1} \left(a_{1,i}-a_{0.i}\right) \prod_{j=0}^{n-1} x_j^{i_j}\\ &= \left(\sum_{i=0}^{2^n-1} \left(a_{0,i}-\frac{a_{1,i}-a_{0.i}}{x_{1,n}-x_{0,n}}x_{0,n}\right) \prod_{j=0}^{n-1} x_j^{i_j}\right) + \left(\sum_{i=0}^{2^n-1}\frac{a_{1,i}-a_{0.i} }{x_{1,n}-x_{0,n}}x_n\prod_{j=0}^{n-1} x_j^{i_j}\right)\\ &= \left(\sum_{i=0}^{2^n-1} a^\prime_i x_n^0 \prod_{j=0}^{n-1} x_j^{i_j}\right) + \left(\sum_{i=0}^{2^n-1}a^\prime_{2^n+i} x_n^1\prod_{j=0}^{n-1} x_j^{i_j}\right) = \sum_{i=0}^{2^{n+1}-1} a^\prime_i \prod_{j=0}^n x_j^{i_j} \end{align*} \]
so that the formula will be correct in \(n+1\) dimensions too. Finally, in one dimension we have
\[ f_1\left(x_0\right) = a_0 + a_1 x_0 = a_0 x_0^0 + a_1 x_0^1 = \sum_{i=0}^{2^1-1} a_i \prod_{j=0}^{1-1} x_j^{i_j} \]
and so, by induction, it must be correct in all dimensions.

Secondly, note that the formula has \(2^n\) unknown coefficients which must be uniquely determined by the \(2^n\) nodes and values at the corners of the \(n\) dimensional cuboid and therefore the interpolation must be equivalent to exactly the same function no matter the order in which the hyperfaces are chosen.

Now you've probably noticed that the formula isn't actually linear but is instead an \(n^\mathrm{th}\) order polynomial. We can at least justify the name of the algorithm with the fact that it's conditionally linear in the sense that if we hold all but one of the \(x_i\) constant then it's linear in that one.

Implementing The Algorithm

Given the description of the multilinear interpolation algorithm, it should come as no surprise that the simplest way to implement it is by using recursion, as illustrated for a single hyperdimensional rectangular cuboid by listing 1.

Listing 1: Recursive Multilinear Interpolation
function interpolate(i, x, lb, ub, values) {
 if(i===lb.length) return values;

 return (x[i]-lb[i])/(ub[i]-lb[i]) * interpolate(i+1, x, lb, ub, values[1])
      + (ub[i]-x[i])/(ub[i]-lb[i]) * interpolate(i+1, x, lb, ub, values[0]);
}

Here i is the index of the axis along which we're interpolating between faces, x is an array representing the point whose value we're approximating, lb and ub are arrays of the lower and upper bounds of the cuboid upon each axis and values is a nested array of the values at the corner nodes. For the latter we require that the value associated with a node \(\mathbf{x}_{i_0,i_1,\dots,i_{n-1}}\) should be found at values[i0][i1]...[in-1], where the subscripts \(i_j\) indicate that the node's \(i^\mathrm{th}\) element equals the \(j^\mathrm{th}\) tick on its axis which for a single cuboid will be either the lower or upper bound, indexed with zero and one respectively.
The values at the corner nodes of the opposing faces along the current axis are consequently stored in values[0] and values[1] and so this function simply linearly interpolates between the results of multilinear interpolation within those faces.

Program 1 uses it to interpolate within a multidimensional cuboid with randomly chosen lower and upper bounds having values at its corner nodes determined by a polynomial of the form that we proved is equivalent to multilinear interpolation and so should be exactly represented.

Program 1: The Recursion Versus The Formula

The differences between the interpolation and the polynomial are not inconsistent with those that we might expect from floating point rounding errors and we can therefore be confident that this recursive approach is correct.

For our ak library's implementation we'll represent the nodes with an ak.grid object, as stated in the introduction, and the values either as a nested array, as above, or use a function to calculate them from the ak.vector objects representing the nodes that the former maps arrays of axis tick mark indices to.
Listing 2 gives the ak.multiLinearInterpolate function which you can see first checks that the first argument is an ak.grid object with at least two ticks per axis, all of which are finite, and then checks that the second argument is either an array or a function.

Listing 2: ak.multiLinearInterpolate
ak.multiLinearInterpolate = function(grid, values) {
 var wide = [];
 var axes, dims, lb, dx, i, axis, n, t, interpolate, f;

 if(ak.type(grid)!==ak.GRID_T) {
  throw new Error('invalid grid in ak.multiLinearInterpolate');
 }
 axes = grid.axes();
 dims = axes.length;

 lb = new Array(dims);
 dx = new Array(dims);
 for(i=0;i<dims;++i) {
  axis = axes[i];
  n = axis.length;
  if(n<2) {
   throw new Error('axis too short in ak.multiLinearInterpolate');
  }
  if(!isFinite(axis[0]) || !isFinite(axis[n-1])) {
   throw new Error('non-finite axis in ak.multiLinearInterpolate');
  }
  if(n>2) wide.push({i:i, lb:updateLB(axes[i])});
  lb[i] = 0;
  dx[i] = axis[1]-axis[0];
 }

 switch(ak.nativeType(values)) {
  case ak.ARRAY_T:
   values = copyValues(values);
   break;

  case ak.FUNCTION_T:
   values = mapValues(0, dims, new Array(dims), axes, values);
   break;

  default:
   throw new Error('invalid values in ak.multiLinearInterpolate');
 }
 t = valueType(values);
 validateValues(0, dims, t, grid, values); 

 interpolate = interpolateType(t, dims, values);
 f = function(x) {return interpolate(transformX(x, wide, axes, lb, dx), lb);};

 f.grid = function() {return grid;};
 f.values = function() {return copyValues(values);};
 f.nodes = function() {return makeNodes(dims, grid, values);};
 return Object.freeze(f);
};

As it does so, it saves a copy of the grid axes in axes, the dimension of the interpolation in dims and initialises the lb and dx arrays with a first guess of zero at the indices of the tick marks that will immediately precede the first point whose value we shall approximate and the distances between them and those following them.
It also keeps track of the indices of axes with more than two tick marks, together with a function for updating the guesses at their lower bounds of the cuboids that will contain the next point of interest, in the wide array, initialises the values with the copyValues helper function if they are passed as a nested array and with mapValues if a function is given to calculate them, validates the values and defines an interpolate function based upon their type before finally creating a function to perform the interpolation and assign it property access methods to return the grid defining the nodes, the nested array of values and the set of nodes themselves.

The updateLB helper function uses exactly the same tricks as our one dimensional, or univariate, interpolation algorithms used to efficiently locate the tick marks immediately preceding each element of a vector representing a point upon its associated axis, as illustrated by listing 3.

Listing 3: The updateLB Helper Function
function updateLB(axis) {
 var x0, dx;

 if(axis.length===3) {
  return function(x) {return x<axis[1]?0:1;};
 }

 x0 = axis[0];
 dx = axis[axis.length-1]-x0;

 if(canJump(axis, x0, dx)) {
  return function(x, lb) {return jumpLB(x, axis, x0, dx, lb);};
 }
 return function(x, lb) {return findLB(x, axis, lb);};
}

Specifically, if there are three ticks on an axis then we can locate the lower bound by simply determining whether the element is less or greater than the second.
If there are more, but they are not too unevenly spaced so that
\[ x_{i-1} \leqslant \bigg\lfloor\frac{x-x_0}{x_{n-1}-x_0} \times n\bigg\rfloor \leqslant x_{i+1} \]
for
\[ x_i \leqslant x \leqslant x_{i+1} \]
where the oddly shaped brackets mean the largest integer that is less than or equal to the term between them, then we can use
\[ \bigg\lfloor\frac{x-x_0}{x_{n-1}-x_0} \times n\bigg\rfloor \]
as a first guess at the index of the tick immediately preceding \(x\) and be sure that in the worst case we'll only be off by one in either direction.
The canJump function given in listing 4 determines whether or not this is the case by checking it for the tick marks themselves, since this will mean that it must be for any \(x\) between them, ignoring the first and last because they must be identified as the lower and upper bounds of the first and last pair respectively.

Listing 4: The canJump Helper Function
function canJump(axis, x0, dx) {
 var n = axis.length;
 var i, j;

 for(i=1;i<n-1;++i) {
  j = ak.floor(n*((axis[i]-x0)/dx));
  if(j<i-1 || j>i) return false;
 }
 return true;
}

The indices of the tick marks upon each axis preceding the relevant elements of a point \(\mathrm{x}\) cannot be less than the first or greater than or equal to the last and so we can restrict the search for them to those within that range. Given that we can easily determine whether those elements lie before their current lower bounds or past their current upper bounds we can follow the decision as to whether we need to update them with a jump that restricts the new indices to just one or the other of those extremes as appropriate, as is done in listing 5.

Listing 5: Jumping To The Lower Bound
function jumpLBFwd(x, axis, x0, dx) {
 var n = axis.length;
 var i = Math.min(ak.floor(n*((x-x0)/dx)), n-2);

 if(x<axis[i]) {if(i>0) --i;}
 else if(x>axis[i+1]) {if(i<n-2) ++i;}
 return i;
}

function jumpLBRev(x, axis, x0, dx) {
 var n = axis.length;
 var i = Math.max(ak.floor(n*((x-x0)/dx)), 0);

 if(x<axis[i]) {if(i>0) --i;}
 else if(x>axis[i+1]) {if(i<n-2) ++i;}
 return i;
}

function jumpLB(x, axis, x0, dx, lb) {
 if(x>axis[lb+1]) {
  if(lb<axis.length-2) lb = jumpLBFwd(x, axis, x0, dx);
 }
 else if(x<axis[lb]) {
  if(lb>0) lb = jumpLBRev(x, axis, x0, dx);
 }
 return lb;
}

Failing this, we can simply resort to a binary search through the sorted tick marks on each axis. Listing 6 uses our private by naming convention, ak._unsafeLowerBound function[5] to do so since we know that we are already satisfying the conditions that the public ak.lowerBound would have checked for.

Listing 6: Searching For The Lower Bound
function compare(l, r) {return l-r;}

function findLBFwd(x, axis, lb) {
 if(++lb<axis.length-2 && x>axis[lb+1]) {
  lb = ak._unsafeLowerBound(axis, x, compare, lb+2, axis.length-1)-1;
 }
 return lb;
}

function findLBRev(x, axis, lb) {
 if(--lb>0 && x<axis[lb]) {
  lb = ak._unsafeLowerBound(axis, x, compare, 1, lb)-1;
 }
 return lb;
}

function findLB(x, axis, lb) {
 if(x>axis[lb+1]) {
  if(lb<axis.length-2) lb = findLBFwd(x, axis, lb);
 }
 else if(x<axis[lb]) {
  if(lb>0) lb = findLBRev(x, axis, lb);
 }
 return lb;
}

Here we are exploiting the fact that we know whether the tick that we're looking for is above or below the current lower bound to reduce the range of indices that we search through.

We'll be storing the values associated with the nodes in a nested array and so to initialise it with another nested array we simply need to recursively slice it, its elements, its elements' elements and so on until we get to a non-array element that we can directly copy, as illustrated by listing 7.

Listing 7: Recursively Copying The Values
function copyValues(values) {
 var n;
 if(ak.nativeType(values)===ak.ARRAY_T) {
  n = values.length;
  values = values.slice(0);
  while(n-->0) values[n] = copyValues(values[n]);
 }
 return values;
}

When using a function of the vectors representing the node to calculate their associated values it's still easiest to use recursion construct the nested array of them, only this time we shall store the values of the tick marks as we iterate over them in an array at the index of the current axis, recursing to the next axis as we do so and finally using the fully populated array to construct an ak.vector object to pass to the function once we've run out of axes, as shown by listing 8.

Listing 8: Recursively Calculating The Values
function mapValues(i, dims, x, axes, f) {
 var axis, n, values;

 if(i<dims) {
  axis = axes[i];
  n = axis.length;
  values = new Array(n);
  while(n-->0) {
   x[i] = axis[n];
   values[n] = mapValues(i+1, dims, x, axes, f);
  }
 }
 else {
  values = f(ak.vector(x));
 }
 return values;
}

Once we've initialised the values we need to check that they're consistent with the nodes and that they are all of the same type since we'll be implementing the interpolation for all of the arithmetic types in the ak library that support addition with themselves and multiplication by numbers, just as we did for ak.linearInterpolate. The valueType function given in listing 9 recursively extracts the first element of each of the nested arrays and finally determines the type of the values it contains with ak.type. The validateValues function recursively checks that each nested array has the same length as the number of tick marks upon their associated axis and that they contain values of type t.

Listing 9: Validating The Values
function valueType(values) {
 return (ak.nativeType(values)===ak.ARRAY_T) ? valueType(values[0])
                                             : ak.type(values);
}

function validateValues(i, dims, t, grid, values) {
 var n;
 if(i<dims) {
  n = grid.length(i);
  if(ak.nativeType(values)!==ak.ARRAY_T || values.length!==n) {
   throw new Error('values size mismatch in ak.multiLinearInterpolate');
  }
  while(n-->0) validateValues(i+1, dims, t, grid, values[n]);
 }
 else {
  if(ak.type(values)!==t) {
   throw new Error('values type mismatch in ak.multiLinearInterpolate');
  }
 }
}

To support our own arithmetic types we'll need to use the ak.add and ak.mul overloads to calculate the interpolated values but since this will be inefficient if the values are numbers it makes sense to treat them as a special case using JavaScript's built in operators. The interpolateType function given in listing 10 switches between using the interpolateNumber function for numbers and the interpolateGeneral function for types that have appropriate overloads of ak.add and ak.mul.

Listing 10: Choosing An Implementation
function interpolateType(t, dims, values) {
 var add, mul;

 if(t===ak.NUMBER_T) {
  return function(x, lb) {
   return interpolateNumber(1, 0, dims, values, x, lb);
  };
 }

 try{add=ak.add[t][t]; mul=ak.mul[ak.NUMBER_T][t];} catch(e){}

 if(ak.nativeType(add)!==ak.FUNCTION_T || ak.nativeType(mul)!==ak.FUNCTION_T) {
  throw new Error('non-arithmetic value type in ak.multiLinearInterpolate');
 }

 return function(x, lb) {
  return interpolateGeneral(add, mul, 1, 0, dims, values, x, lb);
 };
}

We shall slightly rearrange the formula for the interpolation from that implemented in listing 1 by passing down the accumulated products of the ratios of the differences between the elements and the tick marks preceding them and the differences between those enclosing them to cut down on the number of multiplication operations required for the calculation, as done by interpolateNumber in listing 11.

Listing 11: Interpolating Numbers
function interpolateNumber(vol, i, dims, values, x, lb) {
 var xi, j;

 if(i===dims) return vol*values;

 xi = x[i]; j = lb[i];
 return interpolateNumber(vol*xi[0], i+1, dims, values[j+1], x, lb)
      + interpolateNumber(vol*xi[1], i+1, dims, values[j],   x, lb);
}

The interpolateGeneral function given in listing 12 performs the same calculation using our ak.add and ak.mul operators for general arithmetic types.

Listing 12: Interpolating General Arithmetic Types
function interpolateGeneral(add, mul, vol, i, dims, values, x, lb) {
 var xi, j;

 if(i===dims) return mul(vol, values);

 xi = x[i]; j = lb[i];
 return add(interpolateGeneral(add, mul, vol*xi[0], i+1, dims, values[j+1], x, lb),
            interpolateGeneral(add, mul, vol*xi[1], i+1, dims, values[j],   x, lb));
}

The transformX function given in listing 13 first checks that the argument of the interpolation is an ak.vector object with the correct number of dimensions, then updates the lower bounds on each axis with more than two tick marks using the lb function that we chose as being the most efficient and, if it has changed, updates it and the width of the interval between it and the next. It then replaces the elements of the array converted from that argument with the pair of ratios of the differences between them and the lower and upper tick marks with the differences between those lower and upper tick marks.

Listing 13: The transformX Helper Function
function transformX(x, wide, axes, lb, dx) {
 var n, i, axis, lbi, dxi, xi;

 if(ak.type(x)!==ak.VECTOR_T) {
  throw new Error('invalid argument in ak.multiLinearInterpolate');
 }

 if(x.dims()!==axes.length) {
  throw new Error('argument size mismatch in ak.multiLinearInterpolate');
 }

 x = x.toArray();
 n = wide.length;
 while(n-->0) {
  i = wide[n].i;
  axis = axes[i];
  lbi = wide[n].lb(x[i], lb[i]);
  if(lbi!==lb[i]) {
   lb[i] = lbi;
   dx[i] = axis[lbi+1]-axis[lbi];
  }
 }
 i = x.length;
 while(i-->0) {
  axis = axes[i];
  lbi = lb[i];
  dxi = dx[i];
  xi = x[i];
  x[i] = [(xi-axis[lbi])/dxi, (axis[lbi+1]-xi)/dxi];
 }
 return x;
}

The makeNodes function given in listing 14 converts the grid and nested array of values into a set of node/value pairs using the former's facility for mapping from arrays of indices to ak.vector objects and the ak.nextState function to iterate over them.

Listing 14: Converting The Grid And Values Into Nodes And Values
function makeNodes(dims, grid, values) {
 var ub = grid.lengths();
 var pos = new Array(dims);
 var n = 1;
 var i, nodes, y;

 for(i=0;i<dims;++i) {pos[i]=0; n*=ub[i];}
 nodes = new Array(n);

 n = 0;
 do {
  y = values;
  for(i=0;i<dims;++i) y = y[pos[i]];
  nodes[n++] = {x:grid.map(pos), y:y};
 }
 while(ak.nextState(pos, ub));
 return nodes;
}

Finally, you can find this implementation in MultiLinearInterpolate.js.

Using The Interpolator

Program 2 demonstrates ak.multiLinearInterpolate by using it to interpolate between nodes with associated values calculated by the formula representing the results of multilinear interpolation, just as we did in program 1.

Program 2: The Interpolator Versus The Formula

To give us a sense of what the results of multilinear interpolation looks like, program 5 uses ak.multiLinearInterpolate to approximate Rosenbrock's function
\[ f\left(x_0, x_1\right) = 100 \times \left(x_1-x_0^2\right)^2 + \left(1-x_0\right)^2 \]
plotting the function itself at the top and the interpolation beneath.

Program 3: Interpolating Rosenbrock's Function

This illustration of the artefacts of multilinear interpolation, which are most clearly revealed along the curved valley floor, is as good a point as any to conclude.

References

[1] Chalk The Lines, www.thusspakeak.com, 2018.

[2] Cubic Line Division, www.thusspakeak.com, 2018.

[3] We're Not For Turning, www.thusspakeak.com, 2018.

[4] Cuboid Space Division, www.thusspakeak.com, 2018.

[5] I Still Haven't Found What I'm Looking For, www.thusspakeak.com, 2017.

Leave a comment