更快地计算(近似)方差

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()需要浮子向量的向量(即一个向量Points,其中,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.我确信它们之间存在差异,因为我在伪站点中使用这些方法之一进行了时间测量.

请注意,现在,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]？

答:point[i]只是std :: vector的第i个元素,其中第二个[]是在Point类中实现的.

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是std::vector的float值.这看起来有点沉重,但编译器是否应该足够聪明才能正确预测分支始终为真？(我的意思是冷启动后).而且,std::vector.at()应该是恒定时间(如参考文献中所述).我把它改成只有.at()功能的主体和时间测量保持,几乎相同.

该师在compute_variances()是肯定的沉重的东西!然而,Knuth的算法是一个数值稳定的算法,我无法找到另一种算法,即数值稳定且没有除法.

请注意,我现在对并行性并不感兴趣.

[EDIT.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. 从值返回浮点Point::operator[]，而不是通过引用。
2. 使用coords[i]而不是coords.at(i)，因为您已经断言它在范围内。成员at检查边界。您只需检查一次。
3. Point::operator[]用断言替换自制的错误指示/检查。这就是断言的用途。它们在发布模式下名义上是无操作的 - 我怀疑您是否需要在发布代码中检查它。
4. 将重复除法替换为一次除法和重复乘法。
5. 通过展开外循环的前两次迭代，消除浪费的初始化的需要。
6. 为了减少缓存未命中的影响，请交替向前和向后运行内部循环。这至少让您有机会使用一些缓存的avg和var. 事实上，它可能会删除所有缓存未命中avg，并且var如果预取按迭代的相反顺序工作，那么它也应该如此。
7. 在现代 C++ 编译器上，std::fill和std::copy可以利用类型对齐，并且有机会比 C 库memset和更快memcpy。

这Point::operator[]将有机会在发布版本中内联，并且可以减少到两个机器指令（有效地址计算和浮点加载）。那就是你想要的。当然，它必须在头文件中定义，否则只有在启用链接时代码生成（又名 LTO）时才会执行内联。

请注意，的主体仅相当于调试版本中的Point::operator[]单行 。return coords.at(i)在一次发布构建中相当于整个return coords[i]身体，不是 return coords.at(i)。

FT Point::operator[](int i) const {
  assert(i >= 0 && i < (int)coords.size());
  return coords[i];
}

const FT * Point::constData() const {
  return &coords[0];
}

void compute_variances(size_t t, const std::vector<Point>& points, float* avg,
                       float* var, size_t* split_dims)
{
  assert(points.size() > 0);
  const int D = points[0].dim();

  // i = 0, i_n = 1
  assert(D > 0);
#if __cplusplus >= 201103L
  std::copy_n(points[0].constData(), D, avg);
#else
  std::copy(points[0].constData(), points[0].constData() + D, avg);
#endif

  // i = 1, i_n = 0.5
  if (points.size() >= 2) {
    assert(points[1].dim() == D);
    for (int d = D - 1; d >= 0; --d) {
      float const delta = points[1][d] - avg[d];
      avg[d] += delta * 0.5f;
      var[d] = delta * (points[1][d] - avg[d]);
    }
  } else {
    std::fill_n(var, D, 0.0f);
  }

  // i = 2, ...
  for (size_t i = 2; i < points.size(); ) {
    {
      const float i_n = 1.0f / (1.0f + i);
      assert(points[i].dim() == D);
      for (int d = 0; d < D; ++d) {
        float const delta = points[i][d] - avg[d];
        avg[d] += delta * i_n;
        var[d] += delta * (points[i][d] - avg[d]);
      }
    }
    ++ i;
    if (i >= points.size()) break;
    {
      const float i_n = 1.0f / (1.0f + i);
      assert(points[i].dim() == D);      
      for (int d = D - 1; d >= 0; --d) {
        float const delta = points[i][d] - avg[d];
        avg[d] += delta * i_n;
        var[d] += delta * (points[i][d] - avg[d]);
      }
    }
    ++ i;
  }

  /* Find t dimensions with largest scaled variance. */
  kthLargest(var, D, t, split_dims);
}