在数值算法的准确性和稳定性中，等式 1.6a，Higham 给出了以下均值更新公式：
好的，这个更新的一个好处是它不会溢出，当我们天真计算时，这是一个潜在的问题然后除以在任何使用之前。
但是天真的总和需要添加和 1 除法，并且可能很好地矢量化，而 Higham 提出的更新需要补充和部门，我不知道有任何汇编指令可以矢量化所有这些。
那么海厄姆的更新公式值得使用吗？有没有我看不到的好处？
注意：Higham 给出了以下通用建议（第 1.18 节“设计稳定算法”）：
“表达更新公式是有利的，就
newvalue = oldvalue + smallcorrection好像可以用许多正确的有效数字计算小的修正一样。”
更新 1.6a 确实采用了这种形式，但我不清楚是否可以对许多有效数字进行小的修正。
编辑：我发现对各种计算方法的性能进行了实证研究，强烈推荐 1.6a；看
但是，在阅读了那篇论文后，我仍然不清楚更新是否物有所值。无论如何，我希望通过累积舍入误差得到最坏和平均的情况。
Higham 的算法对于在线算法或存储容量有限的算法（例如边缘处理）似乎非常有用。在这些情况下，它可能总是物有所值。
但是，为了解决您提出的问题，我实现了一个 SSE2 版本，我认为它可以捕捉到您的问题：
#include <chrono>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <xmmintrin.h>
int main(){
const int num_count = 100000;
alignas(16) float data[num_count];
for(int i=0;i<num_count;i++)
data[i] = rand()%100000+rand()/(double)RAND_MAX;
{
const auto t1 = std::chrono::high_resolution_clock::now();
float sum=0;
for(int i=0;i<num_count;i++){
sum += data[i];
}
const float mean=sum/num_count;
const auto t2 = std::chrono::high_resolution_clock::now();
const auto time_span1 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
std::cout << "Exec time:\t" << time_span1.count() << " s\n";
std::cout << "Mean = "<<std::setprecision(20)<<mean<<std::endl;
}
{
const auto t1 = std::chrono::high_resolution_clock::now();
__m128 mean = _mm_load_ps(&data[0]);
for (int i=4;i<num_count;i+=4){
const __m128 x = _mm_load_ps(&data[i]);
const __m128 diff = _mm_sub_ps(x,mean);
const __m128 k = _mm_set1_ps(i/4);
const __m128 div = _mm_div_ps(diff, k);
mean = _mm_add_ps(mean, div);
}
float result[4];
_mm_store_ps(result, mean);
const float tmean = (result[0] + result[1] + result[2] + result[3])/4; //I'm suspicious about this step: probably throws away precision
const auto t2 = std::chrono::high_resolution_clock::now();
const auto time_span1 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);
std::cout << "Exec time:\t" << time_span1.count() << " s\n";
std::cout << "Mean = "<<std::setprecision(20)<<tmean<<std::endl;
}
}
}
并观察到
Exec time: 0.000225851 s
Simple Mean = 49891.23046875
Exec time: 0.0003759360000000000002 s
Higham Mean = 49890.26171875
Higham 的平均值需要更长的时间来计算，并且值的差异可能是不可忽略的数量，尽管您需要的准确性实际上取决于您的实施。
假设你有数字和向量长度划分， IE，. 您可以将 Higham 方案应用于每组数字。这很好地矢量化，特别是如果数字是交错的，即所有第一个元素在内存中都是连续的，然后是所有第二个元素，等等。当你计算出你拥有的每个组的平均值时剩下的号码。您应用 Higham 方案的顺序实现，然后您就完成了。如果向量长度不除，然后用零填充数据，将数据拆分为团体与每个数字并像以前一样继续。当你有平均值时每组，你计算并使用简单的求和算法简单地添加这些数字。
海厄姆的公式是否值得使用取决于情况。当溢出是一个问题时，天真的公式就没用了。这两个公式的算术强度都很低，我希望 Higham 的公式运行非常缓慢，这仅仅是因为除法的延迟非常高。就个人而言，我发现与可靠性和准确性相比，速度是非常次要的。
Higham 在第一章的主要目标只是提醒读者注意一个事实，即使用浮点数进行计算与真正的算术有很大的不同。他关于更正的建议被正确引用，但我会写得有点不同。具体而言，只要校正相对于原始值较小，则校正是否以较大的相对误差计算无关紧要。