SSE版本的差平方和算法的累积计算误差

Question

我正在尝试优化以下代码（两个数组的平方差之和）：

inline float Square(float value)
{
    return value*value;
}

float SquaredDifferenceSum(const float * a, const float * b, size_t size)
{
    float sum = 0;
    for(size_t i = 0; i < size; ++i)
        sum += Square(a[i] - b[i]);
    return sum;
}

所以我使用 CPU:

的 SSE 指令进行了优化

inline void SquaredDifferenceSum(const float * a, const float * b, size_t i, __m128 & sum)
{
    __m128 _a = _mm_loadu_ps(a + i);
    __m128 _b = _mm_loadu_ps(b + i);
    __m128 _d = _mm_sub_ps(_a, _b);
    sum = _mm_add_ps(sum, _mm_mul_ps(_d, _d));
}

inline float ExtractSum(__m128 a)
{
    float _a[4];
    _mm_storeu_ps(_a, a);
    return _a[0] + _a[1] + _a[2] + _a[3];
}

float SquaredDifferenceSum(const float * a, const float * b, size_t size)
{
    size_t i = 0, alignedSize = size/4*4;
    __m128 sums = _mm_setzero_ps();
    for(; i < alignedSize; i += 4)
        SquaredDifferenceSum(a, b, i, sums);
    float sum = ExtractSum(sums);
    for(; i < size; ++i)
        sum += Square(a[i] - b[i]);
    return sum;
}

如果数组的大小不太大，此代码可以正常工作。但如果尺寸足够大，则基函数给出的结果与其优化版本之间存在很大的计算误差。所以我有一个问题：SSE优化代码哪里有bug，导致计算错误

Answer 1

错误来自有限精度浮点数。两个浮点数的每次相加都会产生与它们之间的差值成比例的计算误差。在你的标量版本的算法中，结果总和比每一项大得多（当然，如果数组的大小足够大）。所以会导致大的计算误差的积累。

在SSE版的算法中其实有4个和用于结果累加。并且这些和和每一项之间的差相对于标量代码小四倍。所以这导致较小的计算误差。

有两种方法可以解决这个错误：

1) 使用双精度浮点数进行累加

2) 使用卡汉求和算法（也称为补偿求和），与明显的方法相比，它显着减少了通过添加有限精度浮点数序列而获得的总数中的数值误差。

https://en.wikipedia.org/wiki/Kahan_summation_algorithm

使用 Kahan 求和算法，您的标量代码将如下所示：

inline void KahanSum(float value, float & sum, float & correction)
{
    float term = value - correction;
    float temp = sum + term;
    correction = (temp - sum) - term;
    sum = temp; 
}

float SquaredDifferenceKahanSum(const float * a, const float * b, size_t size)
{
    float sum = 0, correction = 0;
    for(size_t i = 0; i < size; ++i)
        KahanSum(Square(a[i] - b[i]), sum, correction);
    return sum;
}

SSE 优化后的代码如下所示：

inline void SquaredDifferenceKahanSum(const float * a, const float * b, size_t i, 
                                      __m128 & sum, __m128 & correction)
{
    __m128 _a = _mm_loadu_ps(a + i);
    __m128 _b = _mm_loadu_ps(b + i);
    __m128 _d = _mm_sub_ps(_a, _b);
    __m128 term = _mm_sub_ps(_mm_mul_ps(_d, _d), correction);
    __m128 temp = _mm_add_ps(sum, term);
    correction = _mm_sub_ps(_mm_sub_ps(temp, sum), term);
    sum = temp; 
}

float SquaredDifferenceKahanSum(const float * a, const float * b, size_t size)
{
    size_t i = 0, alignedSize = size/4*4;
    __m128 sums = _mm_setzero_ps(), corrections = _mm_setzero_ps();
    for(; i < alignedSize; i += 4)
        SquaredDifferenceKahanSum(a, b, i, sums, corrections);
    float sum = ExtractSum(sums), correction = 0;
    for(; i < size; ++i)
        KahanSum(Square(a[i] - b[i]), sum, correction);
    return sum;
}

SSE版本的差平方和算法的累积计算误差

An accumulated computing error in SSE version of algorithm of the sum of squared differences

c++

sse

simd