Normally, this boards sees a lot of traffic about the most efficient and most powerful solvers for huge linear equation systems. But this time, I have the opposite problem:

I need to implement a solver for a linear equation system in a very limited programming system that offers only basic array functionality and no linear algebra whatsoever. Efficiency and scalability is not an issue as the systems have a single-digit number of unknowns. What is way more important is that it's easy to implement and to fully understand with the most basic linear algebra knowledge.

So I'm looking for the simplest algorithm for:

  • Direct solution of the system Ax = b with a full rank A
  • Iterative solution of the system Ax = b with A not necessarily full rank, but with a good initial guess for x.
  • 1
    $\begingroup$ If your system is single-digit small (less than 10 x 10), closed-form expressions (the way you would do it by hand) for the matrix inverse, followed by a matrix-vector product are good enough. This is as simple as you're going to get. In fact, many sophisticated libraries (e.g. Eigen) do this for small fixed sizes such as 4x4, as it is also the fastest method. $\endgroup$ – Charlie S Apr 8 at 14:15
  • 3
    $\begingroup$ There's no point in using an iterative method on a system of equations with fewer than 10 unknowns. $\endgroup$ – Brian Borchers Apr 8 at 15:45

I'd recommend a BLAS1/2-style implementation of LU decomposition with partial pivoting, along with forward/backward triangular solves. If you only have single digit numbers of unknowns, it should be fine to run dense/O(N^3) algorithms even on structured systems (ie banded/sparse/?). All these algorithms can be written in terms of "axpy" operations, which is a particularly easy BLAS1 primitive to understand/implement (axpy literally just means y = alpha* x + (plus) y, a vector-vector addition with scaling)

As a special case, if single digit means just N=2 or N=3 (maybe you just need coordinate transforms or something), those are small/easy enough to just code up the explicit inverse formulae by hand, which would certainly be less work than LU.


This question is as ill-posed as the matrices you may want to invert :-)

Even the most complicated iterative methods, say GMRES to name an example, are not terribly difficult to implement and require only a couple of hundred lines of code. That's probably only a factor of 4 or 5 worse than the arguably simplest method, Richardson defect correction. But it is far more robust.

Where the real effort lies in complicated solver schemes is the preconditioner. But for problems as small as yours, you probably don't need to bother with that to begin with.


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.