[FIXED] numpy dot product with floating point

Issue

I was trying to do a dot product using numpy.dot() function. A MWE is given below.

import numpy as np
a = np.array([1.24912871, 2.45123219, 6.87478761, 8.24982742, 3.31231561, 9.42387436])
b = np.array([2.24912871, 5.45123451, 9.87487986, 7.24982743, 1.31231215, 8.42387427])
ab=np.dot(a,b)
print(ab)

a0, a1 = a[:3], a[3:]
b0, b1 = b[:3], b[3:]

ab_split = np.dot(a0,b0) + np.dot(a1,b1)
print(ab_split)

Ideally, ab and ab_split should give the exact same result. However, when I run this code I get the output as

227.6015443489002
227.60154434890018

Now, the difference between the two is not substantial. But, this dot product is required for an iterative algorithm I am running by splitting data across multiple cores using using mpi4py. You can visualize as a0 and b0 to be in processor-0 and a1 and b1 to be in processor-1. Because of this behaviour, the number of iterations in both the cases are different and sometimes it almost doubles.

What I want to know is that if there is any way such that ab_split can be made to achieve the exact same result as ab?

Solution

You are running into the non-associativity of floating point arithmetic. So, ordinarily, no, you can not get two computations with the same result in exact arithmetic to get the same result in floating point arithmetic.

If your inner product calculation is a minor part of the work, you could contemplate doing it in higher precision. Python has a package for that. There is in fact a result that using about 4000 bits for each number is enough to make computer arithmetic exact.

Of course, using non-standard floating point numbers will lower your efficiency.

Answered By - Victor Eijkhout