In this tutorial we’ll look at some simple matrix multiplication tools in Python. The arrays are constructed using lists within lists, and the multiplication itself is implemented using for loops. The program works well for small arrays, but really starts bogging down at sizes 500×500 and above. Anyone doing serious work will want to check out NumPy, which has matrix operations built-in and runs several thousand times faster. In Part 2 I’ll compare NumPy’s performance to other optimized methods.
Let’s feast our eyes on the function definitions we’ll be using in matrix.py:
import random from time import * import cProfile def zero(m,n): # Create zero matrix new_matrix = [[0 for row in range(n)] for col in range(m)] return new_matrix def rand(m,n): # Create random matrix new_matrix = [[random.random() for row in range(n)] for col in range(m)] return new_matrix def show(matrix): # Print out matrix for col in matrix: print col def mult(matrix1,matrix2): # Matrix multiplication if len(matrix1[0]) != len(matrix2): # Check matrix dimensions print 'Matrices must be m*n and n*p to multiply!' else: # Multiply if correct dimensions new_matrix = zero(len(matrix1),len(matrix2[0])) for i in range(len(matrix1)): for j in range(len(matrix2[0])): for k in range(len(matrix2)): new_matrix[i][j] += matrix1[i][k]*matrix2[k][j] return new_matrix def time_mult(matrix1,matrix2): # Clock the time matrix multiplication takes start = clock() new_matrix = mult(matrix1,matrix2) end = clock() print 'Multiplication took ',end-start,' seconds' def profile_mult(matrix1,matrix2): # A more detailed timing with process information # Arguments must be strings for this function # eg. profile_mult('a','b') cProfile.run('matrix.mult(' + matrix1 + ',' + matrix2 + ')')
Load the script into an interactive Python session:
>>> import matrix >>> a = matrix.zero(1,5) >>> a [[0, 0, 0, 0, 0]] >>> b = matrix.zero(3,5) >>> b [[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]] >>> matrix.show(b) # This function makes arrays more readable. [0, 0, 0, 0, 0] [0, 0, 0, 0, 0] [0, 0, 0, 0, 0]
Try multiplying zero matrices just to check the resulting dimensions:
>>> a = matrix.zero(2,3) >>> b = matrix.zero(4,5) >>> c = matrix.mult(a,b) Matrices must be m*n and n*p to multiply! >>> b = matrix.zero(3,5) >>> c = matrix.mult(a,b) >>> matrix.show(c) # Note that we end up with an n by p matrix [0, 0, 0, 0, 0] # In this case 2x5 [0, 0, 0, 0, 0]
Now do some calculations:
>>> a = [[2,3]] # 1x2 matrix >>> b = [[1],[5]] # 2x1 matrix >>> c = matrix.mult(a,b) >>> matrix.show(c) # [1x2] * [2x1] --> [1x1] matrix [17] >>> a = matrix.rand(3,3) >>> b = matrix.rand(3,3) >>> c = matrix.mult(a,b) >>> matrix.show(c) [1.5756730190235944, 0.70870079587302315, 0.94109431416086442] [0.48198264627485315, 0.32042216800153944, 0.36885894389577889] [1.2931982570350344, 0.85088650318922798, 1.1630298461072877]
We can time the multiplication with the last two functions. Try increasing the matrix sizes and see what happens:
>>> a = matrix.rand(3,3) >>> b = matrix.rand(3,3) >>> matrix.time_mult(a,b) Multiplication took 0.0 seconds # 0.o >>> a = matrix.rand(300,300) # Each matrix has 90000 elements now >>> b = matrix.rand(300,300) >>> matrix.time_mult(a,b) Multiplication took 10.35 seconds # Many more calculations happening >>> a = matrix.rand(600,600) >>> b = matrix.rand(600,600) >>> matrix.time_mult(a,b) Multiplication took 91.05 seconds # Not worth going much larger
Take a look at the matrix multiplication loop, mult(matrix1,matrix2). It’s iterating over i, j, and k — each of which is a column or row of a 600×600 matrix. 600^3 = 216 million cycles through Python for loops. That’s O(n^3) complexity, which is quite the pain for large arrays. Take a look at the results of the profiler:
>>> matrix.profile_mult('a','b') 721811 function calls in 91.598 CPU seconds Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.005 0.005 91.598 91.598 <string>:1(<module>) 1 89.996 89.996 91.593 91.593 matrix.py:20(mult) 1 0.080 0.080 0.088 0.088 matrix.py:5(zero) 360605 0.049 0.000 0.049 0.000 {len} 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 361202 1.468 0.000 1.468 0.000 {range}
3 Comments
Hi. Can you please explain this part of the code in detail:
else:
# Multiply if correct dimensions
new_matrix = zero(len(matrix1),len(matrix2[0]))
for i in range(len(matrix1)):
for j in range(len(matrix2[0])):
for k in range(len(matrix2)):
new_matrix[i][j] += matrix1[i][k]*matrix2[k][j]
return show(new_matrix)
Thank you very much.
A new matrix is created. For multiplication of matrices with dimension n*m and m*p, the new matrix will be dimension n*p:
These three lines set up a brute-force loop through all the matrix elements. Wikipedia page has decent examples of how this calculation is performed
The next line performs the calculation for each new element. It’s tough to explain without lots of notation, but this page does a good job of simple examples. You’ll see that elements in the new matrix are a summation of multiplication of corresponding elements in columns of the first matrix, and rows of the second matrix. That’s a really awful sentence, but I can’t think of a simpler way to explain it. The += lets each of these multiplications accumulate on the proper element.
The last line returns the new matrix, so it can be printed, stored as a variable, or whatever you like!
Thanks to author a blog, I actually saved a lot of time! I would go again!