# 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. Jun 11, 2017 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, 2017 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. Feb 19, 2020 at 20:55
• An even faster method was described in a recent paper, using a randomized algorithm to solve sparse linear systems. Mar 10, 2021 at 18:01
• @fade2black A helpful approach is not to reproduce exactly modulo p, but have some leeway. Like reduce to the range +/- p. For x % p you need floor(x / p) correctly in all cases. If you allow some error it can get a lot faster. Apr 3 at 19:50

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. Jun 14, 2017 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. Jun 14, 2017 at 18:38
• With 1 million x 1 million you will be unhappy. Strassen increases the computation time by a factor 7 instead of 8 when you double the problem size. If they start with 16x16 and 4096 multiplications and double the problem size 16 times, then instead of 1.15 billion billion products the number is divided by 8.5. Apr 3 at 19:13

In real problems when the size is large ($$n>100000$$) the best approach is the numerical solution of the system by an iterative method like GMRES, CG, BiCG, etc. All those algorithms depends strongly on matrix-vector products, so the complexity is as much as $$O(n^2)$$ which performs better than a direct solution by a LU decomposition with complexity of $$O(n^3)$$.

Sometimes the scale of the problem is not large enough to the proposed methods to be an improvement compared to LU decomposition, so there is a work you have to do in balance the applicability of the methods.

• I think the best approach is to check your hardware budget :-) Apr 3 at 19:14