Chalk The Lines

| 0 Comments

Given a set of points \(\left(x_i, y_i\right)\), a common problem in numerical analysis is trying to estimate values of \(y\) for values of \(x\) that aren't in the set. The simplest scheme is linear interpolation, which connects points with consecutive values of \(x\) with straight lines and then uses them to calculate values of \(y\) for values of \(x\) that lie between those of their endpoints.
On the face of it implementing this would seem to be a pretty trivial business, but doing so both accurately and efficiently is a surprisingly tricky affair, as we shall see in this post.

Accurately Interpolating Between Two Points

A mathematical formula for a function representing the straight line that passes through a pair of points \(\left(x_0, y_0\right)\) and \(\left(x_1, y_1\right)\), assuming that \(x_0\) is not equal to \(x_1\), is
\[ f(x) = y_0 + \frac{x-x_0}{x_1-x_0} \times \left(y_1-y_0\right) \]
This follows from the fact that it's a linear function, multiplying \(x\) by the constant \(\frac{y_1-y_0}{x_1-x_0}\) and adding the constant \(y_0-\frac{y_1-y_0}{x_1-x_0} \times x_0\) to the result, and that
\[ \begin{align*} f(x_0) &= y_0 + \frac{x_0-x_0}{x_1-x_0} \times \left(y_1-y_0\right) = y_0 + 0 \times \left(y_1-y_0\right) = y_0\\\\ f(x_1) &= y_0 + \frac{x_1-x_0}{x_1-x_0} \times \left(y_1-y_0\right) = y_0 + 1 \times \left(y_1-y_0\right) = y_1 \end{align*} \]
Unfortunately, whilst this is mathematically correct, it's rather susceptible to the rounding errors of IEEE 754 floating point arithmetic[1]. Specifically, since floating point operations typically introduce multiplicative errors within \(1 \pm \tfrac12 \epsilon\) of the exact results, where \(\epsilon\) is the difference between one and the smallest floating point number greater than one, we have
\[ \begin{align*} f^\ast(x_1) &= \Big(y_0 + \left(y_1-y_0\right) \times \left(1 \pm \tfrac12\epsilon\right)\Big) \times \left(1 \pm \tfrac12\epsilon\right)\\ &= \Big(y_0 + \left(y_1-y_0\right) \pm \left(y_1-y_0\right) \times \tfrac12\epsilon\Big) \times \left(1 \pm \tfrac12\epsilon\right)\\ &= \Big(y_1 \pm \left(y_1-y_0\right) \times \tfrac12\epsilon\Big) \times \left(1 \pm \tfrac12\epsilon\right)\\ &= \Big(y_1 \pm \left(y_1-y_0\right) \times \tfrac12\epsilon\Big) \pm \Big(y_1 \pm \left(y_1-y_0\right) \times \tfrac12\epsilon\Big) \times \tfrac12 \epsilon\\ &= y_1 \pm \left(y_1 - y_0\right) \times \tfrac12 \epsilon \pm y_1 \times \tfrac12 \epsilon + O\left(\epsilon^2\right)\\ &= y_1 \pm \left(y_1 - \tfrac12 y_0\right) \epsilon + O\left(\epsilon^2\right) \end{align*} \]
where \(O\) represents a value of the order of magnitude of its argument. Here we are depending upon the fact that, for finite unequal \(x_0\) and \(x_1\), the fraction \(\frac{x_1-x_0}{x_1-x_0}\) must equal one since even though the numerator and denominator may have rounding errors they will have exactly the same values and any non-zero finite floating point value divided by itself is exactly equal to one.
We can ensure that the floating point result of the function passes exactly through the two points by rearranging it as
\[ f(x) = \frac{x-x_0}{x_1-x_0} \times y_1 + \frac{x_1-x}{x_1-x_0} \times y_0 \]
Now we have
\[ \begin{align*} f(x_0) &= \frac{x_0-x_0}{x_1-x_0} \times y_1 + \frac{x_1-x_0}{x_1-x_0} \times y_0 = 0 \times y_1 + 1 \times y_0 = y_0\\\\ f(x_1) &= \frac{x_1-x_0}{x_1-x_0} \times y_1 + \frac{x_1-x_1}{x_1-x_0} \times y_0 = 1 \times y_1 + 0 \times y_0 = y_1 \end{align*} \]
so that the values of the fractions in both cases will be exactly recovered by floating point arithmetic and consequently so will the values of \(y_0\) and \(y_1\).

Efficiently Finding The Endpoints

If we're interpolating over a set of more than two points then we'll need a way to find the pairs of them that come before and after each point that we wish to approximate. At the very least we should sort the set by increasing \(x_i\), known as nodes in this context, so that we can use our ak.lowerBound[2] algorithm to do so efficiently.
We can slightly improve upon this by keeping track of the last pair of nodes that we used; if the next point that we're approximating also lies between them then we don't need to repeat the search and if it isn't then we can easily check whether we need to look to the left or to the right of them, reducing the number that we need to consider.
Furthermore, if the next point is to the left or to the right of the previous pair then we can first check whether it's between the pair immediately to their left or to their right respectively in order to optimise for the case when it's close to the previous point, which I contend is a reasonably common one.

Now it's not unheard of for the nodes to be equally spaced, in which case we can find the indices of those to the left and right of a point much, much more quickly with
\[ \bigg\lfloor\frac{x-x_0}{x_{n-1}-x_0} \times n\bigg\rfloor\\ \bigg\lfloor\frac{x-x_0}{x_{n-1}-x_0} \times n\bigg\rfloor+1 \]
for \(n\) nodes, where the oddly shaped brackets stand for the largest integer that is less than or equal to the value between them, or at least we could were it not for floating point rounding errors which might throw the result off.
Rather than try to correct for rounding errors numerically, we can deal with them by replacing the nodes with the pair to their left if the point falls below the first and with the pair to their right if it falls above the second. In addition to being extremely simple, a tremendous advantage of this approach is that it doesn't require the nodes to be perfectly evenly spaced but just not so unevenly spaced that it fails to locate the correct nodes, which will certainly be the case if it works for the nodes themselves.

Finally, we must always make sure that we don't try to read past the first or last pair of points in the set.

Not A Number? Not A Problem!

Now we don't actually require that the \(y_i\) are numbers, just that adding them and multiplying them by numbers yield the same type; for example they might be complex numbers, vectors or matrices.
From an implementation perspective we simply need to check that we have appropriate overloads for the ak.mul and ak.add operators[3] if their type is not a Number, as is done by the interpolatePair function given in listing 1.

Listing 1: interpolatePair
function interpolatePair(t) {
 var add, mul;

 if(t===ak.NUMBER_T) {
  return function(x, nodes, lb) {
   var n0 = nodes[lb];
   var n1 = nodes[lb+1];
   return ((x-n0.x)/n0.dx)*n1.y + ((n1.x-x)/n0.dx)*n0.y;
  };
 }

 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 node type in ak.linearInterpolate');
 }

 return function(x, nodes, lb) {
  var n0 = nodes[lb];
  var n1 = nodes[lb+1];
  return add(mul((x-n0.x)/n0.dx, n1.y), mul((n1.x-x)/n0.dx, n0.y));
 };
}

This returns a function that interpolates between nodes at lb and lb+1 using JavaScript's multiplication and addition operators if the type t is a Number and the ak overloads otherwise, provided that they exist. Note that the dx member of each node contains the difference between its \(x\) value and that of the next node which is calculated when we check that the nodes are valid in the prepareNode function given in listing 2.

Listing 2: prepareNodes
function prepareNodes(nodes) {
 var n = nodes.length;
 var i, t;

 if(n<2) throw new Error('too few nodes in ak.linearInterpolate');

 for(i=0;i<n;++i) {
  if(!isFinite(nodes[i].x)) {
   throw new Error('invalid node in ak.linearInterpolate');
  }
 }

 nodes.sort(function(l, r){return l.x-r.x;});

 t = ak.type(nodes[0].y);
 for(i=1;i<n;++i) {
  if(nodes[i].x===nodes[i-1].x) {
   throw new Error('duplicate node in ak.linearInterpolate');
  }
  if(ak.type(nodes[i].y)!==t) {
   throw new Error('node type mismatch in ak.linearInterpolate');
  }
  nodes[i-1].dx = nodes[i].x-nodes[i-1].x;
 }
 return t;
}

Here we first check that there are at least two nodes and that they are finite numbers before sorting them in increasing order along with their associated \(y\) values. Finally we check that no node is equal to the one that precedes it and that their \(y\) values are all of the same type.

Initialising The Interpolator

These functions are used by the ak.linearInterpolate function given in listing 3 after the nodes have been initialised with our usual scheme of dispatching to a constructors object.

Listing 3: ak.linearInterpolate
ak.linearInterpolate = function() {
 var state = {nodes:[], lb:0};
 var arg0  = arguments[0];
 var interPair, f;

 constructors[ak.nativeType(arg0)](state, arg0, arguments);
 interPair = interpolatePair(prepareNodes(state.nodes));

 if(state.nodes.length===2) {
  f = function(x) {
   return interPair(Number(x), state.nodes, 0);
  };
 }
 else if(state.nodes.length===3) {
  f = function(x) {
   x = Number(x);
   return interPair(x, state.nodes, x<state.nodes[1].x ? 0 : 1);
  };
 }
 else if(canJump(state)) {
  f = function(x) {
   return interPair(jumpLB(x, state), state.nodes, state.lb);
  };
 }
 else {
  f = function(x) {
   return interPair(findLB(x, state), state.nodes, state.lb);
  };
 }
 f.nodes = function() {return copyNodes(state.nodes);};
 return Object.freeze(f);
};

var constructors = {};

We shall support construction from an array of nodes and one of their \(y\) values, an array of nodes and a function with which to calculate them or an array of objects having x and y members, as shown by listing 4.

Listing 4: The Constructors
constructors[ak.ARRAY_T] = function(state, nodes, args) {
 var arg1 = args[1];
 constructors[ak.ARRAY_T][ak.nativeType(arg1)](state, nodes, arg1);
};

constructors[ak.ARRAY_T][ak.ARRAY_T] = function(state, x, y) {
 var n = x.length;
 var i;
 if(y.length!==n) throw new Error('node size mismatch in ak.linearInterpolate');
 state.nodes.length = n;
 for(i=0;i<n;++i) state.nodes[i] = {x:Number(x[i]), y:y[i]};
};

constructors[ak.ARRAY_T][ak.FUNCTION_T] = function(state, x, f) {
 var n = x.length;
 var i, xi;
 state.nodes.length = n;
 for(i=0;i<n;++i) {
  xi = Number(x[i]);
  state.nodes[i] = {x:xi, y:f(xi)};
 }
};

constructors[ak.ARRAY_T][ak.UNDEFINED_T] = function(state, nodes) {
 var n = nodes.length;
 var i, x, y;
 state.nodes.length = n;
 for(i=0;i<n;++i) {
  x = nodes[i].x; if(ak.nativeType(x)===ak.FUNCTION_T) x = x();
  y = nodes[i].y; if(ak.nativeType(y)===ak.FUNCTION_T) y = y();
  state.nodes[i] = {x:Number(x), y:y};
 }
};

Now if there are only two nodes then there's no point searching for the pair to interpolate with and so we simply create a function that forwards to the interPair function to calculate the result. If there are three then we can very easily choose between the first and second pair of nodes by comparing the argument to the middle one before calling interPair with the position of the first node in the selected pair.
For more than three nodes we shall directly jump to the correct pair if we can and perform a more computationally expensive search if we can't. The canJump function given in listing 5 tests this by checking whether every node except the first and the last is identified as falling between the pair that has it as an upper bound or as a lower bound.

Listing 5: canJump
function canJump(state) {
 var n = state.nodes.length;
 var x0 = state.nodes[0].x;
 var dx = state.nodes[n-1].x-x0;
 var i, j;

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

Note that the reason why we don't need to check the first and last nodes is because the former is guaranteed to be identified as the lower bound of the first pair of nodes and the latter as the upper bound of the last pair.
If the nodes pass this test then we store the first node and the distance between the first and last nodes in the state object for the benefit of the jump calculation, performed by the jumpLB, jumpLBFwd and jumpLBRev functions given in listing 6.

Listing 6: Jumping To The Correct Pair
function jumpLB(x, state) {
 x = Number(x);
 if(x>state.nodes[state.lb+1].x) {
  if(state.lb<state.nodes.length-2) jumpLBFwd(x, state);
 }
 else if(x<state.nodes[state.lb].x) {
  if(state.lb>0) jumpLBRev(x, state);
 }
 return x;
}

function jumpLBFwd(x, state) {
 var n = state.nodes.length;
 var i = Math.min(ak.floor(n*((x-state.x0)/state.dx)), n-2);

 if(x<state.nodes[i].x) {if(i>0) --i;}
 else if(x>state.nodes[i+1].x) {if(i<n-2) ++i;}
 state.lb = i;
}

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

 if(x<state.nodes[i].x) {if(i>0) --i;}
 else if(x>state.nodes[i+1].x) {if(i<n-2) ++i;}
 state.lb = i;
}

The jumpLB function ensures that we only jump to the left if the argument is to the left of the current lower bound and the current lower bound is greater than zero. Similarly we only jump to the right if it's to the right of the node following the current lower bound and the current lower bound is before the last.
Note that the jumpLBFwd and jumpLBRev functions exploit the fact that the argument must lie to the right of the first node or to the left of the last respectively by only restricting the calculated lower bound node index i in the direction of the jump. After this they both check whether it needs to be decreased or increased to locate the pair of nodes surrounding the argument.

If we can't jump directly to the pair of nodes then we instead perform a search for them with the findLB, findLBFwd and findLBRev functions in listing 7.

Listing 7: Searching For The Correct Pair
function findLB(x, state) {
 x = Number(x);
 if(x>state.nodes[state.lb+1].x) {
  if(state.lb<state.nodes.length-2) findLBFwd(x, state);
 }
 else if(x<state.nodes[state.lb].x) {
  if(state.lb>0) findLBRev(x, state);
 }
 return x;
}

function compare(l, r) {
 return l.x-r;
}

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

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

The findLB function restricts the search for the lower bound node in the same way as the jumpLB function did, deferring to the findLBFwd function if it's after the current lower bound and to findLBRev if it's before it.
These each first check whether the next node represents the new lower bound and, if not and it isn't in the pair at the end of the nodes, uses our ak._unsafeLowerBound function to search for it in the remaining positions. Since it returns the last point before which we can insert the value whilst keeping the array sorted we must search for the index of the upper bound of the pair and then subtract one from it.
Note that we are depending on the fact that we know that the first argument in the comparison function is always taken from the array being searched and the second argument is always the value for which we're searching for the lower bound index.

Finally, the copyNodes function given in listing 8 makes a deep copy of the nodes to be returned from the nodes property accessor that we add to the interpolator before returning it.

Listing 8: copyNodes
function copyNodes(nodes) {
 var n = nodes.length;
 var copy = new Array(n);
 var i;

 for(i=0;i<n;++i) copy[i] = {x:nodes[i].x, y:nodes[i].y};
 return copy;
}

You can find the implementation of ak.linearInterpolate in LinearInterpolate.js.

Using ak.linearInterpolate

Program 1 demonstrates how to use ak.linearInterpolate for real valued functions by using it to approximate the sine function, plotting the function in grey, the interpolation in black and marking the nodes with crosses.

Program 1: Interpolating A Real Valued Function

Program 2 does the same for complex valued functions using our ak.complex type to represent the values at the nodes and plotting the real and imaginary parts of the results in shades of red and green respectively.

Program 2: Interpolating A Complex Valued Function

The first thing to note about the results of both of these programs is that extrapolating beyond the first and last points quickly results in inaccurate values. Nevertheless, sometimes it is appropriate to do so and consequently I'm not prepared to treat extrapolation as an error but instead leave it up to the user to wrap the interpolator in a function that throws an exception when extrapolating should they want to.
To extrapolate with constant values, we can simply add a point to the left of the first with the same value of \(y\) as it has and another to the right of the last with the same value of \(y\) as its.

The next thing to note is that linear interpolation isn't smooth, which is to say that its derivative isn't continuous and, in particular, that it's undefined at the nodes. In some circumstances this isn't acceptable; the angles at the nodes may create points of material stress in engineering schematics, for example, or its graph might not be aesthetically pleasing and so in the next post we shall take a look at a smooth interpolation algorithm.

References

[1] IEEE standard for binary floating-point arithmetic. ANSI/IEEE std 754-1985, 1985.

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

[3] Climbing Mount Impractical, www.thusspakeak.com, 2013.

Leave a comment