Cubic B-spline interpolation

Synopsis

#include <boost/math/interpolators/cubic_b_spline.hpp>

namespace boost { namespace math {

  template <class Real>
  class cubic_b_spline
  {
  public:

      template <class BidiIterator>
      cubic_b_spline(BidiIterator a, BidiIterator b, Real left_endpoint, Real step_size,
                     Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
                     Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
      cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
                     Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
                     Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());

      Real operator()(Real x) const;

      Real prime(Real x) const;
  };

}} // namespaces

The cubic B-spline class provided by boost allows fast and accurate interpolation of a function which is known at equally spaced points. The cubic B-spline interpolation is numerically stable as it uses compactly supported basis functions constructed via iterative convolution. This is to be contrasted to traditional cubic spline interpolation is ill-conditioned as the global support of cubic polynomials causes small changes far from the evaluation point exert a large influence on the calculated value.

There are many use cases for interpolating a function at equally spaced points. One particularly important example is solving ODE's whose coefficients depend on data determined from experiment or numerical simulation. Since most ODE steppers are adaptive, they must be able to sample the coefficients at arbitrary points; not just at the points we know the values of our function.

The first two arguments to the constructor are either:

A pair of bidirectional iterators into the data, or
A pointer to the data, and a length of the data array.

These are then followed by:

The start of the functions domain
The step size

Optionally, you may provide two additional arguments to the constructor, namely the derivative of the function at the left endpoint, and the derivative at the right endpoint. If you do not provide these arguments, they will be estimated using one-sided finite-difference formulas. An example of a valid call to the constructor is

std::vector<double> f{0.01, -0.02, 0.3, 0.8, 1.9, -8.78, -22.6};
double t0 = 0;
double h = 0.01;
boost::math::cubic_b_spline<double> spline(f.begin(), f.end(), t0, h);

The endpoints are estimated using a one-sided finite-difference formula. If you know the derivative at the endpoint, you may pass it to the constructor via

boost::math::cubic_b_spline<double> spline(f.begin(), f.end(), t0, h, a_prime, b_prime);

To evaluate the interpolant at a point, we simply use

double y = spline(x);

and to evaluate the derivative of the interpolant we use

double yp = spline.prime(x);

Be aware that the accuracy guarantees on the derivative of the spline are an order lower than the guarantees on the original function, see Numerical Analysis, Graduate Texts in Mathematics, 181, Rainer Kress for details.

Finally, note that this is an interpolator, not an extrapolator. Therefore, you should strenuously avoid evaluating the spline outside the endpoints. However, it is not an error if you do, as often you cannot control where (say) an ODE stepper will evaluate your function. As such the interpolant tries to do something reasonable when the passed a value outside the endpoints. For evaluation within one stepsize of the interval, you can assume something somewhat reasonable was returned. As you move further away from the endpoints, the interpolant decays to its average on the interval.

Complexity and Performance

The call to the constructor requires 𝑶(n) operations, where n is the number of points to interpolate. Each call the the interpolant is 𝑶(1) (constant time). On the author's Intel Xeon E3-1230, this takes 21ns as long as the vector is small enough to fit in cache.

#include <boost/math/interpolators/cubic_b_spline.hpp>

int main()
{
    // We begin with an array of samples:
    std::vector<double> v(500);
    // And decide on a stepsize:
    double step = 0.01;

    // Initialize the vector with a function we'd like to interpolate:
    for (size_t i = 0; i < v.size(); ++i)
    {
        v[i] = sin(i*step);
    }
    // We could define an arbitrary start time, but for now we'll just use 0:
    boost::math::cubic_b_spline<double> spline(v.data(), v.size(), 0 /* start time */, step);

    // Now we can evaluate the spline wherever we please.
    std::mt19937 gen;
    boost::random::uniform_real_distribution<double> absissa(0, v.size()*step);
    for (size_t i = 0; i < 10; ++i)
    {
        double x = absissa(gen);
        std::cout << "sin(" << x << ") = " << sin(x) << ", spline interpolation gives " << spline(x) << std::endl;
        std::cout << "cos(" << x << ") = " << cos(x) << ", spline derivative interpolation gives " << spline.prime(x) << std::endl;
    }

    // The next example is less trivial:
    // We will try to figure out when the population of the United States crossed 100 million.
    // Since the census is taken every 10 years, the data is equally spaced, so we can use the cubic b spline.
    // Data taken from https://en.wikipedia.org/wiki/United_States_Census
    // An eye
    // We'll start at the year 1860:
    double t0 = 1860;
    double time_step = 10;
    std::vector<double> population{31443321,  /* 1860 */
                                   39818449,  /* 1870 */
                                   50189209,  /* 1880 */
                                   62947714,  /* 1890 */
                                   76212168,  /* 1900 */
                                   92228496,  /* 1910 */
                                   106021537, /* 1920 */
                                   122775046, /* 1930 */
                                   132164569, /* 1940 */
                                   150697361, /* 1950 */
                                   179323175};/* 1960 */

    // An eyeball estimate indicates that the population crossed 100 million around 1915.
    // Let's see what interpolation says:
    boost::math::cubic_b_spline<double> p(population.data(), population.size(), t0, 10);

    // Now create a function which has a zero at p = 100,000,000:
    auto f = [=](double t){ return p(t) - 100000000; };

    // Boost includes a bisection algorithm, which is robust, though not as fast as some others 
    // we provide, but lets try that first.  We need a termination condition for it, which
    // takes the two endpoints of the range and returns either true (stop) or false (keep going),
    // we could use a predefined one such as boost::math::tools::eps_tolerance<double>, but that
    // won't stop until we have full double precision which is overkill, since we just need the 
    // endpoint to yield the same month.  While we're at it, we'll keep track of the number of
    // iterations required too, though this is strictly optional:

    auto termination = [](double left, double right)
    {
       double left_month = std::round((left - std::floor(left)) * 12 + 1);
       double right_month = std::round((right - std::floor(right)) * 12 + 1);
       return (left_month == right_month) && (std::floor(left) == std::floor(right));
    };
    std::uintmax_t iterations = 1000;
    auto result =  boost::math::tools::bisect(f, 1910.0, 1920.0, termination, iterations);
    auto time = result.first;  // termination condition ensures that both endpoints yield the same result
    auto month = std::round((time - std::floor(time))*12  + 1);
    auto year = std::floor(time);
    std::cout << "The population of the United States surpassed 100 million on the ";
    std::cout << month << "th month of " << year << std::endl;
    std::cout << "Found in " << iterations << " iterations" << std::endl;

    // Since the cubic B spline offers the first derivative, we could equally have used Newton iterations,
    // this takes "number of bits correct" as a termination condition - 20 should be plenty for what we need,
    // and once again, we track how many iterations are taken:

    auto f_n = [=](double t) { return std::make_pair(p(t) - 100000000, p.prime(t)); };
    iterations = 1000;
    time = boost::math::tools::newton_raphson_iterate(f_n, 1910.0, 1900.0, 2000.0, 20, iterations);
    month = std::round((time - std::floor(time))*12  + 1);
    year = std::floor(time);
    std::cout << "The population of the United States surpassed 100 million on the ";
    std::cout << month << "th month of " << year << std::endl;
    std::cout << "Found in " << iterations << " iterations" << std::endl;

}

Program output is:

sin(4.07362) = -0.802829, spline interpolation gives - 0.802829
cos(4.07362) = -0.596209, spline derivative interpolation gives - 0.596209
sin(0.677385) = 0.626758, spline interpolation gives 0.626758
cos(0.677385) = 0.779214, spline derivative interpolation gives 0.779214
sin(4.52896) = -0.983224, spline interpolation gives - 0.983224
cos(4.52896) = -0.182402, spline derivative interpolation gives - 0.182402
sin(4.17504) = -0.85907, spline interpolation gives - 0.85907
cos(4.17504) = -0.511858, spline derivative interpolation gives - 0.511858
sin(0.634934) = 0.593124, spline interpolation gives 0.593124
cos(0.634934) = 0.805111, spline derivative interpolation gives 0.805111
sin(4.84434) = -0.991307, spline interpolation gives - 0.991307
cos(4.84434) = 0.131567, spline derivative interpolation gives 0.131567
sin(4.56688) = -0.989432, spline interpolation gives - 0.989432
cos(4.56688) = -0.144997, spline derivative interpolation gives - 0.144997
sin(1.10517) = 0.893541, spline interpolation gives 0.893541
cos(1.10517) = 0.448982, spline derivative interpolation gives 0.448982
sin(3.1618) = -0.0202022, spline interpolation gives - 0.0202022
cos(3.1618) = -0.999796, spline derivative interpolation gives - 0.999796
sin(1.54084) = 0.999551, spline interpolation gives 0.999551
cos(1.54084) = 0.0299566, spline derivative interpolation gives 0.0299566
The population of the United States surpassed 100 million on the 11th month of 1915
Found in 12 iterations
The population of the United States surpassed 100 million on the 11th month of 1915
Found in 3 iterations