How to calculate a polynomial using matrix calculation with Eigen

Question

I want to calculate a polynomial using matrix calculation and not for loops.

Theory

The equation of a polynomial of degree k:

a: Coeficients of the polynomial
t: X value
v: Y value to calculate

We can calutate all Y values for n X values with this matrix calculation:

Question

• I have all coeficients.
• I have a vector with X values using Eigen::VectorXd::LinSpaced(size, start, stop);

How to generate the T matrix with Eigen?

Current solution

For now, I'm using two for loops:

std::vector<double> yVector;
Eigen::VectorXd xVector = Eigen::VectorXd::LinSpaced(40, -20, 20);

for(const double x : xVector)
{
  double y= 0;
  for(size_t i = 0; i < coeff.size(); i++)
  {
    y+= coeff[i] * pow(x, i);
  }
  yVector.push_back(y);
}

Answer 1

The usual approach to evaluate polynomials would be Horner's method. This avoids any complex functions (such as pow) , is fast and numerically stable.

A version in Eigen could look something like this:

/** Coefficients ordered from highest to lowest */
Eigen::VectorXd evaluate_poly(
      const Eigen::Ref<const Eigen::VectorXd>& xvals,
      const Eigen::Ref<const Eigen::VectorXd>& coeffs)
{
    auto coeff = std::begin(coeffs);
    const auto last_coeff = std::end(coeffs);
    assert(coeff != last_coeff && "Empty set of polynomial coefficients");
    Eigen::VectorXd yvals =
          Eigen::VectorXd::Constant(xvals.size(), *coeff);
    for(++coeff; coeff != last_coeff; ++coeff)
        yvals = yvals.array() * xvals.array() + *coeff;
    return yvals;
}

Eigen's unsupported Polynomial module implements a "stabilized" version. I wasn't able to find a reference for this. In any case, if we adapt that code to your input pattern, we get this:

#include <cmath>
// using std::abs, std::pow
#include <iterator>
// using std::make_reverse_iterator

Eigen::VectorXd evaluate_poly(
      const Eigen::Ref<const Eigen::VectorXd>& xvals,
      const Eigen::Ref<const Eigen::VectorXd>& coeffs)
{
    assert(coeffs.size() && "Empty set of polynomial coefficients");
    return xvals.unaryExpr([&coeffs](double x) noexcept -> double {
          auto coeff = std::begin(coeffs);
          const auto last_coeff = std::end(coeffs);
          double y;
          if(! (std::abs(x) > 1.) /*NaN or range [-1, 1] */) {
              // normal Horner method
              for(y = *(coeff++); coeff != last_coeff; ++coeff)
                  y = y * x + *coeff;
              return y;
          }
          const double inv_x = 1. / x;
          auto reverse_coeff = std::make_reverse_iterator(last_coeff);
          const auto last_reverse = std::make_reverse_iterator(coeff);
          for(y = *(reverse_coeff++); reverse_coeff != last_reverse; ++reverse_coeff)
              y = y * inv_x + *reverse_coeff;
          return std::pow(x, coeffs.size() - 1) * y;
    });
}

Wikipedia lists several other methods to evaluate polynomials.