Hacking Matrices And Gauss

I was invited to teach a lesson at a not-so-nearby high school and was told exactly the topic the class would be learning that day: matrices (yay!), specifically, solutions of linear systems by Gaussian elimation (awww...).

My heart sank because I can remember a hot autumn day (or did it last a week?) decades ago, when my brother and I were sitting in a business math class at Merrimack College being led through this exact topic by a teacher-in-training. It ranked among the most boring episodes in my educational experience. And pointless, too. Why were we going through all this? Because we had matrices to solve, I guess.

More recently I've posted about the more interesting, visual, interactive applications of matrices in 3D graphics:

So I worked through a few problems in Gaussian elimination to refresh my memory and predictably enough, I thought about automating this boring stuff with Python. Sure, you can just use Numpy's Linear Algebra module to solve them. Let's solve this example (from Wikipedia):

2x + y - z = 8

-3x - y + 2z = -11

-2x + y + 2z = -3

You have to put the coefficients and constants into matrices like this and type a few words to get the solution:

That means x = 2, y = 3 and z = -1. Check it by plugging those values into the original equations. 

But maybe using modern-day solvers is cheating. This is math class, after all. Best to learn to make a solver yourself. I started by dividing the whole first row by the element in the first column. To use a list's values and their place in the list (their index), Python has a handy function called "enumerate." For example, declare a list like this and print the index and value of each item:

Using enumerate, it was easy to divide all the values in the first row (index 0) by the value of the top left entry.

    #get top left entry to be 1

    #by dividing row by A[0][0]

    divisor = A[0][0]

    for i,v in enumerate(A[0]):

        A[0][i] = v / divisor

The next step was making the other elements in that column 0 by adding the negative of those values. Like in our example, the first item in the second row is -3. So to everything in that row we add 3 times the corresponding entry in the first row. Do the same thing for the third row. So having the index and the value of each entry in the matrix is crucial.

    #add row 2 and 3 to additive inverse to make value 0

    for i in [1,2]:

        addinv = -1*A[i][0]

        for ind,value in enumerate(A[i]):

            A[i][ind] = value + addinv*A[0][ind]

Then you do the same thing to the second row, second column and make the other entries 0, and the same thing for the third row, third column, etc. I could smell a loop. You just need to have a list of rows and take out the one that's on the diagonal:

        rowList = [x for x in range(len(A))]
        rowList.remove(j)

I also needed to make sure we weren't dividing by zero. The whole function ended up being less than 20 lines:

But does it work? Let's use our Wikipedia example. I entered the matrix in augmented form, meaning it was a 3x4 matrix:

My function returned the entire matrix, but if you look at the last values in each row, those are the solutions for x, y and z.

The great thing is this function even works for bigger matrices. If you had a system of 4 unknowns, you'd need 4 equations like this:

4w + x + 2y + 3z = 6

3w - 6x + 5y + 4z = 1

9w + 7x - 7y + 8z = 4

w + x + y - z = 1

That would be cruel to make a student do by hand, but using our function, the worst part is just typing it all in:

The last numbers in each row are the solutions for w, x, y and z. Checked it with numpy:

>>> 

[[-1.01875]

 [ 1.8375 ]

 [ 1.75625]

 [ 1.575  ]]

One of the questions the students asked was to find the equation of a parabola given 3 points like this:

We know it's a parabola, so the solution is in the form

y = ax² + bx + c

and we know 3 pairs of x and y coordinates, so it's a system of equations to solve for a, b and c:

(2) = a(-2)² + b(-2) + c

(-1) = a(1)² + b(1) + c

(-6) = a(2)² + b(2) + c

or

4a - 2b + c = 2

a + b + c = -1 

4a + 2b + c = -6

No problem solving this with gauss():

The solutions are a = -1, b = -2 and c = 2. So the equation of the parabola is 

y = -x² - 2x + 2

There are plenty of computer-based and online tools to solve linear systems using matrices. We should be teaching our math students how to make the tools themselves.

Automata to Pyramids

Was flipping through Steven Wolfram's A New Kind of Science and was struck by the picture of the steps in the evolution of this 2D cellular automaton:

If every step became a level in a 3D shape, it'd look like this:



I'm motivated to try this out in Python. The first challenge is to create an array. It's easy enough to type in an 11x11 array for step 0 manually:

     surface = [[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, 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, 0, 0, 0, 0, 0],
               [0, 0, 0, 0, 0, 1, 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, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

But we eventually want a bigger array. It should be easy to create, right? Here's the code for it:

tableSize = 11
surface = tableSize*[tableSize*[0]]

Which gives us a 11x11 array of all zeroes.



The code for making the middle one a 1 is this:

surface[int(tableSize/2)][int(tableSize/2)] = 1

Strangely, that produces this!

Numpy to the Rescue

I have no idea why this glitch is happening, so let's switch to Numpy, Python's professional-strength array-wrangling package. The code for making the array is easier, and it works:

    tableSize = 101
    surface = numpy.zeros((tableSize,tableSize))
    surface[int(tableSize/2)][int(tableSize/2)] = 1



Now I can make a monster array just by changing the tableSize variable.

Automating the Pattern

Now the question is how to produce the pattern in the first picture. For each cell, the program has to check its four neighbors (left, right, up, down) and see whether exactly 1 or 4 of the neighbors is a 1. If so (or if the cell itself is a 1), the cell is a 1 in the next level. This function will count the number of 1's among a cell's four neighbors:

def nbs(listA,r,c):
    '''counts the 1's in the 4 neighbors of a cell'''
    neighborList = []
    neighborList.append(listA[r-1][c])
    neighborList.append(listA[r][c-1])
    neighborList.append(listA[r][c+1])
    neighborList.append(listA[r+1][c])
    #count the 1's in the list
    return neighborList.count(1)

The function for the level needs to go through every row and column, check the neighbor number of the cell, and create the next level's cell. I called it "wolfram" after the author who inspired the exploration. The "enumerate" function is getting a workout.

def wolfram(level):
    '''returns the array for level'''
    #make the level 0 array
    tableSize = 11
    surface = numpy.zeros((tableSize,tableSize))
    surface[int(tableSize/2)][int(tableSize/2)] = 1
    
    for i in range(level): #repeat for each level
        newSurface = []     #empty list for next level
        for r,u in enumerate(surface): # for each row in array
            newSurface.append([]) #empty list in next level
            for c,v in enumerate(u): # for every item in row
                #if it has exactly 1 or 4 neighbors that are 1's
                if nbs(surface,r,c) in [1,4] or v == 1:
                    newSurface[r].append(1) #put a 1 in next level
                else: #otherwise
                    newSurface[r].append(0) #put a 0 in next level
        #replace the old level with new one for next loop
        surface = newSurface 
    return numpy.array(surface)

Unfortunately, this gives me an error when checking the cells on the sides of the array.

IndexError: index 11 is out of bounds for axis 0 with size 11

The "nbs" function has to be adjusted so it won't panic if there isn't a previous row or next column. The exception keyword in Python comes in handy again. Now I can tell it to "try" searching the neighbors and if it gets a urge to throw an IndexError at me, just move on.

def nbs(listA,r,c):
    '''puts 4 neighbors of a cell into a list'''
    neighborList = []
    try:
        neighborList.append(listA[r-1][c])
        neighborList.append(listA[r][c-1])
        neighborList.append(listA[r][c+1])
        neighborList.append(listA[r+1][c])
        
    except IndexError:
        pass
    return neighborList.count(1)

Now printing "wolfram(3)" gives me an array representing the 4th step in the first picture:

On to Higher Dimensions

So now I can generate an array for any level I want. Time to make a huge array and place blocks for each 1 in that level. I take my file and add the code necessary to make and display a 3D model in Pi3D, my favorite 3D graphics package. Here's the relevant code for taking the levels and making them layers of the pyramid:

cubeList = [] #list for storing cubes

#loop over each level:

for i in range(50,-1,-1): #50 levels, going down to 0

    rowList = wolfram(i)  #get the array for that "layer"

    for row,value in enumerate(rowList): #loop through the rows

for j,v in enumerate(value): #loop through the columns

   if v == 1: #if the value of that cell is 1

                #put a cube there

cubeList.append(pi3d.Cuboid(x=row,

                                            y=-i,

                                            z=j))

Then it's just a matter of iterating over the cube list and drawing all the cubes. Here's the final outcome in pi3d:

Finally I fire up the Raspberry Pi version of Minecraft and use the "setBlock" function to create the pyramid in a 3D world!