Numerically stable explicit solution of small linear system

Question

I have an inhomogeneous linear system

$$ Ax=b $$

where $A$ is a real $n\times n$ matrix with $n\leq 4$. The nullspace of $A$ is guaranteed to be of zero dimension so the equation has a unique inverse $x=A^{-1} b$. Since the result enters the right hand side of an ODE, which I intend to solve using an adaptive method, it is important that the solution is smooth with respect to small variations of the elements of $A$ and $b$. Because of this requirement and the small dimensionality I thought to implement explicit formulas for $A^{-1} b$. The elements can be exactly zero or take very different values. My question is if this makes sense to you and if there are known stable expressions for this. I am coding in C for x86 systems.

Answer 1

Before implementing explicit formulas, I would ask myself the question: "is it worth it?":

• Is it worth to spend the time writing, debugging and validating these explicit formulas while you could easily link to BLAS+LAPACK that use the classical Gaussian elimination?
• Do you expect to gain stability? I don't think programming explicit formulas (like Cramer's rule) will give you better stability, on the contrary.
• Do you expect to gain speed? Did you already profile your whole program? What fraction of time is spent in solving hese 4x4 systems?
• What is the probability that in one year's time, you improve your model and you need 5 equations instead of 4?

My advice: use the BLAS/LAPACK combination first, see if it works, profile the whole program, ask a student to implement explicit formulas (sorry, being sarcastic here) and make a comparison on speed and robustness.

Answer 2

The only explicit inverse result I know of is Cramer's Rule, which has been recently shown to be computable in $\mathcal{O}(n^{3})$ time (like Gaussian elimination; unsure of the constant in front of the leading factor, though).

The matrix inverse of $A$ is a smooth function of $A$ so long as $\det(A) \neq 0$, and the solution $x$ is certainly a smooth function of $b$, so as long as the right-hand side of the ODE is a smooth function of $x$ and you avoid cases where $A$ is rank-deficient, I'd think your right-hand side would be smooth. (Here, I take smooth to mean "at least twice continuously differentiable".)

To be safe, it's probably best to make sure that $A$ is not numerically rank-deficient either (i.e., does not have small singular values).

The problem with Cramer's Rule is that its stability properties are unknown except for $n = 2$ (which is forward stable, but not backward stable). (See Accuracy and Stability of Numerical Algorithms, 2nd edition, by N. Higham.) It is not considered a reliable algorithm; Gaussian Elimination with Partial Pivoting (GEPP) is favored.

I would expect the problem with using BLAS+LAPACK to carry out GEPP in an ODE solve would be any finite differencing used in an implicit ODE method. I know that people have solved linear programs as part of a right-hand side evaluation, and because they did so naïvely (just plugged the linear program solve into the right hand side, calling a simplex algorithm), they greatly reduced the accuracy of their computed solution and substantially increased the time it took to solve the problem. A labmate of mine figured out how to solve such problems in a much more efficient, accurate manner; I'll have to look to see if his publication has been released yet. You may have a similar problem regardless of whether you opt to use GEPP or Cramer's Rule.

If there's any way you can calculate an analytical Jacobian matrix for your problem, you may wish to do that to save yourself some numerical headaches. It will be cheaper to evaluate, and probably more accurate, than a finite difference approximation. Expressions for the derivative of the matrix inverse can be found here if you need them. Evaluating the derivative of the matrix inverse looks like it would take at least two or three linear system solves, but they would all be with the same matrix and different right-hand sides, so it would not be considerably more expensive than a single linear system solve.

And if there's any way you can compare your computed solution to a solution with known parameter values, I would do that, so that you can diagnose whether you've encountered any of these numerical pitfalls.

Answer 3

Not sure that can help but i just think when you talk about stable solution, you're talking about approximation methods. When you compute things explicitle, stability doesn't have sens. That say you have to accept an approximative solution if you want to gain on stabilities.