std::vector 的 C++ 高效插值答案

【问题标题】：C++ Efficient interpolation of a std::vectorstd::vector 的 C++ 高效插值
【发布时间】：2021-08-07 04:24:55
【问题描述】：

我需要通过插值找到一个给定未知函数的值。问题是我创建的效率太低了。

首先，我读取了一个包含 y=g(T) 和 T 的数据文件，但采用离散形式。我将它们的值存储在 std::vector<double> 中。

在此之后，我将 T (std::vector<double> Tgdat) 转换为 x (std::vector<double> xgdat)。这将是伴随 y 轴的 x 轴，(std::vector<double> gdat)。

然后，我创建一个函数来插入我的向量 std::vector<double> gdat，这样，给定一些 x（它的值在向量 std::vector<double> xgdat 的两个元素之间），程序可以为 g(x ）。这个函数通过引用接收向量，不是因为我想修改它们（这就是为什么我也将它们传递为const），而是为了让计算机不必创建它的副本。

double geffx (double x, const std::vector<double> &gdat, const std::vector<double> &xgdat)
{
  //Local variables
  double g;
  int k,l;

  //Find the index of the element of xgdat that is nearest to x
  auto i = min_element(xgdat.begin(), xgdat.end(),
      [x] (double a, double b)
      {
          return abs(x-a)<abs(x-b);
      });
  k = std::distance(xgdat.begin(), i); //Nearest index

  //Find the index of the element of xgdat that is nearest to x
  //and it is not the same index as before
  auto j = min_element(xgdat.begin(), xgdat.end(),
      [x,&xgdat,k] (double a, double b)
      {
          if (a!=xgdat[k]) return abs(x-a)<abs(x-b);
          else return false;
      });
  l = std::distance(xgdat.begin(), j); //Second nearest index

  //Interpolation:
  if(xgdat[k]<xgdat[l]) 
      g = gdat[k]+(x-xgdat[k])*(gdat[l]-gdat[k])/(xgdat[l]-xgdat[k]);
  else 
      g = gdat[l]+(x-xgdat[l])*(gdat[k]-gdat[l])/(xgdat[k]-xgdat[l]);

  return g;
}

这似乎效率非常低，但我无法解决解决同样问题但以更有效的方式解决问题的方法。我已经尝试过 const 的事情并且也通过引用传递，但我想最大的问题是 min_element() 函数/方法，也许它也与最后的 if-else 返回值有关g.

编辑：额外信息

我使用g++作为编译器，元素个数是275。

由于此函数是 EDO 求解器的一部分，因此在每一步（1e4 步）中都会多次调用它，直到收敛。我需要 4 个插值器，每个插值器都被多次调用以进行评估，所以我会说该函数需要被访问超过 1e6 次。

当我用一个常数（不需要插值器）替换 g(x) 时，执行时间大约是 1-10 秒。现在是 45 分钟 - 1 小时。（非常糟糕，我知道，这就是我需要帮助的原因）

【问题讨论】：

缺少信息，例如使用的编译器、您是在运行优化构建还是“调试”构建、向量中元素的典型数量以及实际时序统计信息。
好的，我会添加信息。但是，我不知道您所说的“优化构建”与“调试构建”是什么意思，我只是写了典型的g++ name.cpp -o name，然后在终端中写.\name 运行它。
二分查找，既然数据是排序的？（如果不是，先排序）
@AdriánDavid 你没有给它任何优化标志，所以你的代码会很慢！有你的问题，我想。尝试添加-O2。假设调试信息默认包含在 g++ 中，这是一个“调试构建”。
@JDługosz 你只需要对 x 进行排序，不需要单调函数...

标签： c++ sorting interpolation std stdvector

【解决方案1】：

输min_element，输绝对距离比较。它是一种凸变换，可以有效地搜索，但不能通过 C++ 标准库中存在的任何函数进行搜索。

无论如何，您不想要两个最接近的点，您想将您的评估放在上方和下方。（“最接近的两个”总是给出括号对的唯一情况是当样本是均匀间隔的，如果这是真的你根本不需要搜索，你可以直接使用间隔计算索引）。

使用lower_bound 在排序后的数组中对您关心的x 进行高效的二分搜索。那是括号的一侧，另一个索引低一个。

最后，您的代码的顶部将如下所示：

//Find the index of the element of xgdat that is nearest-above to x
auto i = lower_bound(xgdat.begin(), xgdat.end(), x); 
//If the vector values are in decreasing order use:
//auto i = lower_bound(xgdat.rbegin(), xgdat.rend(), x);
k = xgdat.begin() - i; //Nearest index
if (i == xgdat.end())
  --k;  // extrapolating above
else if (*i == x)
  return gdat[k];

l = k? k - 1: 1; //nearest-below index, except when extrapolating downward

// proceed with linear interpolation/extrapolation using l and k

【讨论】：

lower_bound() 似乎只有在xgdat 以升序而不是降序排序时才有效，有时会出现这种情况。 i 总是指向第一个元素，如果它是按降序排列的。
@AdriEscañuela：迭代器必须按升序查找数据。如果您的向量按降序排列，请使用lower_bound(xgdat.rbegin(), xgdat.rend(), x);
@AdriEscañuela：您知道这可能会使您的求解器快多少吗？
由于 EDO 求解器的 jacobian 中的一些错误，它被放慢了。现在它在几分之一秒内就解决了，就像我之前设置 g(x)=constant 时一样。现在正确的程序在我的问题中给出的实现也运行得非常快，但我接受了你的解决方案，因为min_element() 复杂性增长为 $O(n)$ 而lower_bound() 为 $O(log_2(n) )$ （二进制排序更好）。在我的程序中它根本不明显，但有人可能会发现您的解决方案很有帮助。

【解决方案2】：

关于代码的一些提示：

在需要的地方声明变量，而不是全部在顶部。

if(xgdat[k]<xgdat[l]) 
      g = gdat[k]+(x-xgdat[k])*(gdat[l]-gdat[k])/(xgdat[l]-xgdat[k]);
  else 
      g = gdat[l]+(x-xgdat[l])*(gdat[k]-gdat[l])/(xgdat[k]-xgdat[l]);

除了交换 k 和 l 之外，这两行看起来是相同的。所以不要重复该行：只需交换 k 和 l！

if(xgdat[k]>=xgdat[l]) std::swap(k,l);

而g 只在这里使用？为什么要在函数顶部声明它？现在完全不需要了：

return gdat[k]+(x-xgdat[k])*(gdat[l]-gdat[k])/(xgdat[l]-xgdat[k]);

在调用与迭代器一起使用的标准算法后，您将很难恢复“索引位置”。你应该只使用迭代器。

不过，这些都不是整体效率。

【讨论】：