在分工的在 compute_variances（）可以肯定的是沉重的东西！然而，Knuth的算法是一个数值稳定的一个，我没能找到另一种算法，将取消这两个数值稳定，没有分裂。

请注意，我的没有的在有趣的并行的现在。

最小点类（我想我也没有忘记展示的东西）的例子

 类点{
 上市：

  TYPEDEF浮动FT;

  ...

  / **
   *获取点的尺寸。
   * /
  昏暗的size_t（）const的{
    返回coords.size（）;
  }

  / **
   *运营商返回的坐标给定索引处。
   * @参数我 - 的坐标指数
   * @返回的坐标索引i
   * /
  FT＆安培;运营商[]（const int的我）{
    返回coords.at（ⅰ）;
    //这是同样的，如果我有下面的评论code
    / *如果（ⅰ≤（INT）coords.size（）＆安培;＆安培; I＆GT; = 0）
      返回coords.at（ⅰ）;
    其他 {
      性病::法院＆LT;＆LT; 在点错误:: []＆LT;＆LT;的std :: ENDL;
      出口（1）;
    }
    返回的coords [0]; //清除-Wall警告* /
  }

  / **
   *运营商返回的坐标给定索引处。 （不变）
   * @参数我 - 的坐标指数
   * @返回的坐标索引i
   * /
  常量FT＆放大器;运营商[]（const int的我）const的{
        返回coords.at（ⅰ）;
    / *如果（ⅰ≤（INT）coords.size（）＆安培;＆安培; I＆GT; = 0）
      返回coords.at（ⅰ）;
    其他 {
      性病::法院＆LT;＆LT; 在点错误:: []＆LT;＆LT;的std :: ENDL;
      出口（1）;
    }
    返回的coords [0]; //清除-Wall警告* /
  }

 私人：
  的std ::矢量＆lt; FT＆GT; COORDS;
};
 

解决方案

下面是我会怎么做，在重要性guesstimated顺序：

在点::运算符[] 将得到内联的发行版本的机会，并能减少到2个机器指令（有效地址计算和浮点负载） 。这就是你想要的。当然，它必须是定义的头文件中，否则会内联仅当您启用链接时执行code代（又名LTO）。

注意点::运算符[] 的身体只相当于单行 返回coords.at（我）在调试版本。在发布版本中的全部的身体就等于返回COORDS [I] ，不可以 返回coords.at（我）。

I can see with the CPU profiler, that the compute_variances() is the bottleneck of my project.

  %   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*)
 ...

Here is the body of the function:

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

where kthLargest() doesn't seem to be a problem, since I see that:

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

The compute_variances() takes a vector of vectors of floats (i.e. a vector of Points, where Points is a class I have implemented) and computes the variance of them, in each dimension (with regard to the algorithm of Knuth).

Here is how I call the function:

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

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

The question is, can I do better? I would really happy to pay the trade-off between speed and approximate computation of variances. Or maybe I could make the code more cache friendly or something?

I compiled like this:

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

Notice, that before edit, I had used -o3, not with a capital 'o'. Thanks to ypnos, I compiled now with the optimization flag -O3. I am sure that there was a difference between them, since I performed time measurements with one of these methods in my pseudo-site.

Note that now, compute_variances is dominating the overall project's time!

[EDIT]

copute_variances() is called 40 times.

Per 10 calls, the following hold true:

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

Each call handles different data.

Q: How fast is access to points[i][d]?

A: point[i] is just the i-th element of std::vector, where the second [], is implemented as this, in the Point class.

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 
}

where coords is a std::vector of float values. This seems a bit heavy, but shouldn't the compiler be smart enough to predict correctly that the branch is always true? (I mean after the cold start). Moreover, the std::vector.at() is supposed to be constant time (as said in the ref). I changed this to have only .at() in the body of the function and the time measurements remained, pretty much, the same.

The division in the compute_variances() is for sure something heavy! However, Knuth's algorithm was a numerical stable one and I was not able to find another algorithm, that would de both numerical stable and without division.

Note that I am not interesting in parallelism right now.

[EDIT.2]

Minimal example of Point class (I think I didn't forget to show something):

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

Here's what I would do, in guesstimated order of importance:

The Point::operator[] will have a chance of getting inlined in the release build and can reduce to two machine instructions (effective address computation and floating point load). That's what you want. Of course it must be defined in the header file, otherwise the inlining will only be performed if you enable link-time code generation (a.k.a. LTO).

Note that the Point::operator[]'s body is only equivalent to the single-line return coords.at(i) in a debug build. In a release build the entire body is equivalent to return coords[i], not 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);
}