np.einsum vs np.dot giving different results

Question

Why don't these calculations give identical results?

import numpy as np
M = 1000
N = 500
tab = np.random.random_sample([N,M])
vectors = np.random.random_sample([P,M])
np.einsum('ij,kj->ki',tab,vectors) - np.dot(tab,vectors.T).T

Why is np.einsum('ij,kj->ki',tab,vectors) unequal to np.dot(tab,vectors.T).T?

Note that in terms of run time, np.dot(tab,vectors.T).T is faster than np.einsum('ij,kj->ki',tab,vectors).

Answer 1

It's a precision problem. Lets take a look to the result of np.einsum('ij,kj->ki',tab,vectors) - np.dot(tab,vectors.T).T with smaller dimension

import numpy as np
M = 5
N = 5
P = 2
tab = np.random.random_sample([N,M])

vectors = tab

print np.einsum('ij,kj->ki',tab,vectors) - np.dot(tab,vectors.T).T

>> [[  0.00000000e+00   2.22044605e-16   2.22044605e-16   2.22044605e-16
    0.00000000e+00]
 [  2.22044605e-16   0.00000000e+00   2.22044605e-16   0.00000000e+00
    0.00000000e+00]
 [  2.22044605e-16   2.22044605e-16   0.00000000e+00  -4.44089210e-16
    0.00000000e+00]
 [  2.22044605e-16   0.00000000e+00  -4.44089210e-16   0.00000000e+00
    0.00000000e+00]
 [ -2.22044605e-16   0.00000000e+00   0.00000000e+00   0.00000000e+00
    0.00000000e+00]]

As we can see, it gives a very "small" floats. Let's now do the same thing with int dtype instead of float

import numpy as np
import random as rd
M = 5
N = 5
P = 2
tab = np.array([ rd.randint(-10,10) for i in range(N*M) ]).reshape(N,M)

vectors = tab

print np.einsum('ij,kj->ki',tab,vectors) - np.dot(tab,vectors.T).T

>> [[0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]]

So, what you'r trying to do will never give zeros array for the simple reason that np.einsum has a more precise floating point than np.dot() ( because of the positive sign of the first' result )