comparing Accurrate Dot Product with mpmath

Question

This post is not a repost, which is related solely to mpmath. But this is mainly related to the comparison between different dot product algorithms.

I have coded an accurate dot product in python in order to compute ill-conditioned dot products.

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])      

To check the validity of the algorithm, I compare it with a naive dot product, and the exact answer (DotExact) which is computed by mpmath library.

The problem is that, depending on how the DotExact is coded, it can agree with the accurate dot product until ~16 digits, or disagrees with accurate dot product and perform the same as the naive case. The Dot Exact is computed as

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])           # 1ST  this agrees with Dot2K                                                                                                                                                                            
#        sum +=  mp.mpf(str(x[i])) * mp.mpf(str(y[i]))  # 2ND this does NOT agree with 

    return sum
                                                                                                 
                                             
                                                                                                                                                                                                   

We expect the naive dot product to be wrong, but does anyone here know why the accurate product algorithm only agrees with the Exact summation performed by 1st, but not 2nd? For your info, I also provide here the results of the dot products based on DotExact_1st

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 

---

Edit1:

Here is the algorithm for the accurate dot product.

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

Answer 1

Your sample arrays (as copy-n-pasted)

In [314]: x = np.array([32888447.473062776 ,-254.18174014817862 ,-0.027520952868535176 ,0.0 ,-63157.1445
     ...: 9198614 ,6.547216534966907 ,0.0 ,-10637698.845199142])
     ...: 
     ...: y = np.array([3900984.764412485 ,603797189919.5646 ,-116303296140.4418 ,0.0 ,-202.098063597826
     ...: 9 ,-6.596018188968733 ,0.0 ,-2366458.6427612416])

The straight forward matrix product:

In [315]: x@y
Out[315]: 0.73046875

Your "accurate", with (a*b-x):

In [316]: Dot2s_f64(x,y)
Out[316]: 0.733049432917035

I get the same thing using the math.fsum, which handles this sort of sum better:

In [317]: import math    
In [318]: math.fsum([i*j for i,j in zip(x,y)])
Out[318]: 0.733049432917035

Try it with mpmath objects:

In [319]: import mpmath as mp

Creating mpmath object arrays (note the object dtype):

In [320]: xm = np.array([mp.mpf(i) for i in x]); xm
Out[320]: 
array([mpf('32888447.473062776'), mpf('-254.18174014817862'),
       mpf('-0.027520952868535176'), mpf('0.0'),
       mpf('-63157.144591986143'), mpf('6.5472165349669069'), mpf('0.0'),
       mpf('-10637698.845199142')], dtype=object)

In [321]: ym = np.array([mp.mpf(i) for i in y]);

matmul can work with objects (just more slowly):

In [322]: xm@ym
Out[322]: mpf('0.73046875')

Again, fsum does just as well:

In [323]: mp.fsum([i*j for i,j in zip(xm,ym)])
Out[323]: mpf('0.73304943291703495')

In [324]: Dot2s_f64(xm,ym)
Out[324]: mpf('0.73304943291703495')

magnitude partition

Make a list of the products:

In [325]: alist = [i*j for i,j in zip(x,y)]

In [326]: alist
Out[326]: 
[128297332517598.19,
 -153474220430335.22,
 3200777531.536388,
 0.0,
 12763936.624408364,
 -43.18555935175856,
 0.0,
 25173674371312.79]

Sum them in two parts - the large magnitude ones and small:

In [327]: arr = np.array(alist)    
In [328]: arr[np.abs(arr)>1e10].sum()+arr[np.abs(arr)<=1e10].sum()
Out[328]: 0.7330493927001953