3
$\begingroup$

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}$

Improve this question
$\endgroup$
16
  • 2
    $\begingroup$ Define "absolutely failsafe". Computers represent real numbers with finite precision approximations and therefore cannot give you an exact inverse in almost any case. How much error are you willing to tolerate? $\endgroup$
    – Bill Barth
    Jul 1, 2015 at 23:48
  • 1
    $\begingroup$ It is my understanding that for small matrices, also symbolic calculations are possible (e.g. like WolframAlpha does it). Of course, I don't know their algorithmics, but whatever they do should suffice. If you are looking for a numeric quantification, I would say that the worst tolerable deviation for extremely malformed input should be within a permill of the "true" value. Please note that this relative specification of uncertainty implies intentionally that 0-values should come out as 0. $\endgroup$
    – carsten
    Jul 2, 2015 at 0:01
  • 1
    $\begingroup$ If you are happy with an error of one part in one thousand, you should be fine with a 64-bit representation and almost any LU factorization with pivoting assume your matrix is square and has full rank. There are pathological cases where this could break down, but they aren't that likely. How nearly identical are these potential rows? $\endgroup$
    – Bill Barth
    Jul 2, 2015 at 0:09
  • 2
    $\begingroup$ The condition number of your posted matrix is only about 322 meaning you will only lose about two digits of precision during your matrix inversion. For double precision values this still leaves you with about 14 accurate digits. Once again: are you sure this is really necessary? If you really want to validate your results why not just check $M^{-1}M$? $\endgroup$
    – Doug Lipinski
    Jul 2, 2015 at 13:09
  • 3
    $\begingroup$ I personally agree that standard numeric precision is probably sufficient for almost all cases. However, I have been specifically asked to add this as a cross-check for users who don't trust the numerics... :-) $\endgroup$
    – carsten
    Jul 3, 2015 at 3:52

3 Answers 3

Reset to default
2
$\begingroup$

This may be more of a comment than an answer, but I can't comment.

Yes on arbitrary precision. No on symbolic inverse, because at some point, you have to numerically evaluate it, and it is likely to be extremely numerically unstable. Even for 6 by 6 matrix, the symbolic inverse is starting to get pretty ungainly.

Which leaves me the main item. Why do you "need" the matrix inverse? Quite frequently, there is a numerically stable and accurate (and possibly faster) calculation which avoids the need to invert a matrix when a matrix inverse id "needed".

Improve this answer
$\endgroup$
2
  • $\begingroup$ In the case at hand, the actual inversion of the matrix is unfortunately inevitable. The entries of the inverse matrix play the role of prefactors for a linear combination of algebraic expressions. The purpose of this to allow continuously transforming statistical distributions into one another, roughly speaking. $\endgroup$
    – carsten
    Jul 3, 2015 at 3:48
  • 2
    $\begingroup$ I don't understand your answer well enough to know whether actual inversion is necessary. As I wrote before, quite frequently the inverse is not needed when when someone believes it is needed. If you write out in detail how the inverse is used, we may be able to assess whether the inverse is truly needed, or whether there might be a better alternative to accomplish your end goal. Are you sure you can't accomplish your goal using linear equation solves? $\endgroup$
    – Mark L. Stone
    Jul 3, 2015 at 6:38
3
$\begingroup$

There are a number of symbolic algebra packages that will compute the exact matrix inverse; e.g. Sage, Mathematica, Maple, or Maxima. This becomes computationally expensive for large matrices, but there is a chance that a sparse $100 \times 100$ matrix would be within reach. I wouldn't rule it out without trying, since it's so easy to try.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I'm specifically looking for a piece of c++ code or c++-linkable library to do this. $\endgroup$
    – carsten
    Jul 3, 2015 at 20:24
  • $\begingroup$ So, you don't need to invert one matrix, but to have the capability of inverting matrices. That's not clear in your question. If you can try a CAS, you can later generate C++ code. Otherwise, you can try something like GiNaC or Yacas. $\endgroup$
    – nicoguaro ♦
    Jul 3, 2015 at 21:20
1
$\begingroup$

You may want to check division-free Gaussian elimination as described, e.g., in this document titled "A simplified fraction-free integer Gauss elimination algorithm".

I have been told that these approaches are precise until you eventually evaluate the solution (up to loss of significance, which, however, you can monitor). They rather fail because of overflow, which you can handle using appropriate data types.

Improve this answer
$\endgroup$

Your Answer

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

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