PA 2.2: Love is Sparse¶
CEGM1000 MUDE: Week 2.2 Due: complete this PA prior to class on Friday, Nov 22, 2024.
Overview¶
This assignment will introduce you to the concept of sparse matrices in Python and how they can be useful to speed up computations and reduce file sizes. To this end, we will be using the scipy.sparse
library.
You will pass this assignment as long as your respository fulfills the following criteria:
- You have completed this notebook and it runs without errors
Reading¶
Keep the scipy.sparse
documentation handy. Some of the work you'll do is based off this blog, so you may find it helpful. In addition, if you don't know what a byte is, you may want to read up on Wikipdia here (not all of it, as long as you recognize that it is a measure of storage space on a computer).The concepts you learn here are applied to the Finite Element Method in this book chapter UPDATE LINK, which you are expected to read during Week 2.2.
Note: you probably skipped over all the links in the paragraph above. While we forgive you for your haste, just remember to revisit some of them if you are struggling to finish the questions below!
import numpy as np
import scipy.sparse as sparse
import matplotlib.pyplot as plt
import timeit
Task 1: Why sparse?¶
Some matrices have a lot of zeros, with such an example given below. When this is the case, the way we store the actual information of the matrix (the non-zero elements) can have a big impact on computation speed and storage demands. Formats which handle this by only storing non-zero elements are called sparse, and have very different internal representations of the data to the matrices you have been familiarized with in previous programming assignments.
Task 1.1:
- Create a function (
create_dense
) which returns a square matrix of arbitrary size. - The function will take as input the size N (such that the matrix is N x N) and one float between 0 and 1, which represents the approximate fraction of the elements of the matrix which are non-zero (it doesn't have to be super accurate).
For now it just return a regular Numpy matrix. To do this, you can use numpy.random.rand to create a random set of values between 0 and 1 and then threshold the entries with a simple boolean operator.
# SOLUTION (a few things will be removed)
def create_dense(size: int, percentage: float) -> np.array:
matrix = np.random.rand(size, size)
matrix[matrix > percentage] = 0
return matrix
####
Now, set up a test to check if you set the function up correctly:
test_size = 1000
test_percentage = 0.1
matrix = create_dense(test_size, test_percentage)
assert np.count_nonzero(matrix) < test_percentage*1.1*test_size**2
Task 1.2: Use array.nbytes to find out how much space a 1000x1000 matrix with 10% non-zero elements takes. Try to explain where this number came from! (Hint: the answer is in the assert statement)
# Solution
my_matrix_size = matrix.nbytes
assert my_matrix_size == 8*test_size**2
Next we will explore how to use scipy.sparse
, and how this reduces the data size of the matrix. The documentation gives us many different types of formats to choose from, so we'll explore two of them: BSR (Block Sparse Row) and CSR (Compressed Sparse Row).
Task 1.3:
Complete the code below to make a CSR and BSR matrix from the matrix
variable.
# SOLUTION
csr_matrix = sparse.csr_array(matrix)
bsr_matrix = sparse.bsr_array(matrix)
Let's compare the new storage requirements and see how much of an improvement we got (it should approach the value used above for test_percentage
, but not reach it exactly):
print(f"CSR matrix size: {csr_matrix.data.size} bytes")
print(f"Compared to the normal matrix, CSR uses this fraction of space: {csr_matrix.data.nbytes/my_matrix_size:0.3f}")
print(f"BSR matrix size: {bsr_matrix.data.size} bytes")
print(f"Compared to the normal matrix, BSR uses this fraction of space: {bsr_matrix.data.nbytes/my_matrix_size:0.3f}")
Task 2: What is love?¶
Let's look into a small example of how sparse matrices can also help improve calculation speeds. We'll study the mysterious case of a massive friend group with a concerning love circle and how we can predict how each person feels.
We know there is a certain pecking order in this group, and neighbours in this order have a love-hate relationship which can be quantified with a simple differential equation:
$$ \begin{pmatrix} \cfrac{dn_1}{dt}\\ \cfrac{dn_2}{dt} \\ \end{pmatrix} = \begin{pmatrix} 0 & 1\\ -1 & 0 \\ \end{pmatrix} \begin{pmatrix} n_1\\ n_2 \\ \end{pmatrix} $$
The state of any given person indicates how much they love the group in general. So in this case, person 2 doesn't like it when person 1 is happy. If we extend this to a four case scenario we'd get the following matrix: $$ \begin{pmatrix} \cfrac{dn_1}{dt}\\ \cfrac{dn_2}{dt}\\ \cfrac{dn_3}{dt}\\ \cfrac{dn_4}{dt}\\ \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 & -1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ 1 & 0 & -1 & 0 \\ \end{pmatrix} \begin{pmatrix} n_1 \\ n_2 \\ n_3 \\ n_4 \\ \end{pmatrix} $$
What happens if we extend it to even more people?
Coincidentally this is very similar to how we use individual elements in the Finite Element Method! We can easily operationalize it using the method ix_
, for which a simple example is provided in the code cell below (this example is generic to illustrate ix_
usage and is not related to the love circle!):
blank = np.zeros(shape=(4, 4))
blueprint = np.array([[0, 0.5],
[1, 0.5]])
for i in range(2):
# First argument will be used for rows
# Second for columns
blank[np.ix_([i*2, i*2 + 1], [1, 2])] = blueprint
print(blank)
Task 2.1:
Generate the matrix relationship
for the differential equation for 1000 people. Use the numpy.ix_
function to make your life easier.
N = 1000
relationship = np.zeros(shape=(N, N))
# SOLUTION
general = np.array([[0, 1], [-1, 0]])
for i in range(N):
relationship[np.ix_([i, (i+1)%N], [i, (i+1)%N])] = general
##########
Finally, we are going to use the forward Euler method to simulate this differential equation for a total of 5 seconds over 1000 iterations. This has already been implemented in the test
method.
Task 2.2:
Find the time it takes to evaluate the relationship using timeit
by entering the function you wish to evaluate as a string. HINT: you have already learned how to convert a matrix into a sparse format, and the function is defined for you. Run the code cell and compare the performances of the different matrix formats. Which one is faster? How much space do they take?
N_ITS = 1000
T = 5 # Seconds
dt = T/N_ITS
def test(rel_matrix):
state = np.zeros(N); state[0] = 1
for i in range(N_ITS):
state = state + rel_matrix @ state * dt
# SOLUTION (only the strings in timeit will be removed + matrix definitions)
csr_matrix = sparse.csr_array(relationship)
bsr_matrix = sparse.bsr_array(relationship)
print(f"Standard: {timeit.timeit('test(relationship)', globals=globals(), number=10)/10:.4f}")
print(f"CSR: {timeit.timeit('test(csr_matrix)', globals=globals(), number=10)/10:.4f}")
print(f"BSR: {timeit.timeit('test(bsr_matrix)', globals=globals(), number=10)/10:.4f}")
########
One final consideration when using sparse matrices is that it can take a long time to generate them from a regular matrix. You can test this out by placing the matrix generation inside or outside the timeit code to compare their performances.
End of notebook.