Fritsch-Butland interpolation fails in a simple case

[tomwhoiscontrary]

I am trying to do some cubic interpolation. In some cases, the function i am interpolating has a constant value across its domain except at one point. Like this:

x	y
0	1
1	2
2	1
3	1
4	1

I would expect any monotone interpolation of these values to produce a value that rises between 0 and 1, falls between 1 and 2, and is constant between 2 and 4. Indeed, this is what i see using Kruger, harmonic, monotone spline, or monotone Akima derivative approximations. However, using Fritsch-Butland approximation, the value rises between 0 and 1, falls between 1 and 2, and is NaN between 2 and 4!

Should Fritsch-Butland approximation work in this situation, or does it have limitations i'm not aware of?

Here's the test code i'm using:

#include <iostream>
#include <vector>

#include <ql/math/interpolations/cubicinterpolation.hpp>

int main() {
        std::vector<double> knot_xs = {0, 1, 2, 3, 4};
        std::vector<double> knot_ys = {1, 2, 1, 1, 1};

        QuantLib::FritschButlandCubic interpolation(knot_xs.begin(),
                                                    knot_xs.end(),
                                                    knot_ys.begin());

        std::cout << "x,knot,interpolation\n";

        int resolution = 10;

        std::vector<double>::iterator knot_xs_iterator = knot_xs.begin();
        std::vector<double>::iterator knot_ys_iterator = knot_ys.begin();
        for (int i = 0; i <= knot_xs.back() * resolution; i++) {
                double x = (double)i / resolution;

                std::cout << x << ",";

                if (x >= *knot_xs_iterator) {
                        std::cout << *knot_ys_iterator;
                        knot_xs_iterator++;
                        knot_ys_iterator++;
                }

                std::cout << "," << interpolation(x) << '\n';
        }

        return 0;
}

This produces this output:

x,knot,interpolation
0,1,1
0.1,,1.163
0.2,,1.264
0.3,,1.321
0.4,,1.352
0.5,,1.375
0.6,,1.408
0.7,,1.469
0.8,,1.576
0.9,,1.747
1,2,2
1.1,,2.215
1.2,,2.28
1.3,,2.225
1.4,,2.08
1.5,,1.875
1.6,,1.64
1.7,,1.405
1.8,,1.2
1.9,,1.055
2,1,-nan
2.1,,-nan
2.2,,-nan
2.3,,-nan
2.4,,-nan
2.5,,-nan
2.6,,-nan
2.7,,-nan
2.8,,-nan
2.9,,-nan
3,1,-nan
3.1,,-nan
3.2,,-nan
3.3,,-nan
3.4,,-nan
3.5,,-nan
3.6,,-nan
3.7,,-nan
3.8,,-nan
3.9,,-nan
4,1,-nan

[tomwhoiscontrary]

The (non-monotonic) Akima approximation also produces NaNs in this case, although only half as many. Monotonic Akima is fine.

[lballabio]

Definitely not expected. Thanks for the heads-up.

If you (or anybody else) have time to debug the issue, I'd be grateful.

[tomwhoiscontrary]

I'll try to make time to dig into this - it should be simple to track down. However, like everyone, i have lots of other priorities, so i can't make any promises about when i'll get to it!

[lballabio]

Of course. Thanks!

[Grant6899]

Hi, I took a quick look last night.

I’m not an expert on interpolation but first should we make sure that it’s assumed f(x) is monotone?

[Grant6899]

Seems like Fritsch-Butland interpolation should have at least non-strict monotone assumption. See MONOTONE PIECEWISE CUBIC INTERPOLATION.

So we should make data below work:

x	y
0	10
2	10
3	10
5	10
6	10
8	10
9	10.5
11	15
12	50
14	60
15	85

And the root cause is when calling update() in impl, tmp_ array's calculation involves DIV/0 issue.

So questions:

1. In quantlib, is Fritsch-Butland designed only for monotone data?

2. If so, it should only be strict monotone or equality is acceptable in the sequence?

[tomwhoiscontrary]

@Grant6899 raises a point here: what does "monotonic" mean? Hyman's 1983 paper Accurate Monotonicity Preserving Cubic Interpolation talks about local monotonicity - if there is some range of the input data which is monotonic, then it should also be monotonic in the interpolation. In other words, Hyman's algorithm works on my original data. Or am i confused?

So, i had assumed that this is what we mean by "monotonic" for all the interpolation methods. Is that not the case? Does Fritsch-Butland require the input data to be globally monotonic?

I note that the comment on the CubicInterpolation class mentions:

   While some non-linear schemes (Modified Parabolic, Fritsch-Butland,
   Kruger) are guaranteed to be locally monotonic in their original
   approximation, all other schemes must be filtered according to the
   Hyman criteria at the expense of their linearity.

Which certainly suggests to me that Fritsch-Butland is locally, rather than globally, monotonic.

If so, i would suggest that the first change is to clarify the documentation, and to add asserts checking that the input data is indeed suitable.

[tomwhoiscontrary]

I note that if you jitter the y values in my example very slightly:

        std::vector<double> knot_ys = {1, 2, 1.000002, 1.000001, 1.000003};

Then Fritsch-Butland interpolates this quite happily. Well, i say happily; the values are not NaN, but are very silly indeed.

I think this is a classic numerical stability issue. I bet that somewhere in the algorithm, we are taking the difference between successive y values, and then dividing something by it. When they are the same, that makes a NaN; when they are very close, that makes some oversized values. If we can find where that is happening, we can see if we can reformulate that step in more stable terms.

We can track down the troublesome step by breaking the algorithm into a sequence of elementary arithmetic operations, and checking each one for NaNity. This is something so simple even i can do it - but it's tedious enough that i haven't found time to do it yet!

[lballabio]

I also had assumed that the interpolation preserves local monotonicity, but doesn't require it globally.