0

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
WhatsupAndThanks
  • 35
  • 6
  • I only get your `Naive Dot` (`x@y`) value. I'm limited, of course, by the "accuracy" of the posted `x,y` printed values. With `math.fsum` I get a slightly different value, `0.733049432917035`. The greater decimal places of `mpmath` doesn't help. – hpaulj Jun 25 '23 at 16:54
  • The problem isn't so much with decimal places, as it is with the wide range of product values, like `128297332517598.19` and `-43.18555935175856`, call cancelling each other out to end up with a sum less than one. – hpaulj Jun 25 '23 at 17:43
  • There's a `mpmath.fsum`, though the result is basically the same as using `math.fsum`. – hpaulj Jun 25 '23 at 17:46
  • We can't test your coding of the "accurate" product, since all you provide is a pdf link to an article. – hpaulj Jun 25 '23 at 19:20
  • Thanks for testing it. I updated the post with the accurate product code. And Yes, there is canceling because it is a ill-conditioned dot product. The accurate product algorithm is suppose to handle the ill-condition easily. – WhatsupAndThanks Jun 25 '23 at 23:16
  • What do you think `mp.mpf(str(x[i]))` is doing? – Tim Roberts Jun 25 '23 at 23:31
  • my understanding is `mp.mpf(str(x[i]))` turns the numerical values to string representation for mpf to read, whereas `mp.mpf(x[i])` passes the python float values to mpmath directly. – WhatsupAndThanks Jun 25 '23 at 23:37
  • 1
    @TimRoberts, I tested his numbers as array, as list, and as lists of strings (quoting the cooynpaste). Also with and without `mpmath`. The only thing that makes a difference is the summation order. `fsum` appears to sum the high magnitude products separately. – hpaulj Jun 25 '23 at 23:45
  • I'm not going to import `pyfma` just to replicate your results. Do `mpmath` objects survive a passage through `pyfma`? – hpaulj Jun 26 '23 at 00:01
  • You can try changing pyfma.fma into y=a*b-x but that will degrade the accuracy of the resulting dot product – WhatsupAndThanks Jun 26 '23 at 00:25

1 Answers1

0

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
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • I can reproduce your result. I think the difference between mine (0.7289845395392482) comes from the fact that your used y = a*b -x, which is 0. y is the correction term that should be obtained by fma. So your answer here considers the correction from the summation, which is confirmed by fsum in your answer. So at least your answer helps to confirm that the summation part of my accurate product is correct. – WhatsupAndThanks Jun 26 '23 at 16:41