Checking Your Work With NumPy
Let’s start this off with a minor confession: despite spending over half a decade working in machine learning (and spending three years developing actuarial software before that, another type of mathematical computing), I didn’t do so well in mathematics for my first few years of college. In fact, I even failed differential calculus the first time I took it, at the College of Staten Island, an institution not necessarily known for its academic intensity. Obviously, I sorted this out, but despite my career experience, some of my mathematical grounding isn’t as solid as I would like.
Working in machine learning, and with data scientists, I’ve always wanted to sharpen my linear algebra knowledge (I aced differential calculus when I took it again, as well as integral and multivariate :) ). In university, it always felt like a very “mechanical” subject to me: you do a lot of simple algebra and hope you don’t make a mistake that propagates down, whether that’s an actual computational mistake or just writing a 1 in a weird way and mistaking it for a 7 or something later. I did fine, but I never really internalized it and felt mentally comfortable with it in the way I did calculus, logic, set theory, or probability. I really deeply believe that anyone can learn mathematics, and it’s a matter of time, persistence, and finding instructional material that just “clicks” with you. To that end, I got a copy of Dr. Gilbert Strang’s Introduction to Linear Algebra, and proceeded to look at it once every six months as life, jobs, cats, and everything else jumped ahead on my priority list.
But recently, I’ve felt that that’s no excuse — I can put some time to it, can’t I? To that end, I’ve been working through Dr. Strang’s textbook a bit at a time, doing most of the exercises from each section until I have reached a level of understanding. It’s slow going, but progress is progress.
One important part of this is actually checking my work. The aforementioned small errors propagating is certainly an issue when learning linear algebra, and I don’t want to cop out by just checking the back of the textbook or plugging everything in to any sort of solver. The happy middle ground, to me, is to use NumPy to check that I’ve, e.g., found the right solution for a specific Ax=b. Doing so is actually pretty simple, but it took me a day to realize that I could do it, so I decided to write up some words on that.
One problem in Section 1.3 of Strang gives a few matrices where one element c is undefined and the reader needs to pick a c to make the matrix singular. The way I solved that is fairly simple; take the rows/columns where all elements are defined and see how the defined elements of the remaining, “free” rows/columns can be expressed as a linear combination of them. Let’s take this example:
# c1 c2 c3
[
[1, 3, 5],
[1, 2, 4],
[1, 1, x],
]Here, c3 = 2c1 + c2. This is, of course, not considering the undefined x . So, if c3 really does equal 2c1 + c2, then x=3, c3 = [1, 1, 3] will “solve” the matrix by making it singular. And we can double check that; numpy.linalg.inv should throw an error if it does. By example:
import numpy as np
A = np.array([
[1, 3, 5],
[1, 2, 4],
[1, 1, 3],
])
np.linalg.inv(A)Running this, you’ll get:
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
Cell In[2], line 6
1 A = np.array([
2 [1, 3, 5],
3 [1, 2, 4],
4 [1, 1, 3],
5 ])
----> 6 np.linalg.inv(A)
File <__array_function__ internals>:200, in inv(*args, **kwargs)
File ~/.cache/pypoetry/virtualenvs/linear-algebra-pdghZUB5-py3.10/lib/python3.10/site-packages/numpy/linalg/linalg.py:538, in inv(a)
536 signature = 'D->D' if isComplexType(t) else 'd->d'
537 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 538 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
539 return wrap(ainv.astype(result_t, copy=False))
File ~/.cache/pypoetry/virtualenvs/linear-algebra-pdghZUB5-py3.10/lib/python3.10/site-packages/numpy/linalg/linalg.py:89, in _raise_linalgerror_singular(err, flag)
88 def _raise_linalgerror_singular(err, flag):
---> 89 raise LinAlgError("Singular matrix")
LinAlgError: Singular matrixAnd we can check whether a matrix is invertible, too, of course. In that case, it should not throw that LinAlgError, but if we’re sanity checking, it should give the same inverse matrix we worked out by hand (with maybe some floating-point rounding errors). In problem 12 of the same section, Strang shows a “centered difference matrix” in R4:
[
[ 0, 1, 0, 0],
[-1, 0, 1, 0],
[ 0, -1, 0, 1],
[ 0, 0, -1, 0],
]By hand, you can work out that x2=b1, -x3=b4, and then substitute in to solve for x1 and x4 in terms of the b’s. np.linalg.inv will also invert this matrix. However, Strang also includes a problem to show that the centered difference matrix in R5 is noninvertible, by selecting a nonzero vector b that solves C5*b=0 . By hand and by looking at the matrix, I figured out that c1+c3+c5=0, which you can double check:
C5 = np.array([
[ 0, 1, 0, 0, 0],
[-1, 0, 1, 0, 0],
[ 0, -1, 0, 1, 0],
[ 0, 0, -1, 0, 1],
[ 0, 0, 0, -1, 0],
])
# remember, 0-based indexing with most PLs, 1-based indexing in mathematics
C5[0, :] + C5[2, :] + C5[4, :]This’ll give you array([0, 0, 0, 0, 0]) , which is clearly the answer, but, once this extends into even higher dimensions, visually inspecting an ndarray becomes annoying and error-prone. So you could make it a little easier on yourself by running something like ((C5[0, :] + C5[2, :] + C5[4, :]) == np.zeros(5)).all() . This one-liner is sensitive to differences in matrix dimensions, though, so the next step here is to throw together a few modules of helpers for it.
I find this pretty cool, because many programming language tutorials start off by showing you how to use the language as a basic calculator. But generally, you don’t find yourself doing that — so finding an actual application for it is neat.
Quite a lot of words for something that, looking back, is pretty straightforward. I hope you got some value out of this article.