Polyfit function in C++ using ordinary least squares fit

polyfit.hpp

#include <vector>
#include "Eigen/Dense"

template <int... Orders>
Eigen::ArrayXd polyfit1D(const Eigen::ArrayXd& x, const Eigen::ArrayXd& y) {
  std::vector<int> orders{{Orders...}};
  if (x.size() != y.size()) {
    throw std::runtime_error("Polyfit failed: x and y must have same size!");
  }

  if (x.size() < orders.size()) {
    throw std::runtime_error(
        "Polyfit failed: Number of observations must be at least as many as "
        "the requested number of orders");
  }

  // Setup system of equations
  Eigen::MatrixXd A(x.size(), orders.size());
  for (size_t i = 0; i < x.size(); ++i) {
    for (size_t j = 0; j < orders.size(); ++j) {
      A(i, j) = std::pow(x[i], orders[j]);
    }
  }

  // Solve
  return A.colPivHouseholderQr().solve(y.matrix());
}

// Use like:
//
// auto coeffs = polyfit1D<0,2,4>(x_values, y_values);
//
// Will try to fit a polynomial of orders 0, 2 and 4 to the datapoints given.