Eigen Spline Derivative Example

Eigen Spline Derivative Example.md

The following code snippet showcases how one would use Eigen's spline module to fit x^2 and find the first and second derivatives along the spline.

An important note is to remember to scale the derivatives with:

• (x - x_min) / (x_max - x_min) (first derivative)
• (x - x_min) / (x_max - x_min)^2 (second derivative).

Because when fitting a spline to x^2, the x values you feed into the spline fitter is actually normalized with (x - x_min) / (x_max - x_min). So using the chain rule the first derivative of x^2 would be 2x / (x_max - x_min) and the second derivative would be 2 / (x_max - x_min)**2.

How to compile and run this example:

g++ -I/usr/include/eigen3 spline_deriv_demo.cpp -o spline_deriv_demo
./spline_deriv_demo

Example output:

1st deriv at x = 0: 4.52754e-16
1st deriv at x = 1: 2
1st deriv at x = 2: 4
1st deriv at x = 3: 6
1st deriv at x = 4: 8

2nd deriv at x = 0: 2
2nd deriv at x = 1: 2
2nd deriv at x = 2: 2
2nd deriv at x = 3: 2
2nd deriv at x = 4: 2

spline_deriv_demo.cpp

#include <iostream>
#include <Eigen/Core>
#include <unsupported/Eigen/Splines>

Eigen::VectorXd normalize(const Eigen::VectorXd &x) {
  const double min = x.minCoeff();
  const double max = x.maxCoeff();
  Eigen::VectorXd x_norm{x.size()};

  for (int k = 0; k < x.size(); k++) {
    x_norm(k) = (x(k) - min)/(max - min);
  }

  return x_norm;
}

int main() {
  // Create x-y data
  const Eigen::Vector4d x{0, 1, 2, 4};
  const Eigen::Vector4d y = (x.array().square());  // x^2

  // Fit and create spline
  typedef Eigen::Spline<double, 1, 3> Spline1D;
  typedef Eigen::SplineFitting<Spline1D> Spline1DFitting;
  const auto knots = normalize(x);
  Spline1D s = Spline1DFitting::Interpolate(y.transpose(), 3, knots);

  // First derivative
  printf("First derivative:\n");
  const double scale = 1 / (x.maxCoeff() -  x.minCoeff());
  printf("f'(x = 0): %.2f\n", s.derivatives(0.00, 1)(1) * scale);
  printf("f'(x = 1): %.2f\n", s.derivatives(0.25, 1)(1) * scale);
  printf("f'(x = 2): %.2f\n", s.derivatives(0.50, 1)(1) * scale);
  printf("f'(x = 3): %.2f\n", s.derivatives(0.75, 1)(1) * scale);
  printf("f'(x = 4): %.2f\n", s.derivatives(1.00, 1)(1) * scale);
  printf("\n");

  /**
   * IMPORTANT NOTE: Eigen's spline module is not documented well. Once you fit
   * a spline to find the derivative of the fitted spline at any point u [0, 1]
   * you call:
   *
   *     spline.derivatives(u, 1)(1)
   *                        ^  ^  ^
   *                        |  |  |
   *                        |  |  +------- Access the result
   *                        |  +---------- Derivative order
   *                        +------------- Parameter u [0, 1]
   *
   * The last bit `(1)` is if the spline is 1D. And value of `1` for the first
   * order. `2` for the second order. Do not forget to scale the result.
   *
   * For higher dimensions, treat the return as a matrix and grab the 1st or
   * 2nd column for the first and second derivative.
   */

  // Second derivative
  const double scale_sq = scale * scale;
  printf("Second derivative:\n");
  printf("f''(x = 0): %.2f\n", s.derivatives(0.00, 2)(2) * scale_sq);
  printf("f''(x = 1): %.2f\n", s.derivatives(0.25, 2)(2) * scale_sq);
  printf("f''(x = 2): %.2f\n", s.derivatives(0.50, 2)(2) * scale_sq);
  printf("f''(x = 3): %.2f\n", s.derivatives(0.75, 2)(2) * scale_sq);
  printf("f''(x = 4): %.2f\n", s.derivatives(1.00, 2)(2) * scale_sq);

  return 0;
}