This is Words and Buttons Online — a collection of interactive #tutorials, #demos, and #quizzes about #mathematics, #algorithms and #programming.
When I started with Python, there were not yet NumPy or SymPy. We used to do research in MatLab, rapid prototyping in Delphi and, believe me or not, symbolic computation in Prolog. Python was something like a nicer Perl, not really popular with the researchers. But it was fun to write things in.
So we did. I choose to implement a couple of experiments for my thesis in Python, and I had to bring a bit of linear algebra with it. I started with writing some Pascal-style loop over loop routines, which was not very exciting, but the more I got into the Python, the more pleasurable it got to become. At some point, I felt that there is something more in the language than loops over loops and started researching the language instead of my original topic. Probably, the best wrong decision I ever took.
Python is famous for its gentle learning curve, but this may very well fire back. You can do Pascal in Python, you can do C++ in Python, you can do Lisp in Python, but you really get the leverage only by doing Python in Python. Yet it is too easy to live with the very basic knowledge. You may be missing all the awesome parts and still feel fine about it.
This article is supposed to show you how powerful Python is at its core. To illustrate this, I chose the domain of linear algebra. Of course, Python has very decent packages for it so you probably wouldn't have to reimplement anything described here. However, I only chose linear algebra due to my lack of imagination. The syntax that allows us to put so much sense in so little code is universal.
Disclaimer: for the sake of readability not all the code samples are actually written in a single line each. It is, of course, possible since neither of them contains Python control structures that rely on indentation. All the indentation in the examples is entirely voluntaristic.
List comprehensions are the staple of Python one-liners. It is a special syntax that describes list transformations. Let's say we want to multiply a vector by a scalar. In Pascal-ish Python this would be something like this.
def scaled(A, x):
B = list()
for i in range(len(A)):
B.append( A[i] * x )
return B
|
Which is fine. This style has its benefits too, for example, you always have a line to put a breakpoint to. But this is overly verbose. In Python, you can simply do this.
def scaled(A, x): return [ai*x for ai in A] |
And it works like that.
List comprehension [1.0, 4.0, 3.0] * 2.0 [, ... |
List comprehensions are not endemic to Python, Haskell and Clojure have them too. Even C# LINQ provides a special syntax for ranges.
The less well-known Python features are dictionary comprehensions and tuple comprehensions.
>>> [2*v for v in [1.0, 2.0, 3.0]]
[2.0, 4.0, 6.0]
>>> {k:2*v for k, v in {0:1.0, 1:4.0, 2:3.0}.items()}
{0: 2.0, 1: 8.0, 2: 6.0}
>>> tuple(2*v for v in (1.0, 4.0, 3.0))
(2.0, 8.0, 6.0)
|
This feature lets you iterate several iterables as one making everything into a list of tuples. Well, it says “list” in the title, but it also works with tuples, dictionaries, generators, everything you can iterate through.
>>> A = [1, 2, 3]
>>> B = ('a', 'b', 'c')
>>> C = {1:'a', 2:'b', 3:'c'}
>>> D = xrange(123)
>>> zip(A, B, C, D)
[(1, 'a', 1, 0), (2, 'b', 2, 1), (3, 'c', 3, 2)]
|
You can make a vector addition in one line with it.
def sum_of(A, B): return [ai+bi for (ai, bi) in zip(A, B)] |
On a pair of lists, it figuratively zips data just like a real zipper.
List zipping A = [1.0, 4.0, 3.0] B = [7.0, 3.0, 1.0] zip(A, B) = [... |
This function simply sums everything for you. You can do a vector dot product in a single line by accumulating multiplications of vectors' elements.
def dot_of(A, B): return sum([ai*bi for (ai, bi) in zip(A, B)]) |
Sometimes, however, a simple summation is just too simple. Like with floating-point numbers, when you face that annoying small error now and then. To make thing more comfortable for you, Python has another sum function that operates on partial sums resulting in more precise output.
>>> [0.1]*10 [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] >>> sum([0.1]*10) 0.9999999999999999 >>> import math >>> math.fsum([0.1]*10) 1.0 |
Conditional expressions are the Python's take on ternary operator. They are expressions that depend on conditions. Expressions in a sense that they expected to be simple and pure. We can use one to make an identity matrix.
def identity(n):
return [[1.0 if i==j else 0.0
for j in range(n)]
for i in range(n)]
|
It's just a nested list comprehension with a conditional expression that works like this.
Conditional expression: 0 == 0 j = 0 [[1.0, ... i = 0 |
Conditional expressions are great for brevity and expressiveness, but since they are meant to be expressions they may not be easy to add side effects to. Not that it isn't possible, it just requires a little trick. In Python, a tuple computation implies computation of every tuple element from left to right. So whenever you need to add more things to your expression, just use a tuple.
You also can't do print from expressions, since they should not have side effects, but sys.stdout.write is somehow fine.
def identity(n):
return [[1.0 if i==j else (0.0,
sys.stdout.write("i != j\n")
)[0]
for j in range(n)]
for i in range(n)]
|
The same trick works for lambdas. You can make your lambdas as large as you wish with tuples. Although, please don't do this unless you have a very good reason to.
Ok, I almost hear some guy from Reddit yelling: “your example doesn't need a conditional expression, to begin with!” Well, it's true. In Python, you can explicitly convert Boolean type variable to a float.
def identity(n):
return [[float(i==j) for j in range(n)] for i in range(n)]
|
Let's just call it a matter of style. I personally like to draw a line between facts and numbers, but it is completely fine not to do that.
Let's say we have a matrix in Python as a list of lists where every nested list is a row. Now if we feed these rows to zip from before, it will form a list of tuples where every tuple is a matrix's column. That's a ready-made matrix transposition.
But we still have to settle two things. One is trivial: we want our matrix as a list of lists, not tuples, so we would have to fix that with a simple list comprehension. The other requires special Python syntax that will allow us to make a list of lists into a tuple of function's arguments. And the syntax for that is simply an asterisk before the matrix.
def transposed(A): return [list(aj) for aj in zip(*A)] |
The difference between passing a list of lists and lists as arguments looks like this.
Passing a container >>> A = [[1, 2, 3], ... [4, 5, 6], ... [7, 8, 9]] >>> def print_A(A): ... print A >>> print_A(A) [[1, 2, 3], [4, 5, 6], [7, 8, 9]] |
The same syntax works the other way around. If you want to receive all the arguments as a single tuple, you can.
>>> def print_arguments_as_tuple(*A): ... print A >>> print_A_as_tuple(1, 2, 3) (1, 2, 3) |
When you keep your matrices as lists of lists, you can easily multiply a matrix by a vector using a single list comprehension with a dot product inside.
def matrix_vector_of(A, X): return [dot_of(ai, X) for ai in A] |
This operation is particularly important in projective geometry. With a single matrix multiplication you can do rotation, scale, translation, affine and projective transformations, — all at the same time. If you want to know more, please take a look at the Interactive guide to homogeneous coordinates.
In computational geometry, a plane is often set by a point in space and a plane's normal vector at this point. With this notation, you can project an arbitrary point on a plane as simple as this.
def projected_on(A, Pn, Pd):
return sum_of(A,
scaled(Pn,
(Pd — dot_of(Pn, A)) / dot_of(Pn, Pn)))
|
We don't need any special syntax for that but hold on; we'll get back to it.
Let's find a Euclidean distance between two points. All we have to do is to perform a vector subtraction, take its self-dot product, and then take a square root of what we got.
def distance_between(A, B):
return pow(dot_of( *(sum_of(A, scaled(B, -1.0)), )*2 ), 0.5)
|
The trick here is to make a self-dot product, that implies passing the same argument to a function without repeating itself. Python has a nice syntax for multiplication in the sense of making something replicate several times. And it is... multiplication. You can use integer multiplication to replicate lists, tuples and so on.
>>> (1, 2, 3) * 2 (1, 2, 3, 1, 2, 3) >>> [1, 2, 3] * 2 [1, 2, 3, 1, 2, 3] >>> ([1, 2, 3],) * 2 ([1, 2, 3], [1, 2, 3]) |
In our case, we need to duplicate a list and feed it to the function.
List concatenation >>> A = [1, 2, 3] >>> def print_A_B(A, B): ... print A, B >>> print_A_B(A*2 ,[]) [1, 2, 3, 1, 2, 3] [] |
Solving linear equations is a respected field of engineering on its own. It is conceptually simple, but you face all kind of problems that nullify the simplicity. Computational error, unreasonable execution time, memory exhaustion, etc.
What I want to show here is that you can, using the power of Python, make a linear solver in one single line of code. Please remember though, that you shouldn't. There are libraries that are much better in this.
The solver itself is explained in another article: Interactive introduction to iterative algorithms. It is an iterative solver that jumps from one hyperplane to another slowly but steadily getting to the point of the solution.
Or it hangs. Or exhaust your stack. It all depends on a system.
def solve(A, B, Xi):
return
Xi if distance_between(matrix_vector_of(A, Xi), B) < 1e-5
else
solve(A[1:] + [A[0]], B[1:] + [B[0]],
projected_on(Xi, A[0], B[0]))
|
The only interesting thing here Python-wise is the list syntax we use to rotate our hyperplane in a projection queue. It is very simple once you get it. You just chop the head off a list, make it into a list itself, and then concatenate the tail with a former head back together.
Rotation: step 1 >>> B = [1, 2, 3] >>> print B[1:] + [B[0]] [2, 3, 1] |
Once again, this is not the best way to invert a matrix. But it works, and it's technically one line.
The idea here is that you divide the identity matrix by the input matrix solving linear equations one unit vector at a time. Linear equation system is basically a matrix by vector division since its left side is multiplication and we get the second operand of it by solving the system.
def inverted(A):
return transpose([solve(A, ort, [0.0]*len(A))
for ort in identity(len(A))])
|
So in just 10 lines of Python code, we have vector addition, vector multiplication, dot product, matrix by vector multiplication, identity matrix, matrix transposition, matrix inversion, point projection, Euclidean distance, and a linear system solver.
Python is very expressive and very powerful not only because of its packages. Of course, modern Python has packages for everything, so you don't have to write a lot of basic code yourself. But, I believe, it only has packages for everything because it is so much fun to write everything in Python.
| Index #programming #mathematics | ← there's more. |
+ Github & RSS |