A little disclaimer first. This post is not about linear algebra in Python. If you want exactly that, have a look at python.linalg, it is worth it. This post is about how awesome Python is in its core and how much you can achieve with it writing very little code.
We’ll start slow. Python has list comprehensions. It’s a way you can make one list into another without explicitly specifying loops and indexes. So if you want a function that multiplies a vector by scalar, you can do it as simple as this:
def scale(A, x): return [ai*x for ai in A]
You can do several lists at once with
zip. It takes a couple of lists and makes a single list of tuples out of them. Each tuple then contains exactly one element from every list. So adding vectors is rather plain too.
def add(A, B): return [ai+bi for (ai, bi) in zip(A, B)]
Python also has some list handling functions from the box. Like
sum. Let’s use it to write a function for a
def dot(A, B): return sum([ai*bi for (ai, bi) in zip(A, B)])
List comprehension is great, but you still might want to operate on indexes now and then. Python has a
range function that generates a sequence of indexes for you. There is also conditional expressions in Python, so with these two features you can compose an identity matrix as a nested list.
def ident(n): return [[1.0 if i==j else 0.0 for j in range(n)] for i in range(n)]
zip? We used it to add two vectors, but it can operate on an arbitrary numbers of arrays. Matrix, as a nested list, is an arbitrary list of matrix rows. If we manage to zip it, we will get a list of columns. That’s basically a matrix transposition. Python lets you turn a list into an argument list for a function with a special syntax. You have to write your list with an asterisk in front of it.
def transpose(A): return [list(aj) for aj in zip(*A)]
Of course you can use your own functions to build new ones. So far all of our functions were referential transparent meaning that they always provide the same output for the same input. Like mathematical functions or fine tuned mechanical machines. This is a great trait that enables composition without adding much complexity into the code. We can reuse
dot function from before in a context of matrix by vector multiplication.
def mul(A, X): return [dot(ai, X) for ai in A]
However it can get messy. At some point the code consists more of user defined words than basic syntax. By adding new functions you create your own language, that’s why naming is important.
We already defined 3 words for vector operations, let’s combine them to make a
project function. It projects a point from space onto a plane, or hyperplane, or a line defined by its normal vector and a distance to the center of coordinates.
def project(A, Pn, Pd): return add(A, scale(Pn, (Pd - dot(Pn, A)) / dot(Pn, Pn)))
Python also provides a syntax for multiplying lists by repetitive concatenation. It’s when you write
[1, 2, 3] * 3 and get
[1, 2, 3, 1, 2, 3, 1, 2, 3]. We can use this to define a
distance function from a dot product as a dot product of vector and self is basically a square of its length.
def distance(A, B): return pow(dot( *(add(A, scale(B, -1.0)), )*2 ), 0.5)
Now things can get a little tricky. I want to show you how to solve a linear equation system with a one-liner. Well, it is pretty straightforward with Cramers rule, but you would face serious computational issues using it without Chio condensation, and that might require pivoting, and that is no matter of a one line. I have another idea.
Remember, you can define a plane with a normal vector and a scalar? Well, that pretty much states a linear equation. You can picture a system of linear equations as intersecting planes. When two of them are parallel, there is no intersection and no solution. When two of them are identical, you have a whole line as a solution, so no particular point. But when they intersect, there is the one and only point and you want it.
And this point obviously lies on every single plane at once. Now consider this: a point in space can not be closer to a point on plane, than its own projection on that plane. Because point in space, its projection, and a point on a plane make a right triangle, and hypotenuse can be no shorter than a leg. So if we start from any point in space and will repeatedly project it on every plane from a system, we will eventually come close to the solution. That’s if there is one of course.
That’s how all iterative algorithms are made. You make sure every operation brings you closer to the goal, that’s called convergence, and do these operations until you are done.
But we want this to be one-liner, so no loops are appropriate here. We will use recursion. On every step we will check if we are good with the current solution, and if not — rotate our system by one equation and go on.
As for matrix rotation, Python grants a very convenient syntax for list ranges we will use. You can not only get an element from a list by its index in Python, but also specify a range of elements you want to get with a colon. So all elements of
A apart from the first will be simply
A[1:] and a list of a single first element of an
[A]. Now you just concatenate them like
A[1:] + [A] and you got yourself a rotation.
def solve(A, B, Xi): return Xi if distance(mul(A, Xi), B) < 0.0001 else solve(A[1:] + [A], B[1:] + [B], project(Xi, A, B))
Now let’s finish with computing an inverse matrix. Usually you’d have to find an adjoined matrix of cofactors and divide it to the determinant, but that’s just too much typing. We can find an inverse matrix with one-liner as well.
Consider this. When we multiply a matrix by vector, we have a list of dot products. Now if a vector is an orth, like [1, 0, …], by multiplying matrix by it we will get a list of its first column elements. That the first column. We can get any column by selecting a corresponding orth.
We can get a vector, we would of have by multiplying inverse matrix by something, by solving a linear equation system made of the original matrix, which is itself inverse to the inverse, and that something.
So we will simply solve that system for every orth and then transpose the list of resulting solution to have an inverse matrix.
def invert(A): return transpose([solve(A, ort, [0.0]*len(A)) for ort in ident(len(A))])
That’s it. Only 10 small functions, but so much computation.
Of course Python code is not all that terse. You can pretty much write Fortran-style or Java-style in Python, and by
you I mean everyone else. But as for a mainstream language Python sure has a lot of fun to offer.