[SciPy-User] Matrix by matrix vs. matrix by vector multiplications are not the same

[Alexander Kalinin]

Hello,

Look at this code:


import numpy as np


N = 500

A = np.random.rand(N, N)

B = np.random.rand(N, N)


c1 = np.dot(A, B)[:, 0]

c2 = np.dot(A, B[:, 0])


print np.linalg.norm(c1 - c2)

The output is:
3.92795839192e-13

For me it is little bit strange that the results are not the same. Does it
mean that matrix by matrix and matrix by vector multiplications algorithms
are different?

Sincerely,
Alexander

[Pauli Virtanen]

05.05.2013 14:18, Alexander Kalinin kirjoitti:
[clip]

Post by Alexander Kalinin
For me it is little bit strange that the results are not the same. Does
it mean that matrix by matrix and matrix by vector multiplications
algorithms are different?

Correct, they are different. BLAS libraries typically try to optimize
for speed and the most appropriate algorithms for each case are
different. The speed increase gained by this can be quite significant.

The underlying issue is that in floating point, addition and
multiplication are not associative operations, so that mathematically
equivalent algorithms may produce somewhat different results due to
different accumulation of rounding error.

[Alexander Kalinin]

Thank you for explanation.

Alexander.

Post by Pauli Virtanen
[clip]

Post by Alexander Kalinin
For me it is little bit strange that the results are not the same. Does
it mean that matrix by matrix and matrix by vector multiplications
algorithms are different?

Correct, they are different. BLAS libraries typically try to optimize
for speed and the most appropriate algorithms for each case are
different. The speed increase gained by this can be quite significant.
The underlying issue is that in floating point, addition and
multiplication are not associative operations, so that mathematically
equivalent algorithms may produce somewhat different results due to
different accumulation of rounding error.
--
Pauli Virtanen
_______________________________________________
SciPy-User mailing list
http://mail.scipy.org/mailman/listinfo/scipy-user