Matrix Multiplication

matmul.py

>>> import numpy as np

>>> # Create a random matrix with 3 rows and 4 columns.
>>> matrix_a = np.random.random_sample((3, 4))
>>> matrix_a
array([[ 0.21922347,  0.84313988,  0.41381942,  0.53553901],
       [ 0.35322431,  0.38337327,  0.15964194,  0.30629508],
       [ 0.16188791,  0.55971721,  0.33561351,  0.04709838]])

>>> # Create a matrix from the transpose of A. Note the shapes 
>>> # dimensions.
>>> matrix_b = matrix_a.T
>>> print(matrix_a.shape, matrix_b.shape)
(3, 4) (4, 3)
     
>>> # Next we compute the dot product of our matrices, in this
>>> # case the numpy function for matrix product.
>>> matrix_a.dot(matrix_b)
array([[ 1.21699233,  0.63076825,  0.67151594],
       [ 0.63076825,  0.3910447 ,  0.33976735],
       [ 0.67151594,  0.33976735,  0.45434574]])

>>> # Matrix multiplication is not commutative (AB != BA)
>>> matrix_a.dot(matrix_b) == matrix_b.dot(matrix_a)
False

>>> # We can also do the computation manually (inefficient)
>>> def compute_matrix_product(a, b):
...     result = np.zeros((a.shape[0], b.shape[1]))
...     for i in range(len(a)):
...         for j in range(len(b[0])):
...             for k in range(len(b)):
...                 result[i][j] += a[i][k] * b[k][j]
...     return result
... 
>>> compute_matrix_product(matrix_a, matrix_b)
array([[ 1.21699233,  0.63076825,  0.67151594],
       [ 0.63076825,  0.3910447 ,  0.33976735],
       [ 0.67151594,  0.33976735,  0.45434574]])


>>> # Lets assert that the function does what it should -
>>> # we use np.allclose because the floating points can
>>> # be off by a very very tiny amount.
>>> np.allclose(compute_matrix_product(matrix_a, matrix_b), np.dot(matrix_a, matrix_b)) and \
... np.allclose(compute_matrix_product(matrix_b, matrix_a), np.dot(matrix_b, matrix_a))
True