I need to invert a matrix. Of course, I'm not the first person in this situation, and I know that there's a wealth of powerful libraries out there, of which I only know a couple.
That being said, there is a twist: I know that matrix inversion is (in general) numerically unstable, and I need a solution that is absolutely failsafe, if need be at the cost of large amounts of computing time.
The matrices in question will be relatively moderate in size (not exceeding 100x100 entries), and it's really important to get it right. If the computation takes a couple of minutes on a regular desktop computer, that is completely acceptable. The reason for this is that (a) the matrix inversion will only happen once in the execution of the program, and (b) the analytic solution will only be used as an optional replacement for the already implemented numeric one (mainly for running cross-checks).
I do not think that it matters for the issue at hand, but the matrices will be relatively sparse (10-50% of entries different from zero), and relatively well-behaved (entries in the range 0.001-1000 expected, but nearly identical rows are highly probable).
The solution will be used in C++ code that will be distributed freely, so it is necessary to use either an open (BSD, LGPL, ...) library that can be distributed alongside what I'm doing to a small audience of scientific users (or is ideally already installed on most linux systems), or to implement something from scratch. In that case, a clean, short and less convoluted implementation is to be given preference over an optimized one.
Im happy about:
- pointers to existing libraries that I could use for this
- names or outlines of algorithms that could be useful
- additional ideas or considerations that will help me pursue this endeavour (i.e. will I need to define a custom arbitrary-precision-float type?)
Edit
For illustration purposes, I have added a 15x15 matrix and it's inverse using unsatisfactory numeric inversion.
Here's the input matrix:
$\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0.01513 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0.003782 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.003782 \\ 0 & 0 & 0 & 0.0625 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.0625 & 0 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.01538 & 0 \\ 0 & 0 & 0 & 0.0625 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.01537 & 0 & 0 \\ 0 & 0 & 0.003782 & 0 & 0 & 0.003782 & 0 & 0 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0.003782 & 0 & 0 & 0.003782 & 0 & 0 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0.003782 & 0.003782 & 0 & 0.003782 & 0.003782 & 0 & 0.003782 & 0.003782 & 0 & 0 & 0.003782 & 0 & 0 & 0.003782 & 0 & 0 & 0.003782 \\ 0.0625 & 0 & 0 & 0.0625 & 0 & 0 & 0.0625 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.0625 & 0 & 0.003782 & 0.0625 & 0 & 0.003782 & 0.0625 & 0 & 0.003782 & 0 & 0.01537 & 0 & 0 & 0.01537 & 0 & 0 & 0.01538 & 0 \\ 0.0625 & 0.003782 & 0 & 0.0625 & 0.003782 & 0 & 0.0625 & 0.003782 & 0 & 0.01537 & 0 & 0 & 0.01537 & 0 & 0 & 0.01537 & 0 & 0 \\ 0 & 0 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0.01513 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0.003782 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.003782 & 0 & 0 & 0 \\ 0.25 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.0625 & 0 & 0.003782 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.01537 & 0 & 0 & 0 & 0 \\ 0.25 & 0.01513 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.0615 & 0 & 0 & 0 & 0 & 0 \\ \end{matrix}$
Here's the numeric inverse:
$\begin{matrix} -0 & -0 & -0 & 0 & -0 & 0 & -0 & -0 & -0 & 0 & 0 & -0 & 0 & 0 & -0 & 4 & 0 & 0 \\ -0 & -0 & -0 & 0 & -0 & 0 & -0 & -0 & -0 & 0 & 0 & -0 & 0 & 66.1 & -0 & 0 & 0 & 0 \\ -0 & -0 & -0 & 0 & -0 & 0 & -0 & -0 & -0 & 0 & 0 & -0 & 264.4 & 0 & 0 & 0 & 0 & 0 \\ -0 & 0 & -0 & 16 & -0 & 0 & -0 & 2.935e^{-17} & -0 & -1.805e^{-15} & 0 & 0 & -2.037e^{-31} & -7.338e^{-18} & -1.907e^{-31} & 4.513e^{-16} & 3.944e^{-31} & 0 \\ -0 & 264.4 & -0 & 0 & -0 & 0 & -0 & -0 & -0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 66.1 & 0 & -0 & 0 & -0 & 0 & -0 & -0 & -0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -0 & -16 & -0 & 0 & -0 & 0 & -0 & 16 & 0 & 0 & 1.835e^{-15} & -4.441e^{-16} & 1.718e^{-15} & -4 & -3.553e^{-15} & 0 \\ 0 & -264.4 & -0 & 0 & -0 & 0 & -0 & 264.4 & -0 & 0 & 0 & 0 & 0 & -66.1 & 0 & 0 & 0 & 0 \\ -66.1 & 0 & -0 & 0 & -0 & 0 & 264.4 & 0 & -0 & 0 & 0 & 0 & -264.4 & 0 & 0 & 0 & 0 & 0 \\ 0 & 65.04 & -0 & 65.04 & -0 & -65.04 & 0 & -65.04 & -0 & -65.04 & 0 & 65.04 & -7.338e^{-15} & 16.26 & -6.872e^{-15} & 16.26 & 1.421e^{-14} & -16.26 \\ 16.26 & 0 & -0 & 65.04 & -65.04 & 0 & -65.04 & 0 & -0 & -65.04 & 65.04 & 0 & 65.04 & 0 & 0 & 16.26 & -65.04 & 0 \\ 66.1 & 264.4 & -264.4 & 0 & -0 & 0 & -264.4 & -264.4 & 264.4 & 0 & 0 & 0 & 264.4 & 66.1 & -264.4 & 0 & 0 & 0 \\ 0 & 0 & -0 & 0 & -0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -16.26 & 0 & -16.26 & 0 & 16.26 \\ 0 & 0 & -0 & 0 & -0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -65.04 & -7.105e^{-15} & 1.421e^{-14} & -16.26 & 65.04 & 3.553e^{-15} \\ 0 & 0 & -0 & 0 & -0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -264.4 & -66.1 & 264.4 & 0 & 0 & 0 \\ 0 & -65.04 & -0 & -65.04 & -0 & 65.04 & 0 & -1.193e^{-16} & 0 & 7.338e^{-15} & 0 & 0 & 8.28e^{-31} & 2.983e^{-17} & 7.498e^{-31} & -1.835e^{-15} & -1.578e^{-30} & 0 \\ -16.26 & 0 & -0 & -65.04 & 65.04 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -66.1 & -264.4 & 264.4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{matrix}$