Python matrix class

Question

To prepare for my studies in linear algebra, I've decided to write a simple matrix class. This is the initial program, with functions to find the determinant and decomposition(s) to follow. I wanted to get feedback first on my approach.

from __future__ import annotations
import operator
import copy
# Don't use [[0] * width] * height, TRAP!
# https://stackoverflow.com/a/44382900/8968906
def fromdim(width: int, height: int):
 return Matrix([[0] * width for _ in range(height)])
class Matrix:
 
 def __init__(self, arrays: list[list]) -> None:
 self.arrays = copy.deepcopy(arrays)
 self.width = len(arrays[0])
 self.height = len(arrays)
 # Public member functions
 def cols(self) -> list[list[int | float]]:
 """
 Returns a list of columns of this matrix.
 """
 return [list(x) for x in zip(*self.arrays)]
 
 def tpose(self) -> None:
 """
 Transposes the current Matrix, modifying the underlying arrays object. This
 converts a MxN array into an NxM array.
 >>> Matrix([[1, 2, 3], [4, 5, 6]])
 [
 [1, 4]
 [2, 5]
 [3, 6]
 ]
 """
 res = fromdim(self.height, self.width)
 for idx, c in enumerate(self.cols()):
 res[idx] = c
 self.arrays = res.arrays
 
 # Dunder methods
 def __add__(self, other: Matrix | float | int) -> Matrix:
 if isinstance(other, (int, float)):
 return Matrix([
 list(map(operator.add, self.arrays[rdx], [other] * self.width))
 for rdx in range(self.height)
 ])
 if isinstance(other, Matrix) and ((self.width, self.height) == (other.width, other.height)):
 return Matrix([
 list(map(operator.add, self.arrays[rdx], other.arrays[rdx]))
 for rdx in range(self.height)
 ])
 def __sub__(self, other: Matrix | float | int) -> Matrix:
 if isinstance(other, (int, float)):
 return self.__add__(-other)
 if isinstance(other, Matrix) and ((self.width, self.height) == (other.width, other.height)):
 return Matrix([
 list(map(operator.sub, self.arrays[rdx], other.arrays[rdx]))
 for rdx in range(self.height)
 ])
 
 def __mul__(self, other: Matrix | float | int) -> Matrix:
 if isinstance(other, (int, float)):
 return Matrix([
 list(map(operator.mul, self.arrays[rdx], [other] * self.width))
 for rdx in range(self.height)
 ])
 if isinstance(other, Matrix) and ((self.width, self.height) == (other.width, other.height)):
 result = Matrix(self.arrays)
 for rdx, row in enumerate(self.arrays):
 for cdx, col in enumerate(other.cols()):
 result[rdx][cdx] = sum(r * c for r, c in zip(row, col))
 return result
 
 def __repr__(self) -> str:
 return '[\n' + '\n'.join([f" {arr}" for arr in self.arrays]) + '\n]'
 def __neg__(self) -> Matrix:
 result = fromdim(self.width, self.height)
 for x in range(self.width):
 for y in range(self.height):
 result.arrays[x][y] = -self.arrays[x][y]
 return result
 
 def __getitem__(self, idx: int) -> list:
 return self.arrays[idx]
 
 def __setitem__(self, idx: int, item: int | float) -> None:
 self.arrays[idx] = item

Testing code:

from PyMatrix import Matrix
m = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
n = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
o = Matrix([[1, 2, 3], [4, 5, 6]])
print("M:", m)
print("N:", n)
print("M + N:", m + n)
print("M - N:", m - n)
print("M * N:", m * n)
o.tpose()
print(o)
print(o.cols())

• Stef: Hi. I'm confused by result.arrays[x][y] in __neg__. Why does x come before y? The array is a list of rows, right, not a list of columns?

Answer 1

Transposes the current Matrix, modifying the underlying arrays object.

You can do it that way. However, it commonly happens in linear algebra code (especially when transcribing formulas from math into code in a fairly literal manner) that you need a matrix and its transpose as well, so I don't think having it as only an in-place modification is convenient. And this in-place modification allocates a whole new matrix anyway, so the main advantage of an in-place transpose is already lost.

Another way to arrange this is returning a "transposed view" of the matrix, then it doesn't have to cost a large matrix copy and / or memory shuffle (of course, that's not a significant consideration if you're only planning to work with small matrices) and there's also still easy access to both the original matrix and its transpose.

• Toby Speight: If the purpose of transposing a matrix is to get better memory locality for an upcoming algorithm, then a view is not what you want. But it's useful in other contexts, so perhaps a good idea to provide both?

Answer 2

I like your implementation, most of it looks good to me. The issues I see are only minor:

1. blank lines

from PEP 8:

Surround top-level function and class definitions with two blank lines.

Here, I would count the two comment lines above the fromdim definition as part of the definition and insert a second blank line between imports and comments, as well as one more blank line between fromdim and Matrix

2. Matrix.tpose()

I think this method should be named transpose.

3. dunder methods

Your dunder methods do multiple things at once. For the __mul__ that makes sense to me, but I would have expected Matrix(...) + 1 to raise a TypeError. However, it looks like you wanted that feature, so I don't mind too much. You could move the code for the matrix multiplication into another dunder method __matmul__, which would allow you to write matrix1 @ matrix2. In that case, I would keep matrix1 * matrix2 as an option for people who don't know about @ or don't like it.

4. module name

From your testing code I can see you called your module PyMatrix. This follows the old naming convention from Python 2. The current guideline prefers a name like pymatrix or py_matrix.

• Stef: I agree for the + 1 feature being weird. The Matrix(...) + 1 feature seems like a weird partial carry-over of numpy's broadcasting features. Outside of such a broadcasting context, I think anyone doing matrix maths would expect Matrix(...) + 1 to either raise a TypeError, or add the identity matrix; since the OP didn't implement any other kind of numpy-like broadcasting, it seems illogical and arbitrary that it behaves the way it does.

Answer 3

pre-condition / class invariant

class Matrix:
 
 def __init__ ...
 self.arrays = ...
 self.width = len(arrays[0])
 self.height = len(arrays)

This is nice, it's very sensible.

Consider iterating over the rows, to verify that each has the identical self.width length. Think of it as "keeping the caller honest". If we were called in error, at least it will be a "shallow" bug, quickly diagnosed.

comments are for the "why", not the "how"

Please elide these comments.

 # Public member functions
 ...
 # Dunder methods

Names such as cols and transpose (you regrettably called it tpose) explain "hello! I am a public name." A leading _ underscore introduces the name of a _private() helper.

When the code has already spoken eloquently, avoid restating the obvious in a # comment.

an int is a float

Well, not exactly. But kind of.

 def cols(self) -> list[list[int | float]]:

Thank you for the very nice annotation, and the docstring. We can simplify it slightly.

 def cols(self) -> list[list[float]]:

The documentation explains:

when an argument is annotated as having type float, an argument of type int is acceptable

This, despite the fact that e.g. isinstance(3, float) is False. It's a minor exception in how we annotate functions, which proves quite convenient.

Similarly for the various other method signatures.

doctest

 def tpose(self) -> None:
 """
 Transposes the current Matrix, modifying the underlying arrays object. This
 converts a MxN array into an NxM array.
 >>> Matrix([[1, 2, 3], [4, 5, 6]])
 [
 [1, 4]
 [2, 5]
 [3, 6]
 ]
 """

Thank you for supplying a doctest docstring. Usually I find them much more believable than ordinary free-form docstrings, because I am confident that the author automatically validated them. But this one? I don't understand its meaning. It seems to contradict the code.

$ python -m doctest matrix.py
**********************************************************************
File "matrix.py", line 33, in matrix.Matrix.tpose
Failed example:
 Matrix([[1, 2, 3], [4, 5, 6]])
Expected:
 [
 [1, 4]
 [2, 5]
 [3, 6]
 ]
Got:
 [
 [1, 2, 3]
 [4, 5, 6]
 ]
**********************************************************************
1 items had failures:
 1 of 1 in matrix.Matrix.tpose
***Test Failed*** 1 failures.

And please call it transpose().

Also, mutating the underlying representation is kind of bizarre. Yes, the docstring accurately tells us what's going on. But for folks accustomed to \$A'\$ notation, or a.T, this is still going to be somewhat counter-intuitive behavior. The one saving grace is we have a -> None: result annotation which helpfully explains that this can only be evaluated for side effects. So a caller who mistakenly attempts a functional approach will soon be told by mypy that they're doing the wrong thing. Good.

implicit None returned

In __add__ there are three cases: {scalar, conforming matrix, other}. In the "other" case, we implicitly return None. Please don't do that.

• The signature's annotation promised we would return a Matrix rather than None.
• Explicitly return None if that's what you intend.

It would be much better to raise if caller offered a bad addend.

Also, how did you expect me to pronounce rdx? After casting about, maybe it is a "right index"? Or a "reflected index"? Please don't use such an abbreviation without a # comment to explain it. Prefer boring conventional index names like i or j.

[Later code suggests that maybe we should pronounce it "row index"?]

commutativity

 def __add__(self, other ...

This is lovely, and m + 0 is the identity function, great.

But what of 0 + m? Please def __radd__() so we don't get "TypeError: unsupported operand type(s) for +: 'int' and 'Matrix'".

matrix multiply

 if isinstance(other, Matrix) and ((self.width, self.height) == (other.width, other.height)):

Please have the """docstring""" clarify whether you were shooting for the dot product operation, a @ matrix multiply operation, or something non-standard. The test suite should include non-square inputs to exercise it. OTOH if this is only defined for the square case, drop one argument or enforce equality with raise or assert.

Also, it seems odd that no @ matrix multiply operation is provided.

zeros

def fromdim(width: int, height: int):

The name is odd. Please call it the familiar matlab / numpy zeros(), or make it a private _fromdim().

dimensions

Clearly this works:

 def __getitem__(self, idx: int) -> list:
 return self.arrays[idx]

But it seems like a missed opportunity. In addition to supporting e.g. m[1][2], consider supporting m[(1, 2)].

(Also, please have the signature explain that we return list[float].)

test suite

This is nice, as far as it goes.

print("M + N:", m + n)

It is much better than nothing. But it does not "know the right answer", so it cannot display a Green bar in automated fashion. Rather than require a human to eyeball the result and judge whether it is correct, please encode that in the test code, perhaps with a call to self.assertEqual().

• Stef: "Please have the """docstring""" clarify that this is the dot product operation, and not a @ matrix multiply operation." <<< This does not look like a dot product. Given the content of the for-loop below, my guess would be that the OP meant to implement matrix multiplication (the one implemented as @ in numpy), but mistakenly wrote condition (self.width, self.height) == (other.width, other.height) instead of self.width == other.height, and wrongly defined the result's dimensions with result = Matrix(self.arrays) instead of result = fromdim(self.height, other.width).
• Stef: As it stands, I think the OP's __mul__ correctly implements matrix multiplication for square matrices, but is just a buggy function for non-square matrices. And indeed the test suite provided at the bottom only tests m * n with m and n being square matrices.
• Stef: Also, I would add that fromdim(height: int, width: int) would seem much much much more natural for me than the OP's fromdim(width: int, height: int)

Answer 4

This signature isn't very precise:

class Matrix:
 def __init__(self, arrays: list[list]) -> None:

list[list] means list[list[Any]], but clearly you want a matrix of floats/ints.

Remember, explicit is better than implicit . Always use parametrized generics, even when all type arguments are Any.

class Matrix:
 def __init__(self, arrays: list[list[float]]) -> None:

Also, be sure to specify Matrix as the return type for fromdim():

def fromdim(width: int, height: int) -> Matrix:
 ...

Without this, Mypy would think fromdim() returns Any. Pyright, on the other hand, can infer the correct type. Rely on neither.

Answer 5

Caching results

A very minor tip: you can cache the result of functions and re-qualify some class methods as properties.

Caching is useful when dealing repeatedly with computationally intensive functions.

All you need is something like:

from functools import cached_property
class Matrix:
 def __init__(self, arrays: list[list[float]]) -> None:
 ...
 @cached_property
 def cols(self) -> list[list[int | float]]:
 ...

And then you use self.cols instead of self.cols(). The relevant doc is here. When you don't want any caching you can just use: @property as a decorator to create a read-only property. Indeed, you don't want the user to mess with this property when playing with your code. So this is desirable behavior.

Since you are studying linear algebra, this is a pattern you surely will want to use later.

NB: there are more caching functions in Python, depending on the intended usage.