为什么这段代码这么慢？

Question

所以我有这个函数用于计算统计数据（最小值/最大值/标准/平均值）。现在问题是这通常在 10,000 x 15,000 矩阵上运行。矩阵在类中存储为vector<vector<int> >。现在创建和填充所述矩阵非常快，但是当涉及到统计部分时，它变得非常慢。

例如一次读取一个像素的 geotiff 的所有像素值大约需要 30 秒。 （这涉及到大量复杂的数学运算，以将像素值正确地定位到对应点），计算整个矩阵的统计数据大约需要 6 分钟。

void CalculateStats()
{
    //OHGOD
    double new_mean = 0;
    double new_standard_dev = 0;

    int new_min = 256;
    int new_max = 0;

    size_t cnt = 0;
    for(size_t row = 0; row < vals.size(); row++)
    {
        for(size_t col = 0; col < vals.at(row).size(); col++)
        {
            double mean_prev = new_mean;
            T value = get(row, col);
            new_mean += (value - new_mean) / (cnt + 1);
            new_standard_dev += (value - new_mean) * (value - mean_prev);

            // find new max/min's
            new_min = value < new_min ? value : new_min;
            new_max = value > new_max ? value : new_max;
            cnt++;
        }
    }

    stats_standard_dev = sqrt(new_standard_dev / (vals.size() * vals.at(0).size()) + 1);
    std::cout << stats_standard_dev << std::endl;
}

我在这里做了什么可怕的事情吗？

编辑

为了响应 cmets，T 将是一个 int。

编辑 2

我修正了我的标准算法，这是最终产品：

void CalculateStats(const std::vector<double>& ignore_values)
{
    //OHGOD
    double new_mean = 0;
    double new_standard_dev = 0;

    int new_min = 256;
    int new_max = 0;

    size_t cnt = 0;

    int n = 0;
    double delta = 0.0;
    double mean2 = 0.0;

    std::vector<double>::const_iterator ignore_begin = ignore_values.begin();
    std::vector<double>::const_iterator ignore_end = ignore_values.end();

    for(std::vector<std::vector<T> >::const_iterator row = vals.begin(), row_end = vals.end();  row != row_end; ++row)
    {
        for(std::vector<T>::const_iterator col = row->begin(), col_end = row->end(); col != col_end; ++col)
        {
            // This method of calculation is based on Knuth's algorithm.
            T value = *col;
            if(std::find(ignore_begin, ignore_end, value) != ignore_end)
                continue;
            n++;
            delta = value - new_mean;
            new_mean = new_mean + (delta / n);
            mean2 = mean2 + (delta * (value - new_mean));

            // Find new max/min's.
            new_min = value < new_min ? value : new_min;
            new_max = value > new_max ? value : new_max;
        }
    }
    stats_standard_dev = mean2 / (n - 1);
    stats_min = new_min;
    stats_max = new_max;
    stats_mean = new_mean;

这仍然需要大约 120-130 秒才能完成，但这是一个巨大的改进 :)！

Answer 1

您是否尝试过分析您的代码？

您甚至不需要花哨的分析器。只需在其中粘贴一些调试时间语句即可。

我告诉你的任何事情都只是一个有根据的猜测（而且可能是错误的）

由于您访问向量内容的方式，您可能会遇到很多缓存未命中。您可能想将一些结果缓存到 size() 但我不知道这是否是问题所在。

Answer 2

我刚刚对其进行了分析。 90% 的执行时间都在这一行：

new_mean += (value - new_mean) / (cnt + 1);

Answer 3

您应该在第一个循环中计算值、最小值、最大值和计数的总和， 然后通过将总和/计数除以计算一次操作中的平均值， 然后在第二个循环中计算 std_dev 的总和

那可能会快一点。

Answer 4

我发现的第一件事是您在循环中评估vals.at(row).size()，显然这不应该提高性能。它也适用于vals.size()，但当然内循环更糟。如果vals 是向量的向量，则最好使用迭代器或至少保留外部向量的引用（因为带有索引参数的get() 肯定也会占用相当长的时间）。

这个代码示例应该说明我的意图；-)

for(TVO::const_iterator i=vals.begin(),ie=vals.end();i!=ie;++i) {
    for(TVI::const_iterator ii=i->begin(),iie=i->end();ii!=iie;++ii) {
        T value = *ii;
        // the rest
    }
}

Answer 5

• 首先，将您的 row++ 更改为 ++row。一件小事，但你想要速度，这会有所帮助
• 其次，将您的行
• 中间的“get”方法是什么？那有什么作用？这可能是你真正的问题。
• 我不太确定您的标准偏差计算。看看wikipedia page on calculating variance in a single pass（他们对 Knuth 的算法进行了快速解释，这是递归关系的扩展）。

Answer 6

速度很慢，因为您正在对调试代码进行基准测试。

使用 VS2008 在 Windows XP 上构建和运行代码：

• 具有默认优化级别的发布版本，OP 中的代码运行时间为 2734 毫秒。
• 默认为无优化的调试版本，OP 中的代码运行时间为 398,531 毫秒。

在下面的 cmets 中，你说你没有使用优化，这似乎在这种情况下产生了很大的不同 - 通常它不到十倍，但在这种情况下它慢了一百多倍。

我使用的是 VS2008 而不是 2005，但可能类似：

在 Debug 构建中，每次访问都有两个范围检查，每个检查都使用非内联函数调用 std::vector::size() 并需要分支预测。函数调用和分支都涉及开销。

在 Release 构建中，编译器优化掉了范围检查（我不知道它是直接丢弃它们，还是根据循环的限制进行流分析），向量访问变成了少量的内联没有分支的指针算法。

没有人关心调试构建的速度有多快。无论如何，您都应该对发布版本进行单元测试，因为这是必须正常工作的版本。如果您在尝试单步调试代码时没有提供所需的所有信息，请仅使用 Debug 构建。

---

发布的代码在我的 PC 上运行时间

虽然还有其他建议可以让它更快（例如移动它以使用 SSE，如果所有值都可以使用 8 位类型表示），但显然还有其他一些东西让它变慢了。

在我的机器上，无论是提升对行向量的引用并提升行大小的版本，还是使用迭代器的版本都没有任何可衡量的好处（使用 g++ -O3 使用迭代器重复需要 1511 毫秒；吊装版和原始版都需要 1485 毫秒）。不优化意味着它在 7487ms（原始）、3496ms（提升）或 5331ms（迭代器）中运行。

但是，除非您在非常低功耗的设备上运行，或者正在分页，或者正在运行带有调试器的未优化代码，否则它不应该这么慢，而且任何让它变慢的东西都不太可能成为您发布的代码。

（附带说明，如果您使用偏差为零的值对其进行测试，您的 SD 会显示为 1）

Answer 7

内部循环中的计算太多了：

1. 对于描述性统计（平均值、标准 偏差）唯一需要的就是计算总和 value 的平方和 value 的平方和。从这些 可以计算平均值和标准差的两个总和 在外循环之后（连同第三个值， 样本数量 - n 是您的新/更新代码）。这 方程可以从定义中导出或找到 在网络上，例如维基百科。例如平均值是 只是 value 除以 n 的总和。对于 n 版本（在 与 n-1 版本相比 - 但是 n 在 在这种情况下，使用哪一个都没有关系） 标准差为：
   sqrt( n * sumOfSquaredValue - 总和值 * 总和值）。
   
   所以只有两个浮点数 需要加法和一乘法 内循环。溢出不是这些总和的问题，因为 双打的范围是 10^318。特别是你会 摆脱expensive floating point division 另一个答案中报告的分析已经显示。

2. 一个较小的问题是最小值和最大值是 每次都重写（编译器可能会也可能不会 防止这种情况）。随着最小值迅速变小并且 最大值很快变大，只有两个比较 大多数循环迭代都应该发生：使用 if 语句来确定。可以争论，但是 另一方面，这样做是微不足道的。

Answer 8

我会改变访问数据的方式。假设您将 std::vector 用于您的容器，您可以执行以下操作：

vector<vector<T> >::const_iterator row;
vector<vector<T> >::const_iterator row_end = vals.end();
for(row = vals.begin(); row < row_end; ++row)
{
    vector<T>::const_iterator value;
    vector<T>::const_iterator value_end = row->end();
    for(value = row->begin(); value < value_end; ++value)
    {
        double mean_prev = new_mean;
        new_mean += (*value - new_mean) / (cnt + 1);
        new_standard_dev += (*value - new_mean) * (*value - mean_prev);

        // find new max/min's
        new_min = min(*value, new_min);
        new_max = max(*value, new_max);
        cnt++;
    }
}

这样做的好处是，在您的内部循环中，您不会咨询外部vector，而只是咨询内部循环。

如果您的容器类型是列表，这将明显更快。因为 get/operator[] 的查找时间对于列表是线性的，对于向量是恒定的。

编辑，我将调用 end() 移出循环。

Answer 9

将 .size() 调用移到每个循环之前，并确保您在编译时启用了优化。

Answer 10

如果您的矩阵存储为向量的向量，那么在外部 for 循环中您应该直接检索第 i 个向量，然后在内部循环中对其进行操作。试试看它是否能提高性能。

Answer 11

我不确定vals 是什么类型，但如果vals.at(row).size() 本身遍历集合，可能需要很长时间。将该值存储在变量中。否则它可能会使算法更像 O(n³) 而不是 O(n²)

Answer 12

我认为我会重写它以使用 const 迭代器而不是 row 和 col 索引。我会为 row_end 和 col_end 设置一个 const const_iterator 进行比较，以确保它不会在每个循环结束时进行函数调用。

Answer 13

正如人们所提到的，它可能是get()。例如，如果它访问邻居，您将完全破坏缓存，这将大大降低性能。您应该分析，或者只是考虑访问模式。

Answer 14

在这里聚会有点晚了，但有几点：

1. 您在这里有效地进行了数字工作。我对数值算法了解不多，但我知道参考资料和专家支持通常很有用。这个discussion thread 提供了一些参考；和数字食谱是标准（如果过时）作品。

2. 如果您有机会重新设计您的矩阵，您想尝试使用 valarray 和切片而不是向量的向量；立即浮现在脑海中的一个优势是保证了平面线性布局，这使得缓存预取和 SIMD 指令（如果您的编译器可以使用它们）更有效。

Answer 15

在内部循环中，您不应该测试大小，也不应该进行任何除法，而且迭代器的成本也可能很高。事实上，在那里进行一些展开会很好。 当然，你应该注意缓存的局部性。

如果循环开销足够低，那么在不同的通道中执行它可能是有意义的：一个得到总和（你除以得到平均值），一个得到平方和（你结合得到方差的总和），一个得到最小值和/或最大值。原因是为了简化内部展开循环中的内容，以便编译器可以将内容保存在寄存器中。

我无法编译代码，因此无法确定问题所在。

Answer 16

我已经修改了算法以摆脱几乎所有的浮点除法。

警告：未经测试的代码！！！

void CalculateStats()
{
    //OHGOD

    double accum_f;
    double accum_sq_f;
    double new_mean = 0;
    double new_standard_dev = 0;

    int new_min = 256;
    int new_max = 0;

    const int oku = 100000000;
    int accum_ichi = 0;
    int accum_oku = 0;
    int accum_sq_ichi = 0;
    int accum_sq_oku = 0;

    size_t cnt = 0;
    int v1 = 0;
    int v2 = 0;

    v1 = vals.size();

    for(size_t row = 0; row < v1; row++)
    {

            v2 = vals.at(row).size();
            for(size_t col = 0; col < v2; col++)
            {
                    T value = get(row, col);
                    int accum_ichi += value;
                    int accum_sq_ichi += (value * value);

                    // perform carries
                    accum_oku += (accum_ichi / oku);
                    accum_ichi %= oku;
                    accum_sq_oku += (accum_sq_ichi / oku);
                    accum_sq_ichi %= oku;

                    // find new max/min's
                    new_min = value < new_min ? value : new_min;
                    new_max = value > new_max ? value : new_max;
                    cnt++;
            }
    }

    // now, and only now, do we use floating-point arithmetic
    accum_f = (double)(oku) * (double)(accum_oku) + (double)(accum_ichi);
    accum_sq_f = (double)(oku) * (double)(accum_sq_oku) + (double)(accum_sq_ichi);

    new_mean = accum_f / (double)(cnt);

    // standard deviation formula from Wikipedia
    stats_standard_dev = sqrt((double)(cnt)*accum_sq_f - accum_f*accum_f)/(double)(cnt);        

    std::cout << stats_standard_dev << std::endl;
}