Calculate inverse of dense matrix with entries of very different magnitude

Question

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.

Answer 1

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.