如何使用 Eigen 矩阵计算来计算多项式

Question

我想使用矩阵计算而不是循环来计算多项式。

理论

k次多项式的方程：

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

我们可以通过此矩阵计算计算 n X 值的所有 Y 值：

问题

我有所有系数。
我有一个带有 X 值的向量，使用 Eigen::VectorXd::LinSpaced(size, start, stop);

如何用Eigen生成T矩阵？

当前解决方案

现在，我使用两个 for 循环：

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

评估多项式的常用方法是Horner's method。这避免了任何复杂的函数（例如 pow），速度快且数值稳定。

Eigen 中的一个版本可能如下所示：

/** 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 实现了“稳定”版本。我找不到这方面的参考。无论如何，如果我们根据您的输入模式调整该代码，我们会得到：

#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;
    });
}

维基百科列出 several other methods 来计算多项式。

如何使用 Eigen 矩阵计算来计算多项式

How to calculate a polynomial using matrix calculation with Eigen

c++

eigen

理论

问题

当前解决方案