# Python Matrix Module

Matrices are an important part of linear algebra. By arranging scalars in a rectangular manner, one can elegantly encode vector transformations like scaling, rotating, shearing and squashing, solve systems of linear equations and represent general vector space homomorphisms.
However, as powerful as matrices are, when actually applying the theoretical tools one has to calculate specific values. Doing so by hand can be done, yet gets cumbersome quite quickly when dealing with any matrices which contain more than a few rows and columns.

So, even though there are a lot of other implementations already present, I set out to write a Python matrix module containing a matrix class capable of basic matrix arithmetic (matrix multiplication, transposition, …) together with a set of functions which perform some higher-level matrix manipulation like applying Gaussian elimination, calculating the reduced row echelon form, determinant, inversion and rank of a given matrix.

Module source code can be seen below and downloaded. When saved as matrix.py in the current working directory, one can import the module as follows.

>>> import matrix
>>> A = matrix.Matrix([[13,  1, 20, 18],
...                    [ 9, 24,  0,  9],
...                    [14, 22,  5, 18],
...                    [19,  9, 15, 14]])
>>> print A**-1
-149/1268  -67/634   83/1268 171/1268
51/1268 239/1902 -105/1268 -33/1268
Matrix(    73/634 803/4755  -113/634 -87/3170 )
13/1268  -75/634  197/1268 -83/1268


Matrices are defined over a field, typically $\mathbb{F} = \mathbb{R}$ in theoretical use, though for my implementation I chose not to use a double data structure, as it lacked the conceptual precision in numbers like a third. As one cannot truly represent a large portion of the reals anyways, I chose to use $\mathbb{F} = \mathbb{Q}$, which also is a field though can be — to a certain scalar size and precision — accurately represented using fractional data types (Python’s built-in Fraction is used here).

To simplify working with matrices, the implemented matrix class supports operator overloading such that the following expressions — A[i,j], A[i,j]=l, A*B, A*l, A+B, -A, A/l, A+B, A-B, A**-1, A**"t", ~A, A==B, A!=B — all behave in a coherent and intuitive way for matrices A, B, scalars l and indices i, j.

When working with matrices, there are certain rules that must be obeyed, like proper size when adding or multiplying, invertibility when inverting and others. To minimize potential bug track down problems, I tried to include a variety of detailed exceptions (ValueErrors) explaining the program’s failure at that point.

Apart from basic matrix arithmetic, a large part of my module centers around Gaussian elimination and the functions that follow from it. At their heart lies the implementation of GaussianElimination, a function which calculates the reduced row echelon form rref(A) of a matrix together with the transformation matrix T such that T*A = rref(A), a list of all matrix pivot coordinates, the number of row transpositions used to achieve row echelon form and a product of all scalars used to achieve reduced row echelon form.
From this function, rref(A) simply returns the first, rrefT(A) the second parameter. Functions rcef(A) (reduced column echelon form) and rcefS(A) (A*S=rcef(A)) follow from repeated transposition.
Determinant calculation uses characteristic determinant properties (multilinear, alternating and the unit hypercube has hypervolume one).

• $\det \begin{pmatrix} a_{1 1}&a_{1 2}&\dots&a_{1 n}\\ \vdots&\vdots&\ddots&\vdots&\\ \lambda \cdot a_{i 1}&\lambda \cdot a_{i 2}&\dots&\lambda \cdot a_{i n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{n 1}&a_{n 2}&\dots&a_{n n} \end{pmatrix} = \lambda \cdot \det \begin{pmatrix} a_{1 1}&a_{1 2}&\dots&a_{1 n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{i 1}&a_{i 2}&\dots&a_{i n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{n 1}&a_{n 2}&\dots&a_{n n} \end{pmatrix}$
• $\det \begin{pmatrix} a_{1 1}&a_{1 2}&\dots&a_{1 n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{i 1}&a_{i 2}&\dots&a_{i n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{j 1}&a_{j 2}&\dots&a_{j n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{n 1}&a_{n 2}&\dots&a_{n n} \end{pmatrix} = -\det \begin{pmatrix} a_{1 1}&a_{1 2}&\dots&a_{1 n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{j 1}&a_{j 2}&\dots&a_{j n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{i 1}&a_{i 2}&\dots&a_{i n}\\ \vdots&\vdots&\ddots&\vdots&\\ a_{n 1}&a_{n 2}&\dots&a_{n n} \end{pmatrix}$
• $\det \begin{pmatrix}1&a_{1 2}&a_{1 3}&\dots&a_{1 n}\\0&1&a_{2 3}&\dots&a_{2 n}\\0&0&1&\dots&a_{3 n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&\dots&1\end{pmatrix} = \det \bold{1}_n = 1$

Using these properties, the determinant is equal to the product of the total product of all factors used in transforming the matrix into reduced row echelon form and the permutation parity (minus one to the power of the number of transpositions used).
Questions regarding invertibility and rank can also conveniently be implemented using the above described information.

All in all, this Python module implements basic matrix functionality paired with a bit more involved matrix transformations to build a usable environment for manipulating matrices in Python.

# Python 2.7 code; Jonathan Frech; 9th, 10th, 11th, 12th, 14th of December 2017
# 4th of January 2018: added support for matrix inequality, implemented __ne__

"""
Python 2.7 matrix module. Implements a matrix class, complete with operator

Written by Jonathan Frech, December 2017.
"""

# import
import fractions, numbers, copy, random

"""
Matrix class.
"""
# matrix class
class Matrix():
# initiate matrix
def __init__(A, data):
# parse data
if not isinstance(data, list): raise ValueError("Expected list as matrix data.")
if not all(isinstance(j, list) for j in data): raise ValueError("Given matrix data contains non-lists.")
if not len(data) > 0: raise ValueError("Given matrix data contains too few rows.")
if not len(data[0]) > 0: raise ValueError("Given matrix data contains too few columns.")
if not all(len(j) == len(data[0]) for j in data): raise ValueError("Matrix data column size mismatch.")
try: A.data = [[fractions.Fraction(c) for c in r] for r in data]
except ValueError: raise ValueError("Matrix data contains members which cannot be converted to fractions.")

# matrix size (rows, columns)
else: A.r, A.c = A.size()

# return matrix string representation
def __str__(A):
# convert matrix to nested string list (m represents maximum member length, long fractions are represented in decimal form)
m = 10; s = [[str(A[i, j]) if len(str(A[i, j])) < m else str(float(A[i, j]))[:m] for j in A.indexc()] for i in A.indexr()]

# count string lengths, pad strings
l = [max(len(s[i][j]) for i in A.indexr()) for j in A.indexc()]
s = [[s[i][j].rjust(l[j]) for j in A.indexc()] for i in A.indexr()]

# concatenate row strings
s = [" ".join(v) for v in s]

# add 'Matrix( ... )' wrapper
s = [["        %s  ",
"Matrix( %s )"][i==A.r/2] % v for i, v in enumerate(s)]

# return new line concatenated string
return "\n".join(s) + "\n"

# return either matrix row or matrix element
def __getitem__(A, p):
# matrix row
if isinstance(p, int):
if not 0 <= p < A.r: raise ValueError("Member position outside matrix.")
return tuple(A.data[p])

# matrix element
elif isinstance(p, tuple):
if len(p) != 2: raise ValueError("Member position tuple not of length two.")
i, j = p
if not (0 <= i < A.r and 0 <= j < A.c): raise ValueError("Member position outside matrix.")
return A.data[i][j]

# neither integer nor tuple
else: raise ValueError("Member position neither integer nor tuple.")

# return either matrix row or matrix element
def __setitem__(A, p, v):
# matrix row
if isinstance(p, int):
if not 0 <= p < A.r: raise ValueError("Member position outside matrix.")
if not isinstance(v, list): raise ValueError("Expected list.")
if not len(v) == A.c: raise ValueError("Row size mismatch.")
try: A.data[p] = [fractions.Fraction(m) for m in v]
except ValueError: raise ValueError("Expected list of scalars.")

# matrix element
elif isinstance(p, tuple):
if len(p) != 2: raise ValueError("Member position tuple not of length two.")
i, j = p
if not (0 <= i < A.r and 0 <= j < A.c): raise ValueError("Member position outside matrix.")
try: A.data[i][j] = fractions.Fraction(v)
except ValueError: raise ValueError("Expected scalar.")

# neither integer nor tuple
else: raise ValueError("Member position neither integer nor tuple.")

# multiply matrix either by matrix or scalar
def __mul__(A, B):
# matrix multiplication
if isinstance(B, Matrix):
if A.c != B.r: raise ValueError("Matrix size mismatch.")
return Matrix([[sum(A[i, v]*B[v, j] for v in A.indexc()) for j in B.indexc()] for i in A.indexr()])

# scalar multiplication
elif isinstance(B, numbers.Number): return Matrix([[A[i, j]*B for j in A.indexc()] for i in A.indexr()])

# undefined multiplication
else: raise ValueError("Matrix multiplication operand neither matrix nor scalar.")

# negate matrix, multiply by scalar (-1)
def __neg__(A): return A*-1

# divide matrix by scalar
def __div__(A, l): return A*fractions.Fraction(1, l)

if not isinstance(B, Matrix): raise ValueError("Matrix addition operand is not a matrix.")
if A.c != B.c or A.r != B.r: raise ValueError("Matrix size mismatch.")
return Matrix([[A[i, j]+B[i, j] for j in A.indexc()] for i in A.indexr()])

# subtract matrices, add the negated matrix
def __sub__(A, B): return A+-B

# 'matrix exponentiation', transposition or inversion
def __pow__(A, e):
if e == "t": return transpose(A)
elif e == -1: return invert(A)
else: raise ValueError("Expected either transposition or inversion.")

# invert matrix (~A)
def __invert__(A): return invert(A)

# a matrix is considered false if all elements are zero, else true
def __truth__(A): return any(A[i, j] for i, j in A.index())

# two matrices are considered equal iff they are of equal size and all elements match
def __eq__(A, B):
if not isinstance(B, Matrix) or A.size() != B.size(): return False
return all(A[i, j] == B[i, j] for i, j in A.index())

# two matrices are considered unequal iff they are not equal
def __ne__(A, B): return not A == B

# return matrix size (rows x columns, rows, columns)
def size(A): return len(A.data), len(A.data[0])
def sizer(A): return len(A.data)
def sizec(A): return len(A.data[0])

# index rows, colums, rows and columns, colums and rows, default: rows and columns
def indexr(A): return range(A.r)
def indexc(A): return range(A.c)
def indexrc(A): return ((i, j) for i in range(A.r) for j in range(A.c))
def indexcr(A): return ((i, j) for j in range(A.c) for i in range(A.r))
def index(A): return A.indexrc()

# copy matrix
def copy(A): return copy.deepcopy(A)

"""
Gaussian elimination.
"""
# generate identity matrix
def identity(n):
if not isinstance(n, int) or n < 1: raise ValueError("Invalid identity matrix dimensions.")
return Matrix([[r==c for c in range(n)] for r in range(n)])

# generate matrix E such that E*A has the effect of generating a zero-filled matrix except for the i-th row which is the j-th row in A
def GaussianE(n, i, j):
if not (isinstance(n, int) and isinstance(i, int) and isinstance(j, int)) or not (0 <= i < n and 0 <= j < n): raise ValueError("Invalid matrix dimensions.")
return Matrix([[i==r and j==c for c in range(n)] for r in range(n)])

# generate matrix S such that S*A has the effect of multiplying row i with scalar l in A
def GaussianS(n, i, l):
if not (isinstance(n, int) and isinstance(i, int)) or not 0 <= i < n: raise ValueError("Invalid matrix dimensions.")
if l == 0: raise ValueError("Scalar equals zero.")
return identity(n) + GaussianE(n, i, i)*(l-1)

# generate matrix Q such that Q*A has the effect of adding to the i-th row l times the j-th row in A
def GaussianQ(n, i, j, l):
if not (isinstance(n, int) and isinstance(i, int) and isinstance(j, int)) or not (0 <= i < n and 0 <= j < n): raise ValueError("Invalid matrix dimensions.")
if i == j: raise ValueError("Indices are the same.")
if l == 0: raise ValueError("Scalar equals zero.")
return identity(n) + GaussianE(n, i, j)*l

# generate matrix P such that P*A has the effect of swapping rows i and j in A
def GaussianP(n, i, j):
if not (isinstance(n, int) and isinstance(i, int) and isinstance(j, int)) or not (0 <= i < n and 0 <= j < n): raise ValueError("Invalid matrix dimensions.")
return identity(n) - GaussianE(n, i, i) - GaussianE(n, j, j) + GaussianE(n, i, j) + GaussianE(n, j, i)

# apply Gaussian elimination to matrix
def GaussianElimination(A):
# matrix function
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")

# initialize various variables
M = A.copy()    # copy matrix
n = M.r         # size of invertible matrices T', T'', T''', ..., T such that T*A = T'*T''*T'''*...*A = rref(A)
R = 0           # current row that is getting reduced
T = identity(n) # transformation matrix T
P = []          # list of pivot coordinates
t = 0           # number of transpositions used
l = 1           # product of scalars used

# operate on the entire matrix
while True:
# search pivot
for i, j in M.indexcr():
if M[i, j] != 0 and i >= R: break

# matrix end reached
else: break

# pivot not at the top, move it, record that transposition was used
if i > R: T = GaussianP(n, i, R)*T; M = T*A; t += 1

# set pivot to one, record scalar used and pivot position
l *= M[R, j]**-1
T = GaussianS(n, R, M[R, j]**-1)*T; M = T*A
P.append((R, j))

# clear matrix elements beneath pivot
for r in range(R+1, M.r):
if M[r, j] != 0: T = GaussianQ(n, r, R, -M[r, j])*T; M = T*A

R += 1

# set elements above all pivots to zero
for i, j in P:
for r in range(i):
if M[r, j] != 0: T = GaussianQ(n, r, i, -M[r, j])*T; M = T*A

# return matrix in reduced row echelon form, transformation matrix, pivot coordinate list, number of transpositions and product of scalars used
return M, T, P, t, l

"""
Reduced echelon form functions.
"""
# return matrix reduced row echelon form
def rref(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return GaussianElimination(A)[0]

# return matrix reduced row echelon form transformation matrix
def rrefT(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return GaussianElimination(A)[1]

# return matrix reduced column echelon form
def rcef(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return transpose(GaussianElimination(transpose(A))[0])

# return matrix reduced column echelon form transformation matrix
def rcefT(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return transpose(GaussianElimination(transpose(A))[1])

"""
Transposition, determinant, invertibility, inversion, rank.
"""
# matrix transposition
def transpose(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return Matrix([[A[j, i] for j in A.indexr()] for i in A.indexc()])

# calculate matrix determinant
def det(A):
# matrix function
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")

# non-square matrix implies det(A) = 0
if A.r != A.c: return 0

# apply Gaussian elimination to gather number of transpositions and
# product of scalars needed to generate rref(A)
_, _, _, t, l = GaussianElimination(A)

# calculate determinant
return (-1)**t / l

# determine if matrix is invertible
def invertible(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return det(A) != 0

# return inverted matrix
def invert(A):
if not invertible(A): raise ValueError("Matrix not invertible.")
return rrefT(A)

# count matrix zero rows
def countzerorows(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return sum(all(A[i, j] == 0 for j in A.indexc()) for i in A.indexr())

# count matrix zero columns
def countzerocolumns(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return sum(all(A[i, j] == 0 for i in A.indexr()) for j in A.indexc())

# calculate matrix rank
def rank(A):
if not isinstance(A, Matrix): raise ValueError("Expected matrix.")
return A.r-countzerorows(rref(A))

"""
Various functions.
"""
# generate a randomly filled r x c matrix
def randommatrix(r, c):
if not (isinstance(r, int) and isinstance(c, int) and r > 0 and c > 0): raise ValueError("Expected whole number matrix dimensions.")
return Matrix([[random.randint(-10, 10) for j in range(c)] for i in range(r)])

This site uses Akismet to reduce spam. Learn how your comment data is processed.