I have gotten into writing my own tools, whenever I realize it’s an option.

I will try to document my process of doing this, starting with one of my favourite projects. It’s not pretty, but it did the job for me.

Anyone who has taken a good course on linear algebra knows row operations. They’re fun, they’re easy to program and you can solve systems of equations with them. It’s all great, up until the point where you have to write them.

If you’re an unlucky fellow, you will do this by hand, and no amount of petting the wildebeest will help you. If you’re a, until now, little less unlucky fellow, you will do this with . It’s cumbersome and boring and I grew tired of it within 10 minutes. I guess it’s time for a tool, huh?

Well, first a bit of theory, for the stray reader who has yet to experience the niceties of academia.

### What are row operations?

We have matrices. They are (math people will punch me now) rectangular boxes with numbers in them. You can think of them as proper 2D arrays.

Here’s a matrix for you. It’s a 3×3:

Heres another matrix for you. It’s a 2×4:

If we have a system of linear equations, we can plug them into a matrix.

This will become a 3×4 matrix:

To solve this system, we can use row operations. A row operation is quite simple. You can do a few things:

- Scale a row: Multiply each entry in a row by a common factor
- Add a row to another: For each entry in row , add that entry to the corresponding entry in row
- Subtract: The same as addition, but with a negative attitude
- Swap: Swap two rows.

Why would we do any of this?

Well, if you don’t know linear algebra, you would probably solve that system of equations by substituting terms from one line to the next, mangling stuff around and algebraing all over the place until it hopefully is solved. With row operations, we can do the exact same thing, we just get rid of all the clutter from algebra.

As the purpose of this post is to tell a story of making my computer do my assignments, and not to teach linear algebra, I will get back to the interesting things.

We have our 3×3 matrix from before. This is an identity matrix, but the numbers in the box are completely irrelevant for my example.

Let’s do some row operations!

We will start by scaling. Let us try scaling the first row by a factor of two.

Cool. The entry in is now doubled! Now let’s try adding the new first row to the second row.

This is almost too easy. So what’s the problem?

The problem is, that we should strive to make our look as good as possible, with as little effort as possible. We can write really nice looking row operations with :

What magic did I send to pdflatex to produce this? Cover your face eyes, from hell is coming up.

%Preamble stuff \usepackage[all]{xy} \usepackage{setspace} ... %The code producing nice matrices and arrows \setstretch{1.5} \setlength{\jot}{8pt} \begin{array}{lcl} \left[ \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right] & \xymatrix@C=15ex{ \ar[r]^-{\small \begin{array}{r} \mathbf{r}_1 + \mathbf{r}_2 \to \mathbf{r}_2 \end{array} } & } & \left[ \begin{array}{rrr} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right] \end{array}

Now, to solve a system of equations, you need more than one operation. Probably a lot more. Each of those operations will need a similar 21-line piece of and you need to keep focused and ensure all your matrices have the correct entries. This is way too easy to mess up when you’re copy/pasting.

I can’t explain what every character in those 21 lines do (because I copy/paste!) but I have the grand overview:

We have an array, expecting 3 columns per row. (line 8). Those columns are:

- The original matrix we perform an operation on. (Lines 9-15)
- An xymatrix, a package good for drawing arrows. (lines 17-19)
- The resulting matrix after the row operation. (Lines 20-26)

If we do another row operation, we need to copy the resulting matrix from the previous operation, and make sure that is our starting point.

### Back to writing a tool!

We can now define our problem:

- We want to do row operations. We’d like our computer to do the math, because why not?
- We want to do an arbitrary number of operations starting from one matrix
- We want to generate pretty for all of it

That doesn’t seem too bad. We know what row operations are and we know what the will generally look like.

It’s time to fire up your favourite editor and MacGyver language. Emacs and Python 2.7, I choose you!

First, we need a way to represent a matrix. Python lists are very flexible, and we are not doing anything crazy performancewise. Let’s go with lists.

An empty matrix will be

[[]].

A 1×1 matrix will be [[1]].

The 3×3 identity matrix from earlier will be [[1,0,0],[0,1,0],[0,0,1]].

So far so good, we have a way to make a matrix in Python.

What type should we put in those matrices? I am a fan of fractions, not so much a fan of floating point numbers. It turns out, that python has a built-in fraction type!

We can use it to do things like:

>>> from fractions import Fraction as F >>> F(1,2) Fraction(1, 2) >>> F(1,2) + F(1,2) Fraction(1, 1) >>> F(1,2) + 2 Fraction(5, 2) >>> F(1,2) + 1.2 1.7

This seems reasonable. Python doesn’t care whether we add fractions, floats and ints. If we just remember to convert back to a fraction, we should be good.

*After finishing the project, i realised these conversions are unnecessary. As i try to document my process, and not write a perfect piece of software, i include them anyway. They’re still present in the source as well.
*

How would we go about doing that? Say, we have a matrix full of ints, and we want those ints to be converted to fractions? Easy, list comprehensions.

def intmatrix_to_frac_matrix(intmatrix): """ Converts a matrix of ints to a matrix of fractions.Fraction""" return [[F(ix) for ix in x] for x in intmatrix]

Why does this work? Well, our matrix is not really a matrix. It’s a list of lists.

A list comprehension takes a list as input and applies some function to each element, returning the resulting list. It’s basically a map. Think of it as “For each element ix in the list x apply function f to ix”.

We have a list of lists, so we need a nested list comprehension. This just means that our function to apply is itself a list comprehension. This makes sense since the “first” list only contains lists and those lists contain our fractions.

We now have “For each list x in list intmatrix, do (for each element ix in list x apply f)”.

Our function is here fractions.Fraction, and the result will be a nested list full of fractions.Fraction.

Now we can start defining some row operations. Let’s begin with scale. That seems easy enough.

As input, we need a matrix on which to perform the operation, the row to scale and a scalar to scale with.

def scale_row(matrix, row, multiplier): """Scales a row in supplied matrix""" matrix[row-1] = map(lambda x: x * multiplier, matrix[row-1]) return matrix

I made some low-level decisions here. In math, a matrix is generally 1-indexed. In programming, arrays are generally 0-indexed. Since we are doing this whole project for the sake of math, I will work with the math syntax, designate the first row as row 1, and take care to update my indices whenever necessary.

Since matrix is a list of lists, matrix[0] will give us the first row. We designate the first row with the integer 1, so to manipulate the first row, we subtract 1 from the supplied row.

This gives us a list (of fractions). We can use map on a list. It’s effectively a for loop, but nicer to read. I define a lambda function that takes the element x, which is a fraction, and multiplies it by the supplied scalar. After that, we simply return our matrix.

This is fun! Let’s do another one. Addition!

def add_row(matrix, row1, row2, multiplier): """Add row1 to row2 multiplier times""" for I in range(len(matrix[row1-1])): matrix[row2-1][i] = matrix[row2-1][i] + multiplier * matrix[row1-1][i] return matrix

This is a bit different. We need two lists of equal length, and and for each element in , we need to add the corresponding element from to it. The simplest way to do this I could think of, is simply iterating over the lists.

Subtraction is the same, almost:

def subtract_row(matrix, row1, row2, multiplier): """ Subtract row1 from row2 multiplier times.""" for i in range(len(matrix[row1-1])): matrix[row2-1][i] = matrix[row2-1][i] - multiplier * matrix[row1-1][i] return matrix

We have one operation left, swap rows:

def swap_rows(matrix, row1, row2): """Swap two rows in supplied matrix""" tmp = matrix[row1-1] matrix[row1-1] = matrix[row2-1] matrix[row2-1] = tmp return matrix

It’s about time to test that our functions work. Python 2 has an easy to use testsuite, simply called unittest.

I set up a test skeleton:

import unittest import rowop class add_row_tests(unittest.TestCase): def test_one(self): a = [[1]] expected = [[2]] actual = rowop.add_row(a,1,1,1) self.failUnless(actual == expected) def test_two(self): a = [[1],[1]] expected = [[2],[1]] actual = rowop.add_row(a,2,1,1) self.failUnless(actual == expected) def test_three(self): a = [[rowop.F(-1) for xi in range(3)] for x in range(3)] expected = [[-2,-2,-2],[-1,-1,-1],[-1,-1,-1]] actual = rowop.add_row(a,2,1,1) self.failUnless(actual == expected) def test_four(self): a = [[rowop.F(1,2) for xi in range(10)] for x in range(2)] expected = [[1 for xi in range(10)],[rowop.F(1,2) for x in range(10)]] actual = rowop.add_row(a,2,1,1) self.failUnless(actual == expected) def main(): unittest.main() if __name__ == "__main__": main()

I execute the tests by running

python tests.py and get the output:

.... ---------------------------------------------------------------------- Ran 4 tests in 0.002s OK

It seems the add row operation works as expected. I continue to write similar tests for the other operations, but nobody likes to read 200 lines of test code, so find them on github.

### Making python do the math

from rowop import *

# Matrix from wikipedia

mat = [[ 3, 2, -1, 1],

[ 2, -2, 4, -2],

[-1, F(1,2), -1, 0]]

#Start solving it

scale_row(mat, 3, 2)The scale_row operation updates the values in mat, so after executing these two lines, mat will instead contain:

mat = [[ 3, 2, -1, 1], [ 2, -2, 4, -2], [-2, 1, -2, 0]]

Nice. Let’s add some more operations, and see if we can solve this system!

The way I would go about this, is define a function to print a matrix. I’d type the matrix into a python file, make an operation and print the matrix. Based on the output, i’d add another operation, run the code and print the result, repeating until solved.

Instead of writing up the entire matrix by hand, I now just state the row operation and I get the result. Wonderful!

First, a pretty printer for our matrices:

def print_matrix(matrix): """Prints a matrix""" matrix_string = "" for row in matrix: for col in row: matrix_string += "%5d" % col matrix_string += "\n" print matrix_string

from rowop import * # Matrix from wikipedia mat = [[3,2,-1,1], [2,-2,4,-2], [-1,F(1,2),-1,0]] #Start solving it scale_row(mat, 3, 2) add_row(mat, 3, 2, 1) add_row(mat, 3, 1, 1) add_row(mat, 1, 3, 2) add_row(mat, 2, 3, 4) swap_rows(mat, 3, 2) subtract_row(mat, 2, 1, 1) add_row(mat, 2, 3, F(1,3)) scale_row(mat, 2, F(1,3)) scale_row(mat, 3, F(1,2)) add_row(mat, 3, 1, 3) print_matrix(mat)

This gives the output:

1 0 0 1 0 1 0 -2 0 0 1 -2

Luckily, that system of equations is straight from the wikipedia article on systems of linear equations and they give us the solution: .

### Printing it as LaTeX

Now we need to print it as .

First, since we’re using fractions, we want nice looking fractions, like and . I make a function to convert fractions to .

def frac_to_latex(frac): """Takes a fraction and returns it as a latex-formatted string""" split = str(frac).split("/") if len(split) == 2: return "\\frac{%s}{%s}" % (split[0], split[1]) return split[0]

If the input is a whole number, I want to return that as a string. If it is a fraction, I want to return it as a -formatted fraction. I get the string representation of the fraction and split it. If the fraction is a whole number, the string representation will be the same as the integer representation, and split will only contain one element. If the fraction is indeed a fraction, the string representation is “x/y”. Splitting on “/” gives us each the numerator and demoninator as strings in a two element list, which are then injected into a -formatted string. All backslashes have to be escaped.

Next up, a whole matrix! The first and third part of the -array we want to fill, are matrices, so it’s a good place to start.

def frac_matrix_to_latex(matrix): """Converts a nested list of fractions into a latex-formatted string""" matrix_string = "\\left[\n\\begin{array}{%s}\n" % (len(matrix[0]) * "r") for row in matrix: for col in row: matrix_string += " %s &" % frac_to_latex(col) matrix_string = matrix_string[:-1] matrix_string += " \\\\\n" matrix_string = matrix_string[:-5] matrix_string += "\n\\end{array}\n\\right]\n" return matrix_string

Not much to see here. The number of columns is injected by multiplying the string “r” with the length of the first element of the matrix, i.e. the number of columns.

Then a nested for loop adds each element to the string we’re building.

When the loop has run, we’ve added one too many ampersands to the string. Quick and dirty, it’s removed in line 7.

Then two escaped backslashes, for a newline in , and a newline in the string we’re building are added.

When we’ve added all the entries to the matrix_string, we’ve once again added too much, and it’s removed after the loop. The end of the array is added, and the string is returned.

The next thing that is variable in our , is the operation string above the arrow.

def create_op_string(string_multiplier, row1, row2, operation): """Creates the latex string representing the row operation that will sit atop the arrow between the two matrices""" # fix the multiplier. We don't want' 1r_n or -1r_n if string_multiplier == "1" or string_multiplier == "-1": string_multiplier = "" if operation == "swap": return "\\mathbf{r}_%d \\leftrightarrow \\mathbf{r}_%d" % (row1, row2) if operation == "scale": return "%s\\mathbf{r}_%d \\to \\mathbf{r}_%d" % (string_multiplier, row1, row1) return "%s\\mathbf{r}_%d %s \\mathbf{r}_%d \\to \\mathbf{r}_%d" % ( string_multiplier, row1, operation, row2, row2)

Pretty straight forward. We supply a string representation of the multiplier, i.e. how many times the row operation should be performed, the two rows to operate on, and the operation as a string. If the operation is scale, the second row is just ignored.

Time to use those functions.

def print_row_op(matrix, row1, row2, multiplier, operation): """Perform a row operation on supplied matrix and return it as a latex string""" string_multiplier = frac_to_latex(multiplier) latex_string = ("\\setstretch{1.5}\n" "\\setlength{\\jot}{8pt}\n") latex_string += " \\begin{array}{lcl}\n" # add first matrix latex_string += "%s&\n" % frac_matrix_to_latex(matrix) # Add row operation line latex_string += ("\\xymatrix@C=15ex{\n" " \\ar[r]^-{\\small\n" " \\begin{array}{r}\n") latex_string += create_op_string(string_multiplier, row1, row2, operation) latex_string += "\n \\end{array}\n} & \n}" # add second matrix if operation == "+": add_row(matrix, row1, row2, multiplier) elif operation == "-": subtract_row(matrix, row1, row2, multiplier) elif operation == "swap": swap_rows(matrix, row1, row2) elif operation == "scale": scale_row(matrix, row1, multiplier) latex_string += "\n&\n%s" % frac_matrix_to_latex(matrix) latex_string += "\\end{array}" return (matrix, latex_string)

Here it gets interesting. Since the matrix is updated and not copied on each operation, the function to generate the needs to perform the operation on the matrix. This is not pretty, but neither is duct tape.

First some static is added in lines 6-8.

The -representation of the matrix is generated and added to the string in line 10.

The arrow and row operation string is generated in lines 12-16.

Then the operation is performed on the matrix in lines 18-25.

The for the now updated matrix is added in line 26.

The end of the is added in line 27.

Finally, the matrix and -string is returned in a tuple.

I can now do this:

if __name__ == "__main__": mat = [[3, -1, 5], [-1, 1, 1], [1, -2, -5]] print print_row_op(mat, 1, 3, 1, "swap")[1] print print_row_op(mat, 1, 2, 1, "+")[1] print print_row_op(mat, 1, 3, -3, "+")[1] print print_row_op(mat, 2, 3, 5, "+")[1] print print_row_op(mat, 2, 1, -2, "+")[1]

And have for my assignments printed in my terminal.

I later turned this into a web-app, because I’d never made one. It’s available here.

The full source is on github.