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

A285494

The On-Line Encyclopedia of Integer Sequences gets regularly updated with new integer sequences. One of the recent updates was contributed by me, A285494.

A285494 is the list of all numbers k so that its digit sum equals its number of distinct prime factors.
A number’s digit sum is the sum of all of its decimal digits. The number 62831853, for example, has a digit sum of 6+2+8+3+1+8+5+3 = 36.
A number’s number of distinct prime factors is the number of different prime numbers that multiply together to result in the original number. As an example, 62831853 = 3^2 \cdot 7 \cdot 127 \cdot 7853, so it has five prime factors of which four are distinct.
Thereby one can conclude that 62831853 is not an entry in this sequence, as 36 \neq 4.

The sequence is certainly infinite, as the number k = 2 \cdot 10^n with n \in \mathbb{N}^* has a digit sum of 2 + (0 \cdot n) = 2 and — because k = 2^{n+1} \cdot 5^n — exactly two distinct prime factors.

In the encyclopedia entry, I provided a Mathematica one-liner to compute the first few entries of this sequence. Since then, I have also written a Python two-liner to achieve the same goal.


(* Mathematica *)
Select[Range[2,10000],Total[IntegerDigits[#]]==Length[FactorInteger[#]]&]
Out = {20, 30, 102, 120, 200, 300, 1002, 1200, 2000, 2001, 2002, 3000, 3010}
# Python 2.7
>>> def p(n):exec"i,f=2,set()\nwhile n>1:\n\tif n%i==0:f.add(i);n/=i;i=1\n\ti+=1";return len(f)
>>> print filter(lambda n:p(n)==sum(map(int,str(n))),range(2,10001))
[20, 30, 102, 120, 200, 300, 1002, 1200, 2000, 2001, 2002, 3000, 3010]