c++ - 需要更快地计算(近似)方差

Question

我可以通过 CPU 分析器看到，compute_variances() 是我项目的瓶颈。

  %   cumulative   self              self     total           
 time   seconds   seconds    calls  ms/call  ms/call  name    
 75.63      5.43     5.43       40   135.75   135.75  compute_variances(unsigned int, std::vector<Point, std::allocator<Point> > const&, float*, float*, unsigned int*)
 19.08      6.80     1.37                             readDivisionSpace(Division_Euclidean_space&, char*)
 ...

这是函数体:

void compute_variances(size_t t, const std::vector<Point>& points, float* avg,
                       float* var, size_t* split_dims) {
  for (size_t d = 0; d < points[0].dim(); d++) {
    avg[d] = 0.0;
    var[d] = 0.0;
  }
  float delta, n;
  for (size_t i = 0; i < points.size(); ++i) {
    n = 1.0 + i;
    for (size_t d = 0; d < points[0].dim(); ++d) {
      delta = (points[i][d]) - avg[d];
      avg[d] += delta / n;
      var[d] += delta * ((points[i][d]) - avg[d]);
    }
  }

  /* Find t dimensions with largest scaled variance. */
  kthLargest(var, points[0].dim(), t, split_dims);
}

kthLargest() 似乎不是问题，因为我看到了:

0.00 7.18 0.00 40 0.00 0.00 kthLargest(float*, int, int, unsigned int*)

compute_variances() 采用浮点 vector 的 vector (即 Points 的 vector ，其中 Points 是我实现的类) 并计算它们在每个维度上的方差(关于 Knuth 的算法)。

下面是我调用函数的方式:

float avg[(*points)[0].dim()];
float var[(*points)[0].dim()];
size_t split_dims[t];

compute_variances(t, *points, avg, var, split_dims);

问题是，我能做得更好吗？我真的很乐意在速度和方差的近似计算之间进行权衡。或者也许我可以使代码对缓存更友好什么的？

我是这样编译的:

g++ main_noTime.cpp -std=c++0x -p -pg -O3 -o eg

请注意，在编辑之前，我使用了 -o3，而不是大写的“o”。感谢 ypnos，我现在使用优化标志 -O3 进行编译。我确信它们之间存在差异，因为我使用 these methods 之一进行了时间测量。在我的伪网站中。

请注意，现在 compute_variances 支配着整个项目的时间!

[编辑]

copute_variances() 被调用了 40 次。

每 10 次调用，以下情况成立:

points.size() = 1000   and points[0].dim = 10000
points.size() = 10000  and points[0].dim = 100
points.size() = 10000  and points[0].dim = 10000
points.size() = 100000 and points[0].dim = 100

每个调用处理不同的数据。

问:访问 points[i][d] 的速度有多快？

A:point[i] 只是 std::vector 的第 i 个元素，其中第二个 [] 是这样实现的，在 点类。

const FT& operator [](const int i) const {
  if (i < (int) coords.size() && i >= 0)
     return coords.at(i);
  else {
     std::cout << "Error at Point::[]" << std::endl;
     exit(1);
  }
  return coords[0];  // Clear -Wall warning 
}

其中 coords 是 float 值的 std::vector。这似乎有点沉重，但编译器不应该足够聪明以正确预测分支总是正确的吗？ (我的意思是在冷启动之后)。此外，std::vector.at() 应该是常数时间(如 ref 中所述)。我将其更改为函数主体中只有 .at()，并且时间测量值几乎保持不变。

compute_variances() 中的除法 肯定很重!然而，Knuth 的算法是一个数值稳定的算法，我无法找到另一种算法，既能实现数值稳定又不除法。

请注意，我现在不对并行感兴趣。

[编辑.2]

Point 类的最小示例(我想我没有忘记展示一些东西):

class Point {
 public:

  typedef float FT;

  ...

  /**
   * Get dimension of point.
   */
  size_t dim() const {
    return coords.size();
  }

  /**
   * Operator that returns the coordinate at the given index.
   * @param i - index of the coordinate
   * @return the coordinate at index i
   */
  FT& operator [](const int i) {
    return coords.at(i);
    //it's the same if I have the commented code below
    /*if (i < (int) coords.size() && i >= 0)
      return coords.at(i);
    else {
      std::cout << "Error at Point::[]" << std::endl;
      exit(1);
    }
    return coords[0];  // Clear -Wall warning*/
  }

  /**
   * Operator that returns the coordinate at the given index. (constant)
   * @param i - index of the coordinate
   * @return the coordinate at index i
   */
  const FT& operator [](const int i) const {
        return coords.at(i);
    /*if (i < (int) coords.size() && i >= 0)
      return coords.at(i);
    else {
      std::cout << "Error at Point::[]" << std::endl;
      exit(1);
    }
    return coords[0];  // Clear -Wall warning*/
  }

 private:
  std::vector<FT> coords;
};

Answer 1

<强>1。 SIMD

一个简单的加速方法是使用 vector 指令 (SIMD) 进行计算。在 x86 上，这意味着 SSE、AVX 指令。根据您的字长和处理器，您可以获得大约 x4 甚至更多的加速。这里的代码:

for (size_t d = 0; d < points[0].dim(); ++d) {
    delta = (points[i][d]) - avg[d];
    avg[d] += delta / n;
    var[d] += delta * ((points[i][d]) - avg[d]);
}

可以通过使用 SSE 一次计算四个元素来加快速度。由于您的代码在每次循环迭代中实际上只处理一个元素，因此不存在瓶颈。如果你使用 16 位 short 而不是 32 位 float (当时的近似值)，你可以在一条指令中容纳八个元素。如果使用 AVX，它会更多，但你需要一个最新的处理器。

它不是您的性能问题的解决方案，而只是其中的一个，也可以与其他解决方案结合使用。

<强>2。微并行化

当你有那么多循环时，第二个简单的加速是使用并行处理。我通常使用 Intel TBB，其他人可能会建议使用 OpenMP。为此，您可能必须更改循环顺序。因此，在外循环中对 d 进行并行化，而不是对 i 进行并行化。

您可以将这两种技术结合起来，如果操作正确，在具有 HT 的四核上，您可以将组合的速度提高 25-30，而不会损失任何准确性。

<强>3。编译器优化

首先，也许这只是 SO 上的错字，但它需要是 -O3，而不是 -o3! 作为一般注意事项，如果您在实际使用它们的范围内声明变量 delta, n，编译器可能更容易优化您的代码。您还应该尝试 -funroll-loops 编译器选项以及 -march。后者的选项取决于你的 CPU，但现在通常 -march core2 很好(也适用于最近的 AMD)，并且包括 SSE 优化(但我不相信编译器还没有这样做你的循环)。