Categories

## Computing Large Determinants in Python

Story:
For my seminar work, I had to calculate the determinant of a large integer matrix.
In this case, large meant n>=500. You might say that this isn’t very large and I would agree. However, it is large enough to overflow a float and reach +inf. Reaching +inf is bad: once you do, you lose all the accuracy of the computation. Also, since determinant calculations use addition as well as subtraction, if +inf is reached in the middle of the calculation the value might never drop back, even if the determinant is a “small” number.

(Note: Even if the calculation hadn’t reached +inf, but just some very large floating point number, the same phenomenon could have occurred.)

As I said, my matrix was composed of integers. Since a determinant is a sum of multiples of matrix elements, you would expect such a determinant to be an integer as well. It would be nice to take advantage of that fact.

How should one compute such a determinant? Well, luckily Python has support for arbitrarily large integers. Let’s compute the determinant with them!
Unfortunately, Python’s numpy computes determinants using LU decomposition, which uses division – hence, floating point values.
Well, never mind, I know how to calculate a determinant, right?
Ten minutes later the naive determinant calculation by minors was implemented, with one “minor” setback – it takes O(n!) time.
Googling for integer determinants yielded some articles about division free algorithms, which weren’t really easy to implement. There was one suggestion about using Rational numbers and Gauss Elimination, with pivots chosen to tame the growth of the nominator and denominator.

So, I hitched up a rational number class, using Python’s arbitrarily large integers as nominator and denominator. Then, I got some code to do Gauss Elimination, although my pivoting wasn’t the most fancy. This seemed to work, and I got my desired large determinants.

Epilogue:
Not too long ago, I ran across mpmath, a Python library for arbitrarily large floating point numbers. It is possible I could have used it instead of my own rational numbers. Next time I run into a relevant problem I’ll keep this library in mind.
Also, an even shorter while ago, I became aware that since 2.6 Python boasts a fractions module. This will sure come in handy in the future.

Code:
Available here. If you intend to use it in Python 2.6 and above, I suggest replacing my Rational class with Python’s Fraction.
(Note that I adapted this module to my needs. Therefore it has a timer that prints time estimations for computations, and it uses Psyco.)

Categories

## You know you need to install a new Python version when…

```>>> import decimal >>> decimal.Decimal('0.2') <&nbsp; 0.3 False >>>```

This little gem took me two hours to track down.
It turns out that since my code is using sqlobject, it also uses the decimal module. I had some constants set up, and from within some larger algorithm I wanted to compare values extracted from the db to those constants. However, my code didn’t seem to work, and I was sure the problem lay somewhere else in the algorithm.
Two hours later, and I found this comparison. It might be similar to this issue, but I’m not sure.

In any case, I was using Python 2.5.2, and decimal works as expected in Python 2.6, so I guess it is indeed time for an upgrade.
After some more checking with Python 2.6, it seems that:

```>>> decimal.Decimal('0.6') <&nbsp; 0.3 True```

Not good.

Categories

## The Case Against Floating Point ==

While working on the clipping bug, I had to do some floating point comparisons. To handle this virtual mine-field, I used simple precision comparisons.

Not surprisingly, I wanted to know what’s ‘the right way’ to do this in Python. When I discussed this subject on #python, I was told to take a look at the decimal module, of which I wasn’t aware at the time.

Most experienced programmers know that you shouldn’t compare floating point numbers with ==. If you want to check floating point equality, you usually decide on a precision, and check either the absolute error, or the relative error. Hence, floating point == isn’t used, maybe except for rare circumstances. I would even venture to say that when possible, static checkers should emit warnings on floating point ==.

On the other hand there are the beginner programmers. Those usually use == by mistake, and will be surprised by the strange results they sometimes get back.

My suggestion is to do either of the following:
1. Change floating point == to behave like a valid floating point comparison. That means using precision and some error measure.
2. Change floating point == to raise an exception, with an error string suggesting using precision comparison, or the decimal module.

Since this change is not backwards compatible, I suggest it be added only to Python 3.
Personally, I prefer no. 2. It is clearer, and less confusing.

Arguments against this suggestion are:

1. This change breaks existing programs:
I believe it most likely triggers hidden bugs. Since the suggestion is to change it only in Python 3, Those programs will most likely be broken by other changes as well, and will need to be changed anyway.

2. This change breaks compatibility with C-like languages:
I agree. However, the already agreed on change of the / operator is even a stronger break. Most arguments for changing the / operator apply here as well.

3. Programmers will still need the regular ==:
Maybe, and even then, only for very rare cases. For these, a special function\method might be used, which could be named floating_exact_eq.

What are your thoughts on the subject?