Np.einsum vs np.dot gives different results
Why don't these calculations give the same 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 np.einsum('ij,kj->ki',tab,vectors)
doesn't it match np.dot(tab,vectors.T).T
?
Note that in terms of runtime it is np.dot(tab,vectors.T).T
faster than np.einsum('ij,kj->ki',tab,vectors)
.
source to share
This is a precision issue. Let's see the result np.einsum('ij,kj->ki',tab,vectors) - np.dot(tab,vectors.T).T
with a smaller size
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 very "small" floats. Now let's do the same with int
dtype instead offloat
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 are trying to do will never give an array of zeros for the simple reason that it np.einsum
has a more accurate floating point than np.dot()
(due to the positive sign of the first result)
source to share
The results are the same up to some numerical precision. The differences look like
[ 5.68434189e-14, 0.00000000e+00, 8.52651283e-14, ...,
8.52651283e-14, 0.00000000e+00, -5.68434189e-14],
[ -8.52651283e-14, 0.00000000e+00, -5.68434189e-14, ...,
0.00000000e+00, 5.68434189e-14, -8.52651283e-14],
[ 1.42108547e-13, 5.68434189e-14, 0.00000000e+00, ...,
1.13686838e-13, -5.68434189e-14, 1.13686838e-13]])
As @wflynny points out in the comments, the best way to do this test on arrays a
and b
is
np.allclose(a, b)
Potentially faster method:
from numpy.core.umath_tests import matrix_multiply
matrix_multiply(tab, vectors.T).T
source to share
@YXD
matrix_multiply is very slow!
With the following code:
- sol1 --- 29.6340000629 seconds ---
- sol2 --- 3.78200006485 seconds ---
- sol3 --- 4.25900006294 seconds ---
- sol4 --- 68.1049997807 seconds ---
- sol5 --- 4.06699991226 seconds ---
Code:
import time
import numpy as np
from sklearn.utils.extmath import fast_dot
from numpy.core.umath_tests import matrix_multiply
a =1.1
b = 0.1
d = 2
M = 1000
N = 500
P = 100
tab = np.random.random_sample([N,M])*1000
vectors = np.random.random_sample([P,M])*100
coef = np.random.random_sample([1,P])
#sol1
start_time = time.time()
for i in range(1000):
res1 = np.dot(coef,(a + b*np.einsum('ij,kj->ki',tab,vectors))**d)
print("--- %s seconds ---" % (time.time() - start_time))
#sol2
start_time = time.time()
for i in range(1000):
res2 = np.dot(coef,(a + b*np.dot(tab,vectors.T).T)**d)
print("--- %s seconds ---" % (time.time() - start_time))
print(np.allclose(np.einsum('ij,kj->ki',tab,vectors), np.dot(tab,vectors.T).T))
#sol3
start_time = time.time()
for i in range(1000):
res2 = fast_dot(coef,(a + b*fast_dot(tab,vectors.T).T)**d)
print("--- %s seconds ---" % (time.time() - start_time))
print(np.allclose(np.einsum('ij,kj->ki',tab,vectors), fast_dot(tab,vectors.T).T))
#sol4
start_time = time.time()
for i in range(1000):
res2 = matrix_multiply(coef,(a + b*matrix_multiply(tab,vectors.T).T)**d)
print("--- %s seconds ---" % (time.time() - start_time))
# error ??? print(np.allclose(matrix_multiply('ij,kj->ki',tab,vectors), matrix_multiply(tab, vectors.T).T))
#sol5
start_time = time.time()
for i in range(1000):
res2 = np.dot(coef,(a + b*np.tensordot(tab,vectors.T, axes=([-1],[0])).T)**d)
print("--- %s seconds ---" % (time.time() - start_time))
print(np.allclose(np.einsum('ij,kj->ki',tab,vectors), np.tensordot(tab,vectors.T, axes=([-1],[0])).T))
source to share