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__

Continue reading

Advertisements

Arithmetic Golfing

A recent PCG golfing question When do I get my sandwich? asked to find a mapping between seven input strings (sandwich names) and the seven days of the week (indexed by number).

The first answer was made by a user named i cri everytim and utilized a string of characters which uniquely appear at the same position in all seven input strings, enklact, to perform the mapping in Python 2 requiring 29 bytes. After their answer, a lot of answers appeared using the same magic string in different languages to reduce the number of bytes needed. Yet nobody reduced the byte count in Python.

Trying to solve the problem on my own, my first attempt was using only the input strings’ last decimal digit to perform the mapping, though this approach did not save on bytes (read my PCG answer for more on this 30 byte solution).

After a few more hours of working on this problem, however, I achieved to bring down the byte count by one entire byte.

I did so by using a simple brute-force algorithm to check for Python expressions which can be used to perform the sought after mapping. To do so, I use Python’s apostrophes (`...`) to turn the found expression into a string — str(...) is three whole bytes longer — and index that string with the input strings’ lengths. It sure is not very readable, but only takes 28 bytes — and that is all that matters.

lambda S:`6793**164`[len(S)]

After finding the 28 byte function which uses a 9 byte expression (6793**164), I attempted to find an even shorter expression. And even though I did not yet find one, I did write a more general brute-force Python program (source code shown below; can also be downloaded) than the one I linked to in my PCG answer.

Brute-forcing takes exponentially more time the more digits you have to check, so my brute-forcer still requires the user to decide for themselves which expressions should be tried.
There are three parameters that define the search; a regex pattern that should be contained in the expression’s string, an offset that pattern should ideally have and a target length. If an expression is found that takes as many bytes as or less bytes than the target length, an exclamation point is printed.
Though this program did not prove useful in this case, there may come another challenge where an arithmetic expression golfer could come in handy.

My program may not have found shorter expressions, but definitely some impressive ones (the +... at the end refers to an additional offset from the string index which — unsurprisingly — take additional bytes):

  • 2**2**24+800415
  • 2**2**27+5226528
  • 2**7**9+11719750
  • 7954<<850

I also considered using division to generate long strings of digits which may match; the only problem is that Python floating-point numbers only have a certain precision which does not produce long enough strings. Again, using exponentiation (**) and bitshifting (<<) I could not come up with a working expression that takes less bytes.


# Python 2.7 code; 7th, 8th of September 2017

Continue reading

Mandelbrot Set III

I wrote my first ever Mandelbrot Set renderer back in 2015 and used Python to slowly create fractal images. Over a year later, I revisited the project with a Java version which — due to its code being actually compiled — ran much faster, yet had the same clunky interface; a rectangle the user had to draw and a key they had to press to the view change to the selected region.
In this post, over half a year later, I present my newest Mandelbrot Set fractal renderer (download the .jar), written in Java, which both runs fast and allows a much more intuitive and immersive walk through the complex plane by utilizing mouse dragging and scrolling.
The still time demanding task of rendering fractals — even in compiled languages — is split up into a low quality preview rendering, a normal quality display rendering and a high quality 4K (UHD-1 at 3840×2160 pixels to keep a 16:9 image ratio) rendering, all running in seperate threads.

Rainbow spiral
Rainbow spiral

The color schemes where also updated, apart from the usual black-and-white look there are multiple rainbow color schemes which rely on the HSB color space, zebra color schemes which use the iterations taken modulo some constant to define the color and a prime color scheme which tests if the number of iterations taken is prime.

Zebra spiral
Zebra spiral

Apart from the mouse and keyboard control, there is also a menu bar (implemented using Java’s JMenuBar) which allows for more conventional user input through a proper GUI.

Controls

  • Left mouse dragging: pan view
  • Left mouse double click: set cursor’s complex number to image center
  • Mouse scrolling: zoom view
  • Mouse scrolling +CTLR: pan view
  • ‘p’: render high definition fractal
  • ‘r’: reset view to default
  • ‘w’, ‘s’: zoom frame
  • Arrow keys: pan view
  • Arrow keys +CTRL: zoom view
  • Menu bar
    • Fractal: extra info about current fractal rendering
    • Color Scheme: change color scheme and maximum iteration depth
    • HD: controls for high definition rendering
    • Extra: help and about
Blue spiral
Blue spiral

A bit more on how the three threads are implemented.
Whenever the user changes the current view, the main program thread renders a low quality preview and immediately draws it to the screen. In the background, the normal quality thread (its pixel dimensions match the frame’s pixel dimensions) is told to start working. Once this medium quality rendering is finished, it is preferred to the low quality rendering and gets drawn on the screen.
If the user likes a particular frame, they can initiate a high quality rendering (4K UHD-1, 3840×2160 pixels) either by pressing ‘q’ or selecting HD->Render current frame. This high quality rendering obviously takes some time and a lot of processing power, so this thread is throttled by default to allow the user to further explore the fractal. Throttling can be disabled through the menu option HD->Fast rendering. There is also the option to tell the program to exit upon having finished the last queued high definition rendering (HD->Quit when done).
The high definition renderings are saved as .png files and named with their four defining constants. Zim and Zre define the image’s complex center, Zom defines the complex length above the image’s center. Clr defines the number of maximum iterations.

Another blue spiral
Another blue spiral

Just to illustrate how resource intensive fractal rendering really is.
A 4K fractal at 3840×2160 pixels with a iteration depth of 256 would in the worst case scenario (no complex numbers actually escape) require 3840 \cdot 2160 \cdot 256 \cdot 4 = 8493465600 double multiplications. If you had a super-optimized CPU which could do one double multiplication a clock tick (which current CPUs definitely cannot) and ran at 4.00 GHz, it would still take that massively overpowered machine \frac{8493465600}{4 \cdot 10^9} = 2.123 seconds. Larger images and higher maximum iterations would only increase the generated overhead.
The program’s source code is listed below and can also be downloaded (.java), though the compiled .jar can also be downloaded.

Green self-similarity
Green self-similarity

Unrelated to algorithmically generating fractal renderings, I recently found a weed which seemed to be related to the Mandelbrot Set and makes nature’s intertwined relationship with fractals blatently obvious. I call it the Mandel Weed.

Mandel Weed
Mandel Weed
// Java Code; Jonathan Frech; 22nd, 23rd, 24th, 25th, 26th, 27th of July 2017

Continue reading