精确点积算法与mpmath两种实现的一致性问题探究

本文聚焦不同点积算法的对比，在Python中实现了一种精确点积算法用于计算病态点积，测试数组如下：

import numpy as np
x = np.array([32888447.473062776 ,-254.18174014817862 ,-0.027520952868535176 ,0.0 ,-63157.14459198614 ,6.547216534966907 ,0.0 ,-10637698.845199142])                                                                                                      
y = np.array([3900984.764412485 ,603797189919.5646 ,-116303296140.4418 ,0.0 ,-202.0980635978269 ,-6.596018188968733 ,0.0 ,-2366458.6427612416])      

为验证算法有效性，将其与朴素点积、mpmath库实现的DotExact对比，DotExact有两种实现方式：

import mpmath as mp
mp.dps = 64  
def DotExact(x,y):
    n=np.int64(x.shape[0])                                                                                                                                                                                                                                
    sum = mp.mpf(0.)

    for i in range(n):
         sum +=  mp.mpf(x[i]) * mp.mpf(y[i])           # 第一种实现，与精确点积一致
#        sum +=  mp.mpf(str(x[i])) * mp.mpf(str(y[i]))  # 第二种实现，与精确点积不一致

    return sum

测试结果显示：

Condition number    =  19.765763843548616                                                                                                                                                                                                                 
Dot Exact by Mpmath =  0.7289845395392482978061485937582                                                                                                                                                                                                  
Naive Dot           =  0.73046875          agree digits at -2.6912228296231142978105880502362                                                                                                                                                             
Accurate Dot        =  0.7289845395392482  agree digits at -16.13278356042518496483357616017 

精确点积算法的实现代码如下：

import pyfma
def TwoProdFMA(a,b):                                                                                                                                                  
  x = a*b                                                                                                                                                             
  y = pyfma.fma(a, b, -x)                                                                                                                                             
#  y = 0                                                                                                                                                                
  return x,y                                                                                                                                                            


def TwoSumP(a,b):                                                                                                                                                     
    pi = (a + b);                                                                                                                                                     
    z  = (pi - a);                                                                                                                                                      
    q  = (a - (pi-z) + (b-z));                                                                                                                                         

    return pi, q                                                                                                                                                        


def Dot2s_f64(x,y):                                                                                                                                                   

  # accurate dot product                                                                                                                                              
  N = x.shape[0]                                                                                                                                                        
  p,s = TwoProdFMA(x[0],y[0]);                                                                                                                                         

  for i in range(1,N):                                                                                                                                                  
    h, r = TwoProdFMA(x[i],y[i]);                                                                                                                                      
    p, q = TwoSumP(p,h);                                                                                                                                               
    s   += (q+r);                                                                                                                                                       

  return p+s

核心疑问 #

为何精确点积算法仅与第一种DotExact实现结果一致，与第二种实现结果不一致？

原因解析 #

两种实现的核心差异在于numpy float64数值转换为mpmath高精度浮点数的方式不同：

• 第一种实现mp.mpf(x[i])：直接读取numpy数组中存储的二进制float64原始值，将其转换为mpmath的高精度浮点数。这和你实现的精确点积算法处理的是同一个底层数值——也就是那些显示的十进制数在float64精度下实际存储的二进制表示。
• 第二种实现mp.mpf(str(x[i]))：先将float64数值转换为十进制字符串，再让mpmath解析这个字符串。但float64的字符串表示是其十进制近似值（并非精确的底层二进制值），比如32888447.473062776这个字符串对应的精确值，和它在float64中实际存储的二进制值存在细微差异。mpmath解析该字符串后得到的高精度数，已经偏离了原始float64的真实存储值，后续的乘积与求和计算自然会和你的精确点积结果产生不一致。

简言之：你的精确点积算法基于float64的真实二进制值做精确计算，第一种DotExact复用了同样的底层值，第二种则使用了十进制字符串的近似值，最终结果自然不同。

内容的提问来源于stack exchange，提问作者WhatsupAndThanks