# Fastest way to solve a system of linear equations

I have to solve a system of up to 10000 equations with 10000 unknowns as fast as possible (preferably within a few seconds). I know that Gaussian elimination is too slow for that, so what algorithm is suitable for this task?

All coefficients and constants are non-negative integers modulo p (where p is a prime). There is guaranteed to be only 1 solution. I need the solution modulo p.

A LU decomposition of a $n \times n$ matrix can be computed in $O(M(n))$ time, where $M(n)$ is the time to multiply two $n \times n$ matrices. Therefore, you can find a solution to a system of $n$ linear equations in $n$ unknowns in $O(M(n))$ time. For instance, Strassen's algorithm achieves $M(n) = O(n^{2.8})$, which is faster than Gaussian elimination. See https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion.

Rather than trying to implement this yourself, I would suggest using a library: e.g., a BLAS library.

• Also reduce modulo p at the end of the computation. – fade2black Jun 11 '17 at 5:17
• @fade2black, actually, it's probably going to be far better to use an implementation that is designed for use with mod p arithmetic (i.e., is reducing mod p at each step, not only at the end). – D.W. Jun 11 '17 at 6:39
• In case the Wikipedia link changes, the reference given there for the $O(M(n))$ result can be found e.g. in section 28.2 of the 3rd edition of Cormen et al, Introduction to Algorithms. Specifically they show the "algorithmic equivalence" between matrix multiplication and matrix inversion. But presumably one can then link matrix inversion and LU decomposition. – Chill2Macht Feb 19 at 20:55

There is what you want to achieve, and there is reality, and sometimes they are in conflict. First you check if your problem is a special case that can be solved quicker, for example a sparse matrix. Then you look for faster algorithms; LU decomposition will end up a bit faster. Then you investigate what Strassen can do for you (which is not very much; it can save 1/2 the operations if you multiply the problem size by 32).

And then you use brute force. Use a multi-processor system with multiple threads. Use available vector units. Arrange your data and operations to be cache friendly. Investigate what is the fastest way to do calculations modulo p for some fixed p. And you can often save operations by not doing operations modulo p (result in the range 0 ≤ result < p) but a bit more relaxed (for example result in the range -p < result < p).

The best way to solve big linear equations is to use parallelisation or somehow to distribute computations among CPUs or so.

See CUDA, OpenCL, OpenMP.

A lot of people suggests Strassen's algorithm but it has a very big hidden constant which makes it inefficient.

By the way: your linear equations might be very sparse(a lot of zeros), there are few very neat optimisations to solve them in parallel.

• The matrix size is 10,000 by 10,000 so I assume Strassen would be able to save something. Just not very much. – gnasher729 Jun 14 '17 at 18:16
• @gnasher729 I've some doubts, Alex Stapanov in one of his lectures it is mentioned that Boing was using Strassen's algorithm for really big matrices(1Mx1M afair) and they were unhappy with performance. But I think this info is kinda outdated for modern hardware. – Oleg Kovalov Jun 14 '17 at 18:38