高精度矢量內積計算方法 (單精度浮點, 超長矢量)

對于單精度浮點類型的超長矢量(超過1億元素)內積計算，累加誤差確實是一個重要問題。以下是幾種減少誤差的方法：

1. Kahan求和算法

這是最常用的補償求和算法，可以有效減少累加誤差：

float kahan_dot_product(const float* a, const float* b, size_t n) {float sum = 0.0f;float c = 0.0f; // 補償變量for (size_t i = 0; i < n; i++) {float product = a[i] * b[i];float y = product - c;float t = sum + y;c = (t - sum) - y;sum = t;}return sum;
}

2. 成對求和 (Pairwise Summation)

將求和過程組織為二叉樹形式，減少大數吃小數的機會：

float pairwise_dot(const float* a, const float* b, size_t start, size_t end) {if (end - start <= 128) { // 基礎塊大小可調整float sum = 0.0f;for (size_t i = start; i < end; i++) {sum += a[i] * b[i];}return sum;}size_t mid = start + (end - start) / 2;float left = pairwise_dot(a, b, start, mid);float right = pairwise_dot(a, b, mid, end);return left + right;
}

3. 多精度累加器

使用雙精度或更高精度累加器：

float high_precision_dot(const float* a, const float* b, size_t n) {double sum = 0.0;for (size_t i = 0; i < n; i++) {sum += (double)a[i] * (double)b[i];}return (float)sum;
}

4. 分塊并行計算

對于超長矢量，可以結合并行計算和誤差控制：

float parallel_dot(const float* a, const float* b, size_t n, int num_threads) {std::vector<double> partial_sums(num_threads, 0.0);#pragma omp parallel for num_threads(num_threads)for (int t = 0; t < num_threads; t++) {size_t start = t * n / num_threads;size_t end = (t + 1) * n / num_threads;double local_sum = 0.0;for (size_t i = start; i < end; i++) {local_sum += (double)a[i] * (double)b[i];}partial_sums[t] = local_sum;}double total = 0.0;for (double s : partial_sums) {total += s;}return (float)total;
}

實際應用建議

1. 對于GPU計算：使用Kahan求和或成對求和，因為雙精度運算在GPU上可能性能較差
2. 對于CPU計算：優先考慮雙精度累加器方法，簡單且有效
3. 混合精度：可以考慮使用融合乘加(FMA)指令結合補償算法
4. 內存訪問優化：確保數據連續訪問，減少緩存未命中

誤差分析

對于1億個單精度數的內積：

• 普通累加：相對誤差可能達到1e-4量級
• Kahan求和：誤差可降低到1e-8量級
• 雙精度累加：誤差可降低到1e-16量級

選擇哪種方法取決于您的精度要求和計算環境限制。

---

測試代碼

#include <cmath>
#include <iostream>
#include <iomanip>
#include <omp.h>
#include <cstdlib>
#include <immintrin.h>using namespace std;//編譯: g++ -DAVX2 -DAVX512F -std=c++17 -Ofast -march=native -fopenmp test_dot.cpp/********************************************************<x,y>內積計算用double雙精度做累加類型,保證數值穩定
**********************************************************/template<typename F,int P=0>
F dot(int n, const F *x, const F *y)
{if constexpr (P==0){double s_time=omp_get_wtime();//累加用單精度F s=0;for(int i=0; i<n; i++){s+=x[i]*y[i];}double e_time=omp_get_wtime();cout<<"P="<<P<<": time used "<<e_time-s_time<<endl;return s;}else if constexpr (P==1){double s_time=omp_get_wtime();//累加用雙精度,乘法用單精度double s=0;for(int i=0; i<n; i++){s+=x[i]*y[i];}double e_time=omp_get_wtime();cout<<"P="<<P<<": time used "<<e_time-s_time<<endl;return s;}else if constexpr (P==2){double s_time=omp_get_wtime();//累加用雙精度,乘法用雙精度double s=0;for(int i=0; i<n; i++){double a=x[i];double b=y[i];s+=a*b;}double e_time=omp_get_wtime();cout<<"P="<<P<<": time used "<<e_time-s_time<<endl;return s;}
#ifdef AVX2 else if constexpr(P==3){static_assert(is_same_v<F,float>);/********************************************************************OpenMP多線程,AVX2計算<x,x>,<x,y>累加和乘法都用雙精度	*********************************************************************/double s_time=omp_get_wtime();__m256d xx_sum,xy_sum;xx_sum=xy_sum=_mm256_setzero_pd();#pragma omp parallel{__m256d sum_xx=_mm256_setzero_pd();__m256d sum_xy=_mm256_setzero_pd();#pragma omp for nowaitfor(int i=0; i<n; i+=8){__m256 x8=_mm256_loadu_ps(x+i);__m256 y8=_mm256_loadu_ps(y+i);__m128 lo,hi;__m256d t1,t2,t3,t4;lo=_mm256_extractf128_ps(x8,0);hi=_mm256_extractf128_ps(x8,1);t1=_mm256_cvtps_pd(lo);t2=_mm256_cvtps_pd(hi);lo=_mm256_extractf128_ps(y8,0);hi=_mm256_extractf128_ps(y8,1);t3=_mm256_cvtps_pd(lo);t4=_mm256_cvtps_pd(hi);
#if 0				sum_xx=_mm256_add_pd(sum_xx,_mm256_mul_pd(t1,t1));sum_xx=_mm256_add_pd(sum_xx,_mm256_mul_pd(t2,t2));sum_xy=_mm256_add_pd(sum_xy,_mm256_mul_pd(t1,t3));sum_xy=_mm256_add_pd(sum_xy,_mm256_mul_pd(t2,t4));
#else/**********************************FMA***********************************/sum_xx=_mm256_fmadd_pd(t1,t1,sum_xx);sum_xx=_mm256_fmadd_pd(t2,t2,sum_xx);sum_xy=_mm256_fmadd_pd(t1,t3,sum_xy);sum_xy=_mm256_fmadd_pd(t2,t4,sum_xy);
#endif}#pragma omp critical{xx_sum=_mm256_add_pd(xx_sum,sum_xx);xy_sum=_mm256_add_pd(xy_sum,sum_xy);}}double tmp[4];_mm256_storeu_pd(tmp,xy_sum);double xy=tmp[0]+tmp[1]+tmp[2]+tmp[3];_mm256_storeu_pd(tmp,xx_sum);double xx=tmp[0]+tmp[1]+tmp[2]+tmp[3];for(int i=n&~7; i<n; i++){double a=x[i];double b=y[i];xx+=a*a;xy+=a*b;}double e_time=omp_get_wtime();cout<<"P="<<P<<": time used "<<e_time-s_time<<endl;cout<<"xx="<<xx<<endl;return xy;}//P==3
#endif 	
#ifdef AVX512F 	else if constexpr(P==4){static_assert(is_same_v<F,float>);/********************************************************************OpenMP多線程,AVX512F計算<x,x>,<x,y>累加和乘法都用雙精度	*********************************************************************/double s_time=omp_get_wtime();__m512d xx_sum=_mm512_setzero_pd();__m512d xy_sum=_mm512_setzero_pd();#pragma omp parallel{__m512d sum_xx=_mm512_setzero_pd();__m512d sum_xy=_mm512_setzero_pd();#pragma omp for nowaitfor(int i=0; i<n; i+=16){__m512 x16=_mm512_loadu_ps(x+i);__m512 y16=_mm512_loadu_ps(y+i);__m256 lo,hi;__m512d t1,t2,t3,t4;lo=_mm512_extractf32x8_ps(x16,0);hi=_mm512_extractf32x8_ps(x16,1);t1=_mm512_cvtps_pd(lo);t2=_mm512_cvtps_pd(hi);lo=_mm512_extractf32x8_ps(y16,0);hi=_mm512_extractf32x8_ps(y16,1);t3=_mm512_cvtps_pd(lo);t4=_mm512_cvtps_pd(hi);
#if 0				sum_xx=_mm512_add_pd(sum_xx,_mm512_mul_pd(t1,t1));sum_xx=_mm512_add_pd(sum_xx,_mm512_mul_pd(t2,t2));sum_xy=_mm512_add_pd(sum_xy,_mm512_mul_pd(t1,t3));sum_xy=_mm512_add_pd(sum_xy,_mm512_mul_pd(t2,t4));
#else/***********************************FMA************************************/sum_xx=_mm512_fmadd_pd(t1,t1,sum_xx);sum_xx=_mm512_fmadd_pd(t2,t2,sum_xx);sum_xy=_mm512_fmadd_pd(t1,t3,sum_xy);sum_xy=_mm512_fmadd_pd(t2,t4,sum_xy);
#endif}#pragma omp critical{xx_sum=_mm512_add_pd(xx_sum,sum_xx);xy_sum=_mm512_add_pd(xy_sum,sum_xy);}}double tmp[8];_mm512_storeu_pd(tmp,xy_sum);double xy=tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];_mm512_storeu_pd(tmp,xx_sum);double xx=tmp[0]+tmp[1]+tmp[2]+tmp[3]+tmp[4]+tmp[5]+tmp[6]+tmp[7];for(int i=n&~15; i<n; i++){double a=x[i];double b=y[i];xx+=a*a;xy+=a*b;}double e_time=omp_get_wtime();cout<<"P="<<P<<": time used "<<e_time-s_time<<endl;cout<<"xx="<<xx<<endl;return xy;}
#endif return(0);}const int N=50000000;void test()
{using FLOAT=float;FLOAT *x=new FLOAT[N];FLOAT *y=new FLOAT[N];for(int i=0; i<N; i++){FLOAT t=0.001*sqrtf(FLOAT(i));FLOAT s=sqrt(sqrt(FLOAT(i)));x[i]=(rand()<RAND_MAX/2)?t:-t;y[i]=(rand()<RAND_MAX/2)?s:-s;}cout<<setprecision(15)<<endl;cout<<(double)dot<FLOAT,0>(N,x,y)<<endl;cout<<(double)dot<FLOAT,1>(N,x,y)<<endl;cout<<(double)dot<FLOAT,2>(N,x,y)<<endl;#ifdef AVX2	cout<<(double)dot<FLOAT,3>(N,x,y)<<endl;
#endif #ifdef AVX512F	cout<<(double)dot<FLOAT,4>(N,x,y)<<endl;
#endif }int main(int argc, char **argv)
{test();return(0);
}

---

資料

Intel /intrinsics-guide