0
$\begingroup$

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.
Improve this question
$\endgroup$
2
  • 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
    Commented Apr 8, 2021 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
    Commented Apr 8, 2021 at 15:45

2 Answers 2

Reset to default
4
$\begingroup$

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.

Improve this answer
$\endgroup$
3
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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