SIMD拡張命令を用いた関数と普通の関数とで結果が異なってしまう

Question

質問内容

私は現在simd拡張命令(AVX 256bit)を使用して
大きなベクトル同士の内積を高速化しようとしており、下記コードを作成し
通常の内積関数とsimd版内積関数とで検算を行いました。

しかし結果は下記実行結果のようになり
結果が合いませんでした。

なぜ計算結果が合わないのでしょうか。

色々考えたり調べてみましたがそれらしい情報に
有りつけなかったのでご質問させていただきます。

どうか御回答よろしくお願いいたします.

問題の発生するコード

#include <iostream>
#include <cstdint>
#include <immintrin.h>
#include <iomanip>

inline size_t round_up(const size_t close_to_num, const size_t divisible_by)noexcept{//    divisible_byで割り切れる最もclose_to_numに近い整数を算出
    return divisible_by == 0 ? close_to_num : close_to_num + divisible_by - (close_to_num % divisible_by);
}

float dot_avx(const int n, const float* const a, const float* const b)noexcept{//   simdでの内積計算
    __m256 ymm_reg = _mm256_setzero_ps();
    
    for(size_t i = 0; i < n; i += 8) ymm_reg = _mm256_add_ps(ymm_reg, _mm256_mul_ps(_mm256_load_ps(a + i), _mm256_load_ps(b + i)));
    
    alignas(32) float array[8];
    _mm256_store_ps(array, ymm_reg);
    
    return array[0] + array[1] + array[2] + array[3] + array[4] + array[5] + array[6] + array[7];
}

float dot(const int n, const float* const a, const float* const b)noexcept{//   simdを使わず内積計算
    float sum = 0;
    for(int i = 0; i < n; ++i) sum += a[i] * b[i];
    return sum;
}

int main(){
    
    size_t apparent_size = 1234;// 見かけ上のサイズ
    size_t true_size;// 本当の配列サイズ
    
    apparent_size % 8 ? true_size = round_up(apparent_size, 8) : true_size = apparent_size;
    
    auto* mem1 = static_cast<float*>(aligned_alloc(32, sizeof(float)*true_size));
    auto* mem2 = static_cast<float*>(aligned_alloc(32, sizeof(float)*true_size));
    
    for(unsigned i = 0; i<true_size; ++i){
        mem1[i] = 0.f;
        mem2[i] = 0.f;
    } 
    
    for(unsigned i = 0; i<apparent_size; ++i){
        mem1[i] = i;
        mem2[i] = i;
    } 
   
    auto result_avx = dot_avx(true_size, mem1, mem2);
    auto result_nomal = dot(true_size, mem1, mem2);

    std::cout << "simd利用時の結果:" << std::fixed << std::setprecision(20) << result_avx << std::endl;
    std::cout << "simdを利用せず計算した場合の結果:" << std::fixed << std::setprecision(20) << result_nomal << std::endl;
    std::cout << "結果の差 = " << std::fixed << std::setprecision(20) << result_avx - result_nomal << std::endl;
    result_avx == result_nomal ? std::cout << "計算結果が一致" << std::endl : std::cout << "計算結果が違う" << std::endl << std::endl;  

    std::free(mem1);
    std::free(mem2);
    
    return 0;
}

上記コードの実行結果

simd利用時の結果:625598848.00000000000000000000
simdを利用せず計算した場合の結果:625598656.00000000000000000000
結果の差 = 192.00000000000000000000
計算結果が違う

開発環境の備考

種類	名前	バージョン	備考
os	Linux Mint	20	-
コンパイラ	glang++	10	本文中のコードをコンパイルする時オプションに-mavxと-std=c++2aを指定
cpu	Intel© Core™ i7-7700 CPU @ 3.60GHz × 4	-	-
ビルドツール	cmake	3.16.3	-

追記(2020/10/29)

raccyさんが回答内で教えて下さった__カハンの総和アルゴリズム__を用いて
内積関数を作成してみたところ誤差が減少したので追記しておきます。

#include <iostream>
#include <cstdint>
#include <immintrin.h>
#include <iomanip>

inline size_t round_up(const size_t close_to_num, const size_t divisible_by)noexcept{//    divisible_byで割り切れる最もclose_to_numに近い整数を算出
    return divisible_by == 0 ? close_to_num : close_to_num + divisible_by - (close_to_num % divisible_by);
}

float dot_avx(const size_t n, const float* const a, const float* const b)noexcept{//   simdでの内積計算(カハンの加算アルゴリズム使用)
    //  ここからa[i]*b[i]の内容を8つに分けてymmレジスタ(変数sum_avx)格納
    __m256 sum_avx, loss_avx, y_avx, t_avx;
    sum_avx = _mm256_mul_ps(_mm256_load_ps(a), _mm256_load_ps(b));
    loss_avx = _mm256_setzero_ps();
    for(size_t i = 8; i < n; i += 8){
        y_avx = _mm256_sub_ps(_mm256_mul_ps(_mm256_load_ps(a + i), _mm256_load_ps(b + i)), loss_avx);
        t_avx = _mm256_add_ps(sum_avx, y_avx);
        loss_avx = _mm256_sub_ps(_mm256_sub_ps(t_avx, sum_avx), y_avx);
        sum_avx = t_avx;
    }
    
    alignas(32) float array[8];
    _mm256_store_ps(array, sum_avx);
    
    //  ここから通常メモリ(変数array)の内容総和
    float sum, loss, y, t;
    sum = array[0];
    loss = 0.f;
    for(size_t i = 1; i < 8; ++i){
        y = array[i] - loss;
        t = sum + y;
        loss = (t - sum) - y; 
        sum = t; 
    }
    return sum;
}

float dot(const size_t n, const float* const a, const float* const b)noexcept{//   simdを使わず内積計算(カハンの加算アルゴリズム使用)
    float sum, loss, y, t;
    sum = a[0] * b[0];
    loss = 0.f;                            //  情報落ちによって失われてしまう数値を保管するための変数
    for(size_t i = 1; i < n; ++i){
        y = a[i] * b[i] - loss;            //  a[i]*b(i)に前回の変数tの計算で失われてしまった部分を反映                       
		t = sum + y;                       //  sumにyを足している(sumが大きな値でyが小さな値ならば情報落ちにより大きい誤差が発生)                        
		loss = (t - sum) - y;              //  tを計算する時に失われた部分                  
		sum = t;       
    }
    return sum;
}

int main(){
    size_t apparent_size = 1234;// 見かけ上のサイズ
    size_t true_size = (apparent_size % 8 ? round_up(apparent_size, 8) : apparent_size);// 本当の配列サイズ

    auto* mem1 = static_cast<float*>(aligned_alloc(32, sizeof(float)*true_size));
    auto* mem2 = static_cast<float*>(aligned_alloc(32, sizeof(float)*true_size));

    for(unsigned i = 0; i<true_size; ++i){
        mem1[i] = 0.f;
        mem2[i] = 0.f;
    } 

    for(unsigned i = 0; i<apparent_size; ++i){
        mem1[i] = i;
        mem2[i] = i;
    } 

    auto result_avx = dot_avx(true_size, mem1, mem2);
    auto result_nomal = dot(true_size, mem1, mem2);

    std::cout << std::fixed << std::setprecision(20) << std::endl;
    std::cout << "simd利用時の結果:" << result_avx << std::endl;
    std::cout << "simdを利用せず計算した場合の結果:" << result_nomal << std::endl;
    std::cout << "結果の差 = " << result_avx - result_nomal << std::endl;
    std::cout << (result_avx == result_nomal ? "計算結果が一致" : "計算結果が違う") << std::endl; 

    std::free(mem1);
    std::free(mem2);

    return 0;
}

結果

simd利用時の結果:625599104.00000000000000000000
simdを利用せず計算した場合の結果:625599104.00000000000000000000
結果の差 = 0.00000000000000000000
計算結果が一致

正確な内積結果:625599129
本コードの内積結果:625599104
差:25

Answer 1

dot()関数の代わりに下のdot2()関数を試してみてください。

float dot2(const int n, const float *const a, const float *const b) noexcept
{
	float array[8] = {0};
	for (int i = 0; i < n; i += 8)
		for (int j = 0; j < 8; j++)
			array[j] += a[i + j] * b[i + j];
	return array[0] + array[1] + array[2] + array[3] + array[4] + array[5] +
	       array[6] + array[7];
}

実際に実行すればわかるとおり、SIMDを使った場合と同じ結果になります。では、dot()とdot2()の違いは何でしょうか？

そう、足し算する順番が異なるだけで、二つの関数は数学的には同じ計算をしています。しかし、プログラミングの世界では同じ計算とは言えません。それは、途中の計算がfloatであるため、有効桁数からはみ出る部分は近似値計算になるからです。各計算で結果は丸められることになり、その結果には不確かさ生じます。それは足し算の順番が異なることで不確かな分だけ、差が生じる場合が出ると言うことです。dot2はSIMDと同じたし方をしているから、同じになったと言うことです。

なお、正確に整数のみで計算すれば625599129になると思います。floatの有効桁数は十進数で約7桁ですので、62559[89]...という結果であればfloatの精度において計算は正しいと言えるでしょう。

Comments

[退会済みユーザー]

raccyさん
ご回答誠にありがとうございます。

> dot2はSIMDと同じたし方をしているから、同じになったと言うことです。
計算精度が違うだけで計算方法自体は同じだったのですね...
SIMDと同じたし方をすれば精度を揃えることが出来るということは
思いつきもしなかったのでとても助かりました。

的確なご回答誠にありがとうございます。

[raccy]

いえ、精度は全てfloatの単精度浮動少数点数です。精度が異なるのではなく、計算の順序が異なるとか同じであるとかと言うことです。計算の順序が異なると丸め誤差の箇所が変わって、精度によって保証されない部分(6～7桁目以降)は異なっていくという浮動小数点数の特性の話です。

[退会済みユーザー]

raccyさん
お返事誠にありがとうございます。
返信に時間がかかってしまい申し訳ございません。

> 計算の順序で丸め誤差の箇所が変わってしまう

丸め誤差については二進数で循環小数となってしまう部分をある桁数で
切って四捨五入したときに発生する誤差だと私は思っていますが
なぜその誤差が計算の順番で変化するのかがわかりません。

ここがわかれば貴方の御回答の理解に繋がりそうなので教えて頂けたら幸いです。

お手数お掛けしますがお返事頂けましたら幸いです。

[raccy]

floatの精度では、`100000000.0f`に`1.0f`を1000000回足しても`100000000.0f`です。しかし、あらかじめ`1.0f`を1000000回足し合わせて`1000000.0f`にしたものを`100000000.0f`に足すと`101000000.0f`になります。つまり、加算の仕方によって、誤差の蓄積が無視できないほど増幅される場合があると言うことです。実際の結果は下記リンク先を見てください。
https://wandbox.org/permlink/LNxWrZl2kHdyeRoS
浮動小数点数の和についてなるべく誤差を少なくしたいという場合は「カハンの加算アルゴリズム」を使う等の工夫が必要になるでしょう。
https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%8F%E3%83%B3%E3%81%AE%E5%8A%A0%E7%AE%97%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0

[退会済みユーザー]

raccyさん
お返事誠にありがとうございます。

結果が異なる原因がわかったような気がするので解釈の確認をお願いします。

dot関数もdot_avx関数もdot2関数も誤差があるので正確ではないが

dot関数ではどんどん加算していくうちに変数sumの中の数がどんどん大きくなっていき、
有効桁数より小さい数が足されても誤差が発生していた。
------------------------------------------
正確な計算値:625599129
dot関数の計算値:625598656
差:473
------------------------------------------

dot_avx関数ではfloat変数が8つに別れておりfor文が終わった時点では
dot関数のsumにたされた数値より大きくなっていたので
それらを総和するときに発生する誤差がdot関数より小さくなった。
------------------------------------------
正確な計算値:625599129
dot_avx関数の計算値:625598848
差:281
------------------------------------------

浮動小数点型では大きい数に小さい数を足しても誤差が発生するので
貴方の作成したdot2関数ではdot_avx関数同様に総和用の変数を8つに分割して
大きな数値に何度も小さな数値を足していくのではなく
for文でたされた中規模の数値8つを最後に足しているので
dot_avx関数とdot2関数の誤差を合わせることが出来た。

ということですね。

お手数おかけしますがどうかお返事お待ちしております。

[raccy]

だいたいそんなところです。

[退会済みユーザー]

raccyさん。

浮動小数点型では大きい数に小さい数を足しても誤差が発生するということは
全く存じ上げていなかったのでとてもためになりました。

教えていただけとても光栄に感じております。

ご回答誠にありがとうございました。