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.

  • $\begingroup$ Also reduce modulo p at the end of the computation. $\endgroup$ – fade2black Jun 11 '17 at 5:17
  • 2
    $\begingroup$ @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). $\endgroup$ – D.W. Jun 11 '17 at 6:39

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.

  • $\begingroup$ The matrix size is 10,000 by 10,000 so I assume Strassen would be able to save something. Just not very much. $\endgroup$ – gnasher729 Jun 14 '17 at 18:16
  • 1
    $\begingroup$ @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. $\endgroup$ – Oleg Kovalov Jun 14 '17 at 18:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.