Eigen で N 次の多項式の係数を求める

#include <Eigen/Dense>
#include <cmath>
#include <iostream>

bool PolyVal(const Eigen::VectorXd& x, const Eigen::VectorXd& y, int order,
             Eigen::VectorXd& p, double* r2) {
  using namespace Eigen;

  if (order < 1)
    return false;

  int N = order + 1;

  MatrixXd W(y.rows(), N);
  VectorXd w(N);
  for (int row = 0; row < W.rows(); row++) {
    for (int i = 0; i < N; i++)
      w(i) = std::pow(x(row), order - i);
    W.row(row) = w;
  }

  p = W.jacobiSvd(ComputeThinU | ComputeThinV).solve(y);

  VectorXd pred(y.rows());
  pred = W * p;

  VectorXd mean = VectorXd::Constant(y.rows(), y.mean());
  double ss_res = ((y - pred).transpose() * (y - pred)).sum();
  double ss_tot = ((y - mean).transpose() * (y - mean)).sum();

  *r2 = 1.0 - ss_res / ss_tot;

  return true;
}

int main() {
  using namespace Eigen;

  VectorXd x(13);
  x <<  0.1,
        1.0,
        5.0,
       10.0,
       20.0,
       40.0,
       60.0,
       80.0,
      100.0,
      150.0,
      200.0,
      250.0,
      300.0;

  VectorXd y(13);
  y <<  7569.9,
        7609.6,
        7684.4,
        7778.6,
        8007.3,
        8619.8,
        9461.2,
       10482.6,
       11572.4,
       15021.0,
       19210.4,
       23777.2,
       28426.9;

  VectorXd p;
  double r2;

  if (PolyVal(x, y, 2, p, &r2)) {
    std::cout << "p = " << p.transpose() << std::endl;
    std::cout << "R2 = " << r2 << std::endl;
  }
  if (PolyVal(x, y, 1, p, &r2)) {
    std::cout << "p = " << p.transpose() << std::endl;
    std::cout << "R2 = " << r2 << std::endl;
  }

  return 0;
}

$ clang++ -I/usr/local/Cellar/eigen/3.3.1/include/eigen3 -o main main.cc
$ ./main
p = 0.139124  29.4058  7420.85
R2 = 0.99904
p = 67.3193 6411.87
R2 = 0.972701