Re-learning something about scientific computing

For maybe a decade I’ve been wondering why Mathematica could calculate the determinant of a 1,000 by 1,000 matrix with random elements between 0 and 1 in much less than a second when all other systems couldn’t calculate is at all. The value of such a determinant has the order of \(10^{743}\) or more, so a ridiculous large number which, of course, is part of the problem. In Mathematica that calculation is as easy as

In[2]:= A = RandomReal[{0, 1}, {1000, 1000}];
In[3]:= Det@A
Out[3]= 5.413208058576518*10^744

It takes Mathematica about 0.02s to calculate it. Every other system I have tried has failed this task.

The embarrassing thing is that instead of exploring why this is, I’ve used the calculation as a test on how other platforms perform regarding computation on numbers.

That is until now.

So, this is the story of how I went from a naïve Mathematica fanboy believing that Mathematica is just hands down the best tool for scientific computing to re-discovering something I learnt many years ago, shouldn’t have forgotten, and realized something about my own presumptions and assumptions – once again (see Twisting Again).

The story unfolds over some weeks and not just one enjoyable evening of fun hackery. I’m a huge fan of Hammock Driven Development (coined by Rich Hickey in a talk he gave at the first Clojure Conj, in Durham, North Carolina on October 23rd, 2010), and I like to have riddles and challenges such as this rummaging in the back of my mind for weeks.

The story begins

For a while now, I’ve been using LLMs for building small tools and rekindling the fun in programming. I might have forgotten how to program in C, the details of awk scripts, never gotten deep into Python, or macOS Swift programming at all, and so forth. LLMs changed all that, as they let me declare a goal which they then write their code to achieve. Sometimes the first draft is enough, sometimes it takes a few or many iterations – and sometimes I just give up on the task accepting that there’s no substitute for hard work; there’s no free lunch – and that’s okay. Anyway, the point is that I now have my own Little Helper, my own Jarvis. One that has read all the manuals and – or so it thinks – knows it all. So, if I don’t use it for anything important or, if needed, close reads the code, it’s productive – and fun. Suddenly, I can use whatever programming language I want if I only have a bit of experience with it: LISP, Tcl, C, C++, FORTRAN, Visual Basic, AppleScript, JavaScript, Python, Wolfram Language, POSIX tools, etc.

So, the other day, I yet again was thinking about alternatives to Mathematica and once again tried to calculate the determinant of a 1,000 by 1,000 matrix, this time in Julia, which is one of the serious contenders in the scientific computing domain and in some areas a viable alternative to Mathematica.

A new realization

I turned to ChatGPT 5.1 for an implementation of the calculation of the determinant, and for the first time I got much better answers than previously. I don’t know whether that was due to a better LLM, me being better at writing a good prompt, or just spending a bit more time with the technology.

Previously, computational experiments had shown me that naïve determinant calculations break when the dimension of the random matrix with elements \(e \in \lbrack 0;1\rbrack\) exceeds 500, as the absolute value of the determinant becomes bigger than the maximum double precision number, which is around \(10^{308}\). So, the question was: what did Mathematica do and how could I get Julia to do the same?

Anyway, during a longer back and forth with the AI, I learned several things:

When calculating the determinant, you would use LU (or QR) factorization; that the determinant of such large matrices is nonsensical and has no real purpose; if one anyway needs a value for such a determinant, one will always calculate the logarithm of the determinant instead:

\[\log\left| \det(A) \right| = \sum_{i}^{}{\log\left| U_{ii} \right|},\]

where U is the upper part of the LU factorization and then keep track of the sign of the determinant on the side.

The LU factorization and the calculation of log|det(A)| has many years ago been implemented in the grand old LAPACK library and is made available from Julia via the LinearAlgebra package using one simple function call: logabsdet(A). O.M.G! That’s straightforward

ChatGPT also guessed that that’s what Mathematica does. So, no magic is involved!

ChatGPT also suggested another factor in the Julia slowdown as Julia relies on Just in Time compilation (JIT), which can make code run slower on the first pass before the JIT compiler gets to compile the code into machine or byte code, and as I only perform one run, I don’t get the advantage of the JIT.

So, ChatGPT wrote a Julia solution based on the log-magnitude of the determinant, and now, for the first time, I could calculate a determinant of a 1,000 by 1,000 matrix in something other than Mathematica – and I understood why.

What about Python? Completely same story: just use LAPACK through NumPy and calculate the log|det|: np.linalg.slogdet(A).

An old book

Then, in the following days, the LU factorization kept bugging me. I knew it wasn’t new to me. I’d seen it before. Then one morning, before going to the office, I went into an old part of my library (a.k.a. a box in the garage) and found my old Numerical Recipes (FORTRAN) book from my university days: Yes! it had a chapter on LU factorization, and lo and behold! I had been introduced to this method and even written notes in the margins and commented the source code!

A New Dawn – or TIL

After running the Julia and Python solutions, I also explored JavaScript: Same story – just use npm install nlapack and lapack.dgetrf(n, n, A, lda, ipiv); I then also implemented the solution in FORTRAN and C. To no one’s surprise, C and FORTRAN win the speed test hands down, but it is a bit surprising that C is a bit faster than FORTRAN (See Future Work). The Run time column in the below table shows the accumulated time used for the overall computation, including loading the data, re-shuffling it to match the requirements of the LAPACK algorithm but excluding the start-up of the run time – this would add seconds to both Wolfram Language and Julia and be a bit unfair. Still, Mathematica is so slow to start and unusable on the command line (if you don’t use a permanently running WSTP server – another story waiting to be written)

Brevity is another story. Still, the differences in lines of code mostly account for code needed to read a CSV file into memory – easy in Mathematica, Julia, and Python, while you must write everything yourself in JavaScript, FORTRAN, and C. The actual matrix manipulations and calculations are for all languages performed by functions implemented in LAPACK.

As can be seen in the table, all solutions calculate almost the same value (the “random” matrix has been calculated beforehand as described in the code referenced in the Appendix). In the absolute number, which for all solutions is \(\left( 4.35647369541 \pm 3 \cdot 10^{- 12} \right) \cdot 10^{745}\), the difference is at the pico-scale (\(10^{- 12}\)) (se Future Work).

Now, I’ve known for quite some time that Mathematica is definitely not the best tool for integer factorization, but experience from this story here described has in a way re-opened the whole domain of scientific computing systems to me. Instead of being frustrated by why Mathematica leaves the others in the dust, I am now more aware of the fact that it’s all just a bag of tricks. If it seems too good to be true, it probably is. So, explore the tricks and always uses whatever tool “makes the most sense” for any given task.

Future Work

• I’m pretty sure of my results, as the six methods correspond, but still, I would like to go deep into the code to ensure the correctness.

• Why do Mathematica, FORTRAN, and C give the same result, Python and JavaScript the same result, and Julia lies in between these two results?

• How fast would the solution from Numerical Recipes (FORTRAN) (or C) be? I.e., create a solution in C or FORTRAN without using any third-party libraries.

• Really grok the LU factorization again ☺

• Expand the benchmark to be a true benchmark by e.g., averaging multiple calculations etc. As it stands, it just gives an idea of the performance.

• Explore why C is faster than FORTRAN

• How parallel is the LAPACK implementation?

• 1,000 by 1,000 has been my test for a decade, but what about 10,000 by 10,000? ChatGPT has given some hints about going to the GPU for huge matrices. So, that could be fun to try, e.g., using Apple Metal or a NVIDIA H100 on Interactive HPC.

Appendix

On the companion GitHub project 1000x1000, I share the actual codes, the dependencies, and some scripts for running the comparison of the six solutions.