1
$\begingroup$

I have a very similar problem to the one described in Calculate inverse of dense matrix with entries of very different magnitude. The reason why I am opening a new question is because as far as I understood the answer given to the problem in that question was too specific for the problem considered.

My issue is that I have a matrix C which is fairly large (say typically 500x500 entries) and the entries are spanning some orders of magnitude, typically from $10^3$ to $10^{-16}$ and I need to INVERT this matrix. I really need to do this, as in fact in a following stage I need to use the inverse matrix in some multiplications with other matrices. This is also the reason why I am more interested in the PRECISION of the inversion rather than the performance in terms of speed.

Therefore my obvious question would be: what is the best algorithm to use to get the inverse matrix in this case?

I tried Singular Value Decomposition, but I noticed that it is producing unstable results.

Improve this question
$\endgroup$
4
  • 10
    $\begingroup$ "I really need to do this, as in fact in a following stage I need to use the inverse matrix in some multiplications with other matrices." You should show the details of what you need/intend to do with the result. It sounds like you do not need to invert the matrix, rather, that you can do linear equation solve instead. $\endgroup$ – Mark L. Stone Jun 23 '17 at 16:27
  • $\begingroup$ To echo the previous comment, in the vast majority of cases you'd want to use a solver, not matrix inversion. Do the matrices in question have any particular properties, e.g. diagonally dominant? Have you tried performing the inversion using twice the precision of the data (so if the data is double precision, perform inversion in quadruple precision or double-double, then round final results back to double precision)? I had reasonable success with that approach in one case where the matrix inversion would have been difficult to avoid. $\endgroup$ – njuffa Jun 25 '17 at 1:58
  • $\begingroup$ Mishchenko and Travis. "Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers". J. Quant. Spectrosc. Radiat. Transfer, Vol. 60, No. 3, Sep. 1998, pp 309-324: Since different elements of the matrix Q can differ by many orders of magnitude, the calculation of the inverse matrix Q⁻¹ is an ill-conditioned process [..] We have shown in Ref. 4 that an efficient way of dealing with this numerical instability is to compute and invert the matrix Q using extended precision [..] instead to double-precision $\endgroup$ – njuffa Jun 25 '17 at 6:56
  • $\begingroup$ The stability of SVD may dependent with the threshold setting. Try to add a small regularization parameter. $\endgroup$ – lennon310 Jun 27 '17 at 19:40