Matrix Operations in Python

Matrices are a major part of math, however they aren't part of regular python. So, I created an easy to use matrix class in python. The code can be found here. It can do a variety of functions, such as addition, subtraction, multiplication, division (multiplying by inverse of another matrix), and solving a system of equations.

I started with two useful functions:

def dot(a,b):
return sum([a[i] * b[i] for i in range(len(b))])
def addPart(array):
return [sum([array[k][i] for k in range(len(array))]) for i in range(len(array[0]))]

The first function returns the dot product of two lists so dot([a,b,c],[d,e,f]) returns [ad, be, cf].
The second function is harder to read, but essentially, given a two dimensional array, it returns an array of the sum of the columns. 

Within the class, I started with the __init__, and __repr__ functions:

def __init__(self, rows):
    self.row = rows
def __repr__(self):
    return '\n'.join([str(i) for i in self.row])

The second function is the result of  printing a matrix, and it returns a row on each line.

I defined the determinant of a matrix as the abs of it, and I wrote it recursively, meaning it could find the determinant of any N x N array.

def __abs__(self):
        if len(self.row) > 2:
            a = 0
            for replace, place in enumerate(self.row[0]):
                sec = Matrix([[i for n,i in enumerate(self.row[j]) if (n != replace)] for j in range(1, len(self.row))])
                a += (-1)**replace * place * abs(sec)
            return a
        elif len(self.row) == 2:
            return (self.row[0][0] * self.row[1][1]) - (self.row[0][1] * self.row[1][0])
        else:
            return self.row[0][0]

This gives three scenarios for determinants: when it's 1 x 1, just return the cell, when it's 2 x 2, it's easy to type out, and anything above that is done recursively. The algorithm for finding a determinant is taking sum of the cofactors of each of the elements in the top row. 

Multiplying, adding, subtracting, negating, and raising to a power are fairly simple, so I'll skip over those, but taking the inverse and solving a system of equations are interesting problems. 

def __invert__(self):
        try:
            determinant = abs(self)
        except:
            return self #Returns, if not invertible
        news = Matrix.transpose(self)
        news = [[(-1) ** (column + row) * abs(Matrix([[news.row[k][i] for i in range(len(self.row[0])) if i != column] for k in range(len(self.row)) if k != row])) for column in range(len(self.row[0]))] for row in range(len(self.row))]
        return Matrix.scalar( Matrix(news),1/determinant)

It is important to realize that not every matrix cannot be inverted, if the determinant of a matrix is 0, it is singular, and it doesn't have an inverse. The way one inverts a matrix is taking the transpose, then taking the matrix of the cofactors. So that long line with " new = [[(-1) ... etc.", is essentially taking the determinant of all of the possible (n-1) x (n-1) matrices (removing one row and one column each time), and multiplying each of them by -1 ** (row + column), in order to negate them when appropriate. Then it multiplies that matrix by 1/determinant. 

def ref(self): #row-eschilon form
        if len(self.row) == 2:
            return Matrix(twobytwo.func(self.row))
        else:
            array = twobytwo.toPos(self.row, 0)
            Barray = Matrix.ref(Matrix([i[1:] for i in array[1:]]))
            array = [array[0]]  + [[0] + Barray[i] for i in range(len(Barray.row))] #For n x m array, it solves the bottom right (n-1) x (m-1)
            for i in range(1,len(array)): #Can't do it inline, 
                array[0] = addPart([array[0], twobytwo.multList(-array[0][i],array[i])]) #Solves for one, but it still won't be simplified
            array[0] = twobytwo.multList(1/array[0][0], array[0]) #simplify
        return Matrix(array)

Many of you may remember I wrote a post about solving systems of equations through row-eschilon form, and in retrospect, I did it very poorly. This way is much better. Once again, it's recursive. When it's a system of two equations, I just used my old algorithm for systems of two equations. For anything else, it takes out the first position of all of the other equations, and it solves the last (n-1) x (m-1) of the array.  Later, it back substitutes, by multiplying the (i+1)'th (i starts w/ 0) array by the value of -array[0][i], which is the the (i+1)th element of the first row. This works because it always eliminates a specific element, because if the matrix is [0,1,0,5], it would multiply it by the negated y value in the first row, and adding it back would remove y.

This is way better than my old way of doing it, and eventually I'll update that post, but for now, this, possibly the biggest computer science innovation of the 21st century, can do all of the Matrix operations very easily.