We're Not For Turning


We have seen how it is possible to smoothly interpolate between a set of points \(\left(x_i, y_i\right)\), with the \(x_i\) known as nodes and the \(y_i\) as values, by specifying the gradients \(g_i\) at the nodes and calculating values between adjacent pairs using the uniquely defined cubic polynomials that match the values and gradients at them[1].
We have also seen how extrapolating such polynomials beyond the first and last nodes can yield less than satisfactory results, which we fixed by specifying the first and last gradients and then adding new first and last nodes to ensure that the first and last polynomials would represent straight lines.
Now we shall see how cubic spline interpolation can break down rather more dramatically and how we might fix it.

The Problem Of Monotonicity

Program 1 illustrates the problem with the set of points
\[ \left\{(0,10),(2,10),(3,10),(5,10),(6,10),(8,10),(9,10.5),(11,15),(12,50),(14,60),(15,85)\right\} \]
chosen by Akima[2].

Program 1: Interpolating Akima's Points

Here you can clearly see that even though the values of \(y_i\) are non-decreasing the interpolated values are not. Now this is not just an in issue of aesthetics since there are many kinds of functions that are, by definition, either non-decreasing or non-increasing, known as monotonic functions. For example, the cumulative distribution function \(F_X\) of a random variable \(X\) must be non-decreasing or else we would be able to find a pair of points \(x_0\) and \(x_1\) such that
\[ \begin{align*} x_0 &< x_1\\ F_X\left(x_0\right) &> F_X\left(x_1\right) \end{align*} \]
and consequently
\[ \Pr\left(x_0 < X \leqslant x_1\right) = F_X\left(x_1\right) - F_X\left(x_0\right) < 0 \]
which is impossible.

A simple way to ensure that such problems don't occur is to insist that each spline is monotonic so that the interpolation \(f\) will be non-decreasing between adjacent pairs of nodes with non-decreasing values and non-increasing between those with non-increasing values or, more specifically, that
\[ x_i \leqslant a < b \leqslant x_{i+1} \rightarrow \begin{cases} f(a) < f(b) &\quad\quad y_i < y_{i+1}\\ f(a) = f(b) &\quad\quad y_i = y_{i+1}\\ f(a) > f(b) &\quad\quad y_i > y_{i+1} \end{cases} \]
where \(\rightarrow\) means implies and the left curly brace is followed by expressions that are conditional upon those to their right.

Tweaking The Gradients

We can exert a lot of control over the shape of a given spline by adjusting the gradients at its endpoints, as demonstrated by program 2.

Program 2: The Effect Of The Endpoint Gradients

This moves from a non-monotonic curve, plotted in red, through a monotonic curve, plotted in green, and back to a non-monotonic curve as those gradients increase.

The obvious question is whether or not this gives us sufficient control to set the gradients so that every spline is monotonic and, thankfully, there's a very simple way to show that it does; set them all to zero.
Having done so, if the values at a pair of adjacent nodes are equal then the spline connecting them must be a horizontal line since the splines are uniquely defined by the values and gradients at the nodes and a horizontal line meets that definition. If they are not then we need simply note that a cubic polynomial cannot have a derivative of zero at more than two points unless it is a horizontal line and since we have set it to zero at both of the endpoints it can't be zero anywhere between them and so, by the intermediate value theorem, can't change sign.
In both cases the spline is monotonic and so if we set the gradients to zero we will have an interpolation that is monotonic between the nodes, referred to as being piecewise monotonic, although it probably won't be very natural looking, as demonstrated by program 3.

Program 3: Setting The Gradients To Zero

Fortunately for our aesthetic sensibilities we are not limited to setting the endpoint gradients to zero to ensure a cubic spline's monotonicity, as was demonstrated by program 2. Specifically, given left and right endpoints of \(\left(x_0,y_0\right)\) and \(\left(x_1,y_1\right)\), if we define
\[ \begin{align*} \Delta x &= x_1 - x_0 & \Delta y &= y_1 - y_0\\ g_0 &= s_0 \times \tfrac{\Delta y}{\Delta x} & g_1 &= s_1 \times \tfrac{\Delta y}{\Delta x} \end{align*} \]
then it will be monotonic if both \(s_0\) and \(s_1\) are no less than zero and no greater than three, as shown by Fritsch and Carlson[3].

Proving Monotonicity

If \(\Delta y\) equals zero then both \(g_0\) and \(g_1\) will equal zero, as we have already seen is required, for any finite \(s_0\) and \(s_1\), so we might as well just set them both to one in such cases.
Otherwise, that they cannot be negative is trivially true since if they were the sign of their associated gradients would be the opposite of that of \(\Delta y\) and the curve would have to change from increasing/decreasing to decreasing/increasing at some point to get from the one endpoint to the other.
To prove that they are bounded above by three we have three shapes of the spline to consider; linear, quadratic and cubic curves.

If the curve is linear then it is trivially monotonic and the endpoint gradients must equal the finite difference \(\frac{\Delta y}{\Delta x}\) and so both \(s_0\) and \(s_1\) will equal one.

Next, if it's quadratic then it must take the form
\[ y = at^2 - \frac{s_0}{s_0-1}at + y_0 \]
for some non-zero value of \(a\), where
\[ t = \frac{x-x_0}{\Delta x}\\ \]
so that
\[ \begin{align*} t_0 &= \frac{x_0-x_0}{\Delta x} = \frac{0}{\Delta x} = 0\\ t_1 &= \frac{x_1-x_0}{\Delta x} = \frac{\Delta x}{\Delta x} = 1 \end{align*} \]
and it must be the case that
\[ s_0+s_1=2 \]
as proven in derivation 1.

Derivation 1: The Quadratic Case
Firstly, we have
\[ \begin{align*} y &= at^2 + bt + c& y_0 &= at_0^2 + bt_0 + c = c& y_1 &= at_1^2 + bt_1 + c = a+b+c\\ g &= \frac{\mathrm{d}y}{\mathrm{d}t} \times \frac{\mathrm{d}t}{\mathrm{d}x} = \frac{\mathrm{d}y}{\mathrm{d}x} = \frac{2at+b}{\Delta x}& g_0 &= \frac{2at_0+b}{\Delta x} = \frac{b}{\Delta x}& g_1 &= \frac{2at_1+b}{\Delta x} = \frac{2a+b}{\Delta x} \end{align*} \]
and therefore
\[ \Delta y = y_1 - y_0 = (a+b+c)-c = a+b \]
Recalling that
\[ \begin{align*} g_0 &= s_0 \times \frac{\Delta y}{\Delta x}& g_1 &= s_1 \times \frac{\Delta y}{\Delta x} \end{align*} \]
we also have the identities
\[ \begin{align*} g_0 \Delta x &= b = s_0 \Delta y = s_0 (a+b)& g_1 \Delta x &= 2a+b = s_1 \Delta y = s_1 (a+b) \end{align*} \]
Rearranging yields
\[ \begin{align*} s_0a + \left(s_0-1\right)b &= 0 & b = -\frac{s_0}{s_0-1}a\\ \left(s_1-2\right)a + \left(s_1-1\right)b &= 0 & b = -\frac{s_1-2}{s_1-1}a \end{align*} \]
so that
\[ \begin{align*} y = at^2 - \frac{s_0}{s_0-1}at + y_0 \end{align*} \]
\[ \begin{align*} \frac{s_0}{s_0-1} &= \frac{s_1-2}{s_1-1}\\ s_0\left(s_1-1\right) &= \left(s_0-1\right)\left(s_1-2\right)\\ s_0s_1-s_0 &= s_0s_1-2s_0-s_1+2\\ s_0 + s_1 &= 2 \end{align*} \]

The turning point \(x^\ast\) at which a quadratic curve
\[ ax^2+bx+c \]
changes from increasing to decreasing, or vice versa, occurs when its derivative equals zero
\[ 2ax^\ast + b = 0 \]
so that
\[ x^\ast = -\frac{b}{2a} \]
which in this case means that
\[ t^\ast = \frac{s_0}{2s_0-2} \]
This can only lie between the nodes of the spline if
\[ 0 < \frac{s_0}{2s_0-2} < 1 \]
Now we know that \(s_0\) cannot be negative and if it is less than one then the denominator will be negative and so
\[ \frac{s_0}{2s_0-2} \leqslant 0 \]
If it is greater than one then the denominator is positive and we have
\[ \begin{align*} \frac{s_0}{2s_0-2} &> 0\\ \frac{s_0}{2s_0-2} &< 1 \rightarrow s_0 < 2s_0-2 \rightarrow s_0 > 2 \end{align*} \]
which have already shown cannot be the case since neither \(s_0\) nor \(s_1\) can be negative and their sum must equal two.
If it equals one then the turning point will occur at infinity and so we may conclude that, if it's quadratic, the curve will be monotonic for all valid values of \(s_0\) and \(s_1\).

Finally, if the curve is cubic then it must take the form
\[ y = at^3 - \frac{2s_0+s_1-3}{s_0+s_1-2}at^2 +\frac{s_0}{s_0+s_1-2}at + y_0 \]
as proven in derivation 2

Derivation 2: The Cubic Case
This time we have
\[ \begin{align*} y &= at^3 + bt^2 + ct + d& y_0 &= d& y_1 &= a+b+c+d\\ g &= \frac{3at^2+2bt+c}{\Delta x}& g_0 &= \frac{c}{\Delta x}& g_1 &= \frac{3a+2b+c}{\Delta x} \end{align*} \]
and so
\[ \Delta y = y_1-y_0 = (a+b+c+d) - d = a+b+c \]
Multiplying the endpoint gradients by \(\Delta x\) yields
\[ \begin{align*} g_0 \Delta x &= c = s_0 \Delta y = s_0 (a+b+c)& g_1 \Delta x &= 3a+2b+c = s_1 \Delta y = s_1 (a+b+c) \end{align*} \]
and consequently
\[ \begin{align*} s_0a + s_0b + \left(s_0-1\right)c &= 0& c &= -\frac{s_0a + s_0b}{s_0-1}\\ \left(s_1-3\right)a + \left(s_1-2\right)b + \left(s_1-1\right)c &=0& c &= -\frac{\left(s_1-3\right)a + \left(s_1-2\right)b}{s_1-1}\\ \end{align*} \]
Solving for \(c\) gives us
\[ \begin{align*} \frac{s_0a + s_0b}{s_0-1} &= \frac{\left(s_1-3\right)a + \left(s_1-2\right)b}{s_1-1}\\ \left(\frac{s_0}{s_0-1} - \frac{s_1-2}{s_1-1}\right)b &= \left(\frac{s_1-3}{s_1-1} - \frac{s_0}{s_0-1}\right)a\\ \left(s_0\left(s_1-1\right) - \left(s_0-1\right)\left(s_1-2\right)\right)b &= \left(\left(s_0-1\right)\left(s_1-3\right) - s_0\left(s_1-1\right)\right)a\\ \left(\left(s_0s_1-s_0\right) - \left(s_0s_1-2s_0-s_1+2\right)\right)b &= \left(\left(s_0s_1-3s_0-s_1+3\right) - \left(s_0s_1-s_0\right)\right)a\\ \left(s_0+s_1-2\right)b &= \left(-2s_0-s_1+3\right)a\\ b &= -\frac{2s_0+s_1-3}{s_0+s_1-2}a \end{align*} \]
and so the spline must take the form
\[ \begin{align*} y &= at^3 - \frac{2s_0+s_1-3}{s_0+s_1-2}at^2 -\frac{s_0a - \frac{2s_0+s_1-3}{s_0+s_1-2}s_0a}{s_0-1}t + y_0\\ &= at^3 - \frac{2s_0+s_1-3}{s_0+s_1-2}at^2 -\frac{\left(\left(s_0+s_1-2\right) - \left(2s_0+s_1-3\right)\right)s_0}{\left(s_0-1\right)\left(s_0+s_1-2\right)}at + y_0\\ &= at^3 - \frac{2s_0+s_1-3}{s_0+s_1-2}at^2 +\frac{\left(s_0-1\right)s_0}{\left(s_0-1\right)\left(s_0+s_1-2\right)}at + y_0\\ &= at^3 - \frac{2s_0+s_1-3}{s_0+s_1-2}at^2 +\frac{s_0}{s_0+s_1-2}at + y_0 \end{align*} \]
Note that it cannot be the case that
\[ s_0+s_1=2 \]
since this would lead to infinite coefficients and thus contradict our assumption that the curve is cubic.

and its turning points will be
\[ \begin{align*} t^- &= \frac{\left(2s_0+s_1-3\right) - \sqrt{\left(s_0+s_1-3\right)^2-s_0s_1}}{3s_0+3s_1-6}\\ t^+ &= \frac{\left(2s_0+s_1-3\right) + \sqrt{\left(s_0+s_1-3\right)^2-s_0s_1}}{3s_0+3s_1-6} \end{align*} \]
as proven in derivation 3.

Derivation 3: The Cubic Turning Points
The derivative of the curve is given by
\[ \frac{\mathrm{d}y}{\mathrm{d}t} = 3at^2 - \frac{4s_0+2s_1-6}{s_0+s_1-2}at +\frac{s_0}{s_0+s_1-2}a \]
and its turning points must therefore satisfy
\[ \begin{align*} 3t^{\ast2} - \frac{4s_0+2s_1-6}{s_0+s_1-2}t^\ast +\frac{s_0}{s_0+s_1-2} &= 0\\ \left(3s_0+3s_1-6\right)t^{\ast2} - \left(4s_0+2s_1-6\right)t^\ast + s_0 &= 0 \end{align*} \]
By the formula for solving quadratic equations we have
\[ \begin{align*} t^\ast &= \frac{\left(4s_0+2s_1-6\right) \pm \sqrt{\left(4s_0+2s_1-6\right)^2 - 4\left(3s_0+3s_1-6\right)s_0}}{2\left(3s_0+3s_1-6\right)}\\ &= \frac{\left(2s_0+s_1-3\right) \pm \sqrt{\left(2s_0+s_1-3\right)^2 - \left(3s_0+3s_1-6\right)s_0}}{3s_0+3s_1-6}\\ &= \frac{\left(2s_0+s_1-3\right) \pm \sqrt{\left(4s_0^2+s_1^2+9+4s_0s_1-12s_0-6s_1\right) - \left(3s_0^2+3s_0s_1-6s_0\right)}}{3s_0+3s_1-6}\\ &= \frac{\left(2s_0+s_1-3\right) \pm \sqrt{s_0^2+s_1^2+9+s_0s_1-6s_0-6s_1}}{3s_0+3s_1-6}\\ &= \frac{\left(2s_0+s_1-3\right) \pm \sqrt{\left(s_0+s_1-3\right)^2-s_0s_1}}{3s_0+3s_1-6} \end{align*} \]

The turning points at the extremes of \(s_0\) and \(s_1\) are
\[ \begin{matrix} s_0 & s_1 & t^+ & t^-\\ 0 & 0 & 0 & 1\\ 0 & 3 & 0 & 0\\ 3 & 0 & 1 & 1\\ 3 & 3 & \frac12 & \frac12\\ \end{matrix} \]
Note that since we know that the curve is cubic, the last case must be a point of inflection where it goes from increasing to flat to instantaneously back to increasing, or decreasing to flat to decreasing. In the other cases the turning points occur at the endpoints of the curve and so it must be monotonic between them in all four cases.
Derivation 4 goes further by showing that the curve must be monotonic if either \(s_0\) or \(s_1\) are equal to zero or three.

Derivation 4: Monotonicity At The Edges
Firstly, if \(s_0\) equals zero then
\[ t^\ast = \frac{\left(s_1-3\right) \pm \sqrt{\left(s_1-3\right)^2}}{3s_1-6} = \frac{\left(s_1-3\right) \pm \left(s_1-3\right)}{3s_1-6} \]
and so
\[ \begin{align*} t^- &= 0& t^+ &= \frac{2s_1-6}{3s_1-6} = \frac{6-2s_1}{6-3s_1} \end{align*} \]
Since \(t^-\) is always zero we require that \(t^+\) cannot be greater than zero or less than one if \(s_1\) lies between zero and three which we can verify with the cases
\[ \begin{align*} s_0 = 0 &\rightarrow t^+ = 1\\ s_0 < 2 &\rightarrow \left(6-2s_1>0 \wedge 6-3s_1>0 \wedge 6-2s_1 > 6-3s_1\right) \rightarrow t^+ > 1\\ s_0 = 2 &\rightarrow t^+ = -\infty\\ s_0 > 2 &\rightarrow \left(6-2s_1>0 \wedge 6-3s_1<0\right) \rightarrow t^+ < 0\\ s_0 = 3 &\rightarrow t^+ = 0 \end{align*} \]
where \(\wedge\) means and.

Secondly, if \(s_0\) equals three we have
\[ t^\ast = \frac{\left(s_1+3\right) \pm \sqrt{s_1^2-3s_1}}{3s_1+3} \]
and the term inside the square root will be negative if \(s_1\) is greater than zero or less than three so that the curve won't have any turning points. If it equals zero or three then that term will equal zero and so \(t^-\) will equal \(t^+\) and will be a point of inflection.

Thirdly, if \(s_1\) equals zero then
\[ t^\ast = \frac{\left(2s_0-3\right) \pm \sqrt{\left(s_0-3\right)^2}}{3s_0-6} = \frac{\left(2s_0-3\right) \pm \left(s_0-3\right)}{3s_0-6} \]
\[ \begin{align*} t^- &= \frac{s_0}{3s_0-6} & t^+ &= 1 \end{align*} \]
This time \(t^+\) always equals one and we can confirm that \(t^-\) lies between zero and one with
\[ \begin{align*} s_0 = 0 &\rightarrow t^- = 0\\ s_0 < 2 &\rightarrow t^- <0\\ s_0 = 2 &\rightarrow t^- = +\infty\\ s_0 > 2 &\rightarrow t^- > 1\\ s_0 = 3 &\rightarrow t^- = 1 \end{align*} \]
Finally, \(s_1\) equal to three yields
\[ \begin{align*} t^\ast &= \frac{2s_0 \pm \sqrt{s_0^2-3s_0}}{3s_0+3} \end{align*} \]
and monotonicity follows by the same argument as made for \(s_0\) equal to three.

The easiest way to demonstrate that the conditions for monotonicity are met between the specified bounds is to draw a diagram showing where they are and are not. Program 4 does so by plotting turning points below zero in red, those above one in green and those that don't exist in yellow for \(s_0\) and \(s_1\) between \(-1\frac14\) and \(4\frac14\), with \(t^-\) at the top and \(t^+\) beneath, marking the bounds with black lines.

Program 4: The Graphical Solution

Whilst this clearly illustrates that the bounds are sufficient to ensure monotonicity, it also proves that they are not necessary since the yellow ellipses extend beyond them equally for both \(t^-\) and \(t^+\). For example, if \(s_0\) equals one and \(s_1\) equals three and a half, we have
\[ \left(s_0+s_1-3\right)^2-s_0s_1 = \left(1+3\tfrac12-3\right)^2-1 \times 3\tfrac12 = \left(1\tfrac12\right)^2-3\tfrac12 = 2\tfrac14 - 3\tfrac12 = -1\tfrac14 \]
so that the square roots in the formulae for \(t^-\) and \(t^+\) are imaginary and consequently the curve doesn't have any turning points.

The Piecewise Monotonic Cubic Spline Implementation

A significant benefit of the stricter bounds is that for any \(s_0\) within them every \(s_1\) within them must yield a monotonic curve and vice-versa.
We can therefore iterate over the nodes and change their gradients one at a time whenever they would lead to a non-monotonic spline; if their gradients would need to have different signs to meet the requirements for the splines to their left and to their right then we just set them to zero and, if not, the bound with the smaller magnitude must fall within the bound with the larger and so we simply choose the smaller.
The toMono function given in listing 1 updates the gradients of the nodes of an ak.cubicSplineInterpolate interpolator f so that their ratios to the finite differences between them lie between zero and ub, ensuring that they fall within the bounds required by the splines to their left and right by iterating over the splines and modifying the gradients at their endpoints so that each will meet the conditions of both.

Listing 1: Updating The Gradients
function toMono(ub, f) {
 var nodes = f.nodes();
 var changed = false;
 var n = nodes.length;
 var x0 = nodes[0].x;
 var y0 = nodes[0].y;
 var g0 = nodes[0].g;
 var i, x1, y1, g1, g, ubg;

 for(i=1;i<n;++i) {
  x1 = nodes[i].x;
  y1 = nodes[i].y;
  g1 = nodes[i].g;
  g = (y1-y0)/(x1-x0);
  ubg = ub*g;

  if(g0!==0) {
   if(g0*g<=0) {changed=true; nodes[i-1].g = 0;}
   else if(Math.abs(g0)>Math.abs(ubg)) {changed=true; nodes[i-1].g = ubg;}

  if(g1!==0) {
   if(g1*g<=0) {changed=true; nodes[i].g = g1 = 0;}
   else if(Math.abs(g1)>Math.abs(ubg)) {changed=true; nodes[i].g = g1 = ubg;}

  x0 = x1;
  y0 = y1;
  g0 = g1;

 return changed ? ak.cubicSplineInterpolate(nodes) : f;

Note that we're keeping track of whether any of the gradients were changed and only create a new ak.cubicSplineInterpolate interpolator if they were, returning the original one if not.

Defined in MonoCubicSplineInterpolate.js, the ak.monoCubicSplineInterpolate function given in listing 2 creates a monotonic cubic spline interpolator with an optional first argument, defaulting to \(0.89\) and which we determine has been specified by checking whether or not it's a number, defining how well it should preserve the gradients that an ak.cubicSplineInterpolate interpolator has after it's been constructed with the remaining arguments. Specifically, the upper bound of the ratio of the gradients to the finite differences passed to toMono as ub will be zero if it is zero and three if it's equal to one.

Listing 2: ak.monoCubicSplineInterpolate
ak.monoCubicSplineInterpolate = function() {
 var args = arguments;
 var arg0 = args[0];
 var p = 0.89;
 var f;

 if(ak.nativeType(arg0)===ak.NUMBER_T) {
  p = Number(arg0);
  if(!(p>=0 && p<=1)) {
   throw new Error('invalid gradient preservation in ak.monoCubicSplineInterpolate');
  args = Array.prototype.slice.call(args, 1);

 try {f = ak.cubicSplineInterpolate.apply(null, args);}
 catch(e) {throw new Error(e.message.replace('ak.cubic', 'ak.monoCubic'));}
 return toMono(3*p, f);

Program 5 shows how this removes the turning points from the interpolation whilst preserving a relatively natural looking shape that setting the gradients to zero did not.

Program 5: Monotonically Interpolating Akima's Points

Whilst this certainly yields an aesthetically pleasing curve, there's no mathematical justification for any particular upper bound on the ratio of the gradients and the finite differences and the choice has a visible effect upon its shape, as demonstrated by program 6.

Program 6: Changing The Upper Bound

Fritsch and Carlson recommended constraining the ratios at the endpoints of a spline to a quarter circle of radius three so that
\[ \begin{align*} s_0 &\geqslant 0\\ s_1 &\geqslant 0\\ s_0^2+s_1^2 &\leqslant 9 \end{align*} \]
Unfortunately this introduces a dependency between them which means that the final gradients of the cubic splines will depend upon the order in which we adjust them to meet this condition. That is unless we add a further condition, like choosing the largest ratios that do, which will make them more difficult to calculate.
It was because of this that I chose to constrain them to a square that is some fraction of the largest in which the splines are monotonic, although I decided on the default of \(0.89\) because it resulted in a square with an area approximately equal to the quarter circle that they had recommended.

This is indicative of the problem with cubic spline interpolation that I mentioned last time as the reason why I chose not to generalise our implementation to other arithmetic types, like vectors and matrices; the final shape of the interpolation depends upon choices that have nothing whatsoever to do with the points that we are interpolating between. This is in stark contract to linear interpolation which represents, in a mathematically justifiable sense, the best guess at the function that passes through the points.


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

[2] Akima, H., A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures, Journal of the ACM, Vol. 17, No. 4, 1970.

[3] Fritsch, F. & Carlson, R., Monotone Piecewise Cubic Interpolation, SIAM Journal on Numerical Analysis, Vol. 17, No. 2, 1980.

Leave a comment