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
	overloading and basic matrix operations.

	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)

	# add two matrices
	def __add__(A, B):
		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

		# advance row
		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)])
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s