C++ Eigen : spline derivatives() gives strange derivatives

Question

我正在尝试了解如何在 Eigen 中使用样条曲线，特别是我想在某个点上找到样条插值及其一阶和二阶导数的值。找到内插值很容易，但是当我尝试计算导数时，我得到了奇怪的值。

我尝试按照手册 ( http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1Spline.html#af3586ab1929959e0161bfe7da40155c6 ) 中关于 derivatives 命令的说明进行操作，这是我在代码中的尝试:

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

using namespace Eigen;
using namespace std;

double scaling(double x, double min, double max)  // for scaling numbers
{
    return (x - min)/(max - min);
}

VectorXd scale(VectorXd xvals)  // for scaling vectors
{
    const double min = xvals.minCoeff();
    const double max = xvals.maxCoeff();

    for (int k = 0; k < xvals.size(); k++)
        xvals(k) = scaling(xvals(k),min,max);

    return xvals;
}

int main()
{
    typedef Spline<double,1,3> spline;

    VectorXd xvals = (VectorXd(4) << 0,1,2,4).finished();
    VectorXd yvals = xvals.array().square();  // x^2

    spline testspline = SplineFitting<spline>::Interpolate(yvals.transpose(), 3,
scale(xvals).transpose());

    cout << "derivative at x = 0: " << testspline.derivatives(0.00,2) << endl;
    cout << "derivative at x = 1: " << testspline.derivatives(0.25,2) << endl;
    cout << "derivative at x = 2: " << testspline.derivatives(0.50,2) << endl;
    cout << "derivative at x = 3: " << testspline.derivatives(0.75,2) << endl;
    cout << "derivative at x = 4: " << testspline.derivatives(1.00,2) << endl;
}

输出

derivative at x = 0:  0  0 32
derivative at x = 1:  1  8 32
derivative at x = 2:  4 16 32
derivative at x = 3:  9 24 32
derivative at x = 4: 16 32 32

也就是说，插值是正确的(c.f. x = 3)，但导数不是，它们以系统的方式关闭，所以我想我做错了什么。由于这些遵循 x^2，导数应为 0,2,4,6,8，二阶导数应为 2。

关于如何解决这个问题有什么想法吗？

编辑 1

将 x^2 更改为 x^2 + 1 会产生相同的导数，因此至少可以检查。但是将 x^2 更改为 x^3 是错误的，但错误的方式略有不同，输出将是:

derivative at x = 2:  8 48  192
derivative at x = 3: 27 108 288
derivative at x = 4: 64 192 384

哪个是错误的，应该是 6, 9, 12。

同时运行 x^2 案例，但将输入 vector 更改为 0,1,...9 会产生与使用原始输入 vector 相同的导数，但是二阶导数变为稳定的 200，这也是错误的。我不明白为什么二阶导数应该取决于输入点的数量。

Answer 1

解决了。你们非常亲密。您所要做的就是扩展衍生品 与

• 1/(x_max - x_min)(一阶导数)
• 1/(x_max - x_min)^2(二阶导数)。

TLDR:您在拟合样条曲线时将 x 值标准化为介于 0 和 1 之间，但没有缩放 y 值。

您实际拟合的不是样条拟合 x^2，而是:

x_norm = (x - x_min) / (x_max - x_min)
y = x_norm**2

因此使用链式法则，y = x_norm**2 的一阶导数为 2x/(x_max - x_min)，二阶导数为 2/(x_max - x_min)**2。

完整示例代码:

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

using namespace Eigen;
using namespace std;

VectorXd normalize(const VectorXd &x) {
  VectorXd x_norm;
  x_norm.resize(x.size());

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

  return x_norm;
}

int main() {
  typedef Spline<double, 1, 3> Spline1D;
  typedef SplineFitting<Spline1D> Spline1DFitting;

  const Vector4d x{0, 1, 2, 4};
  const Vector4d y = (x.array().square());  // x^2

  const auto knots = normalize(x);  // Normalize x to be between 0 and 1
  const double scale = 1 / (x.maxCoeff() -  x.minCoeff());
  const double scale_sq = scale * scale;
  Spline1D spline = Spline1DFitting::Interpolate(y.transpose(), 3, knots);

  cout << "1st deriv at x = 0: " << spline.derivatives(0.00, 1)(1) * scale << endl;
  cout << "1st deriv at x = 1: " << spline.derivatives(0.25, 1)(1) * scale << endl;
  cout << "1st deriv at x = 2: " << spline.derivatives(0.50, 1)(1) * scale << endl;
  cout << "1st deriv at x = 3: " << spline.derivatives(0.75, 1)(1) * scale << endl;
  cout << "1st deriv at x = 4: " << spline.derivatives(1.00, 1)(1) * scale << endl;
  cout << endl;

  /**
   * 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.
   */

  cout << "2nd deriv at x = 0: " << spline.derivatives(0.00, 2)(2) * scale_sq << endl;
  cout << "2nd deriv at x = 1: " << spline.derivatives(0.25, 2)(2) * scale_sq << endl;
  cout << "2nd deriv at x = 2: " << spline.derivatives(0.50, 2)(2) * scale_sq << endl;
  cout << "2nd deriv at x = 3: " << spline.derivatives(0.75, 2)(2) * scale_sq << endl;
  cout << "2nd deriv at x = 4: " << spline.derivatives(1.00, 2)(2) * scale_sq << endl;

  return 0;
}

示例输出:

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