# Calculate inverse of dense matrix with entries of very different magnitude

I need to calculate the inverse of a dense matrix, with some elements taking values as high as 1e9 and some around 1e2. What would be the best method to do it?

Note:

• I am more concerned about the precision of the answer rather than minimizing the computational time.

• I explicitly need to find the inverse and not the solution of any linear system Au=b.

Size: The size of the matrix is below $100 \times 100$.

Formula The formula for Pade Approximation of order 3,3 can be given as:

$\left(120I_{n\times n}+60A+12A^2+A^3\right)\times\left(120I_{n\times n}-60A+12A^2-A^3\right)^{-1}$

where,

A is the matrix whose exponential is to be found out.

This is the Pade Approximation for order 3,3. I am trying to implement a Pade approximation of order 13,13 in which the coefficient of $I_{n\times n} \text{is}$ 17643225600

Implementation: I am implementing the Pade Approximation for finding the exponential of a matrix in which I need to find the inverse of the "denominator" matrix.

• @ChristianClason Thanks, I have edited the question. I hope that makes it more relevant and catchty. – user3496912 Mar 4 '16 at 14:05
• @ChristianClason I tried my best, hope that clears it. – user3496912 Mar 4 '16 at 14:27
• @ChristianClason I know, its kind of worthless to find out the inverse of the matrix. But the above problem is my term project and I have to go ahead with it. – user3496912 Mar 4 '16 at 14:42
• No worries, it's just important to be as explicit as possible about any side conditions so you can get the answer need. – Christian Clason Mar 4 '16 at 14:46

Depending on the condition number of the matrix $A$, computing Padé approximations of the matrix exponential $e^A$ is horribly ill-conditioned, especially for such high orders. If a problem is ill-conditioned, no amount of mathematical trickery will allow a stable computation; however, often you can change the problem to one that is better conditioned (and still leads to the desired result). What to do depends on the precise nature of your assignment:

• If the goal is an approximation of the matrix exponential, look at one of the better ways in Clive Moler and Charles van Loan's seminal paper Nineteen dubious ways to compute the exponential of a matrix, SIAM Rev., 45(1), 3–49. (In particular, #2 on page 10ff discusses Padé approximations, and you should consider the scaling and squaring method #3.)

• If it has to be a Padé approximation, there are more robust ways to computing it via SVDs instead of direct inversion; see, e.g. Pedro Gonnet, Stefan Güttel, and Lloyd N Trefethen's paper Robust Padé approximation via SVD, SIAM Rev., 55(1), 101–117. (This paper deals with Padé approximation of real functions, not matrix functions, so some further work is required.)

• Otherwise, for the matrix size and type you are looking at, and since you don't care so much about performance, I don't see a better way than using a direct solver (e.g., LU decomposition) from a high-quality library to solve for each unit vector $e_i$, $i=1,\dots N$ the linear system $B m_i = e_i$ and set $M=[m_1,\dots,m_N]$. Depending on the software framework, you should use the highest precision data type available to minimize round-off errors.

• Regarding the 3rd point is it better to use Gauss-Jordan Method or solve the Gauss Elimination with pivoting n-times. – user3496912 Mar 4 '16 at 15:37
• Mathematically, all methods implemented by direct solvers are equivalent to Gauss elimination with pivoting (possibly exploiting the structure of the matrix to reduce effort). Since you are dealing with such a difficult problem, you really should not implement your own method and rather use a library written by experts (which are tuned down to the specific machine-language instructions to give the best possible stability). – Christian Clason Mar 4 '16 at 15:49