(good) numerical inversion of an almost singular matrix: is it possible?

Question

Ok, so I know that if I have a system that looks like Ax=b it is foolish to solve it by solving for the inverse of A and one should instead use something like Gaussian elimination, or something like that. However, I am trying to calculate H(s) = B*(sI-A)^-1*C, where s is some scalar and A,B,C are matrices (of dimensions such that H(s) is a scalar too, if that helps). Now, this is a large system and I can't do this calculation symbolically, but when I try and plug in the s values that I am interested in, the (sI-A) matrix has really bad condition numbers (something like 10^10) and I highly doubt that my numeric results are in any way accurate --- it's not just the condition numbers, the answer comes out completely unlike what I would expect it to be. Do you have any ideas on either some algorithm that could help me out, or a math trick that would scale things in a such a way that inversion would not run into machine precision problems?

Just so everyone knows, this comes up when trying to calculate a transfer function from the state space representation of a linear system.

Answer 1

What I'm doing in the routine case is either to SVD-decompose or to LDU-decompose the parenthese-matrix (sI-A) and consider the inverses of the triangular and diagonal factors. Then looking at that dot-products gives sometimes hints, how to improve the numerical computation too.
Also I recommend to consider to use Pari/GP, which is a free math software and which allows arbitrary precision over complex numbers and also in matrices.

Answer 2

Let $M = (sI-A) \in \mathbb{R}^{n\times n}, b = B \in \mathbb{R}^{1\times n}, c = C \in \mathbb{R}^{n\times 1}, h = H(s) \in \mathbb{R}.$ You are looking for $$ h = bM^{-1} c.$$

In a Krylov subspace fashion, let $$f(x) = f_0 + f_1 x + \ldots + f_n x^n \in \mathbb{R}[x]$$ be the characteristic polynomial of $M$. By Cayley-Hamilton, $f(M) = 0$. From that, you can easily derive: $$ M^{-1} = -\frac{1}{f_0}(f_{n} M^{n-1} + f_{n-1} M^{n-2} + \ldots + f_1 I),$$ and hence $$ h = bM^{-1}c = -\frac{1}{f_0}(f_{n} bM^{n-1}c + f_{n-1} bM^{n-2}c + \ldots + f_1 bIc).$$

Computing the sequence $$\{ b M^i c \}_{i=0}^{i=n-1}$$ takes a $O(n)$ matrix vector products (iteratively $M^i c$, no matrix powers) $+$ $O(n)$ vector dot products, and $O(n)$ storage. It involves no inverses, no factorization, and no extra space.

However, personally, I am not that familiar with good numerical methods to compute the characteristics polynomial of a matrix over $\mathbb{R}$.

Answer 3

So you want to compute $$H(s) = B(sI-A)^{-1}C$$ and $(sI-A)$ has a bad condition number. There are a few things to try.

First, you still want to avoid actually forming the inverse. You can use an iterative method like Krylov subspaces as the user J.D. mentioned above, but just because a method is iterative doesn't mean it will perform well when the matrix has a bad condition number.

In order for $H(s)$ to work out to a scalar, the right and left-hand matrices must really be vectors. That is, if $(sI-A)^{-1}$ has dimension $M$ by $N$, then $B$ must be $1$ by $M$ and C must be $N$ by $1$. So really, we're looking first for the solution $x$ to the problem $$(sI-A)x = C$$ with a badly scaled matrix on the left. After that, pre-multiplying by $B$ should be easy.

There are a couple of techniques you can use, but far and away the most popular is preconditioning. The basic idea is that you can use a cleverly chosen matrix to pre-multiply both sides of your problem which will reduce the condition number.

Your problem seems like it could benefit from a simple Jacobi preconditioner, but below that in the link there are other methods that make use of partial matrix factorizations.