# Numerically stable explicit solution of small linear system

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.

• I know it comes very late, but here is my suggestion: as Gaussian elimination with total pivoting is known to be stable, it can make sense to hard-code the algorithm for the tiny sizes. Pivoting complicates matter as there are $(n!)^2$ ways to choose the successive pivots, leading to $(n!)^2$ different sets of formulas; you can reduce this complexity by swapping what needs to be swapped, reducing the number of cases to $1^2+2^2+\cdots n^2$. – Yves Daoust Sep 30 '15 at 12:30

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.

• The effort it takes me to implement it is about 15 minutes, because I simply enter a general 1x1, 2x2, 3x3 and 4x4 matrix into a CAS (Maple for me) and invert it. It shall return an explicit (C-like) result (supposedly based on Cramer's rule). Your second point is exactly my concern. In the result there will be higher order products of the matrix elements. Obviously this could introduce errors due to 'almost cancellation' of the different terms. But the question is, if it is possible to write the result in a form where this does not occur. Speed is not the main concern in this place. – highsciguy Mar 14 '13 at 11:13

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.

• When you write smooth here, do you mean that it is also smooth when evaluated with finite precision, i.e. stable (that is what I tried to say). See also my comment to GertVdE's answer. I think that I can exclude almost singular matrices (I suppose that in such cases the analysis of my physical problem must be reformulated). – highsciguy Mar 14 '13 at 11:21
• I mean "is at least twice continuously differentiable". I think think the matrix inverse map is infinitely continuously differentiable for all $A$ such that $\det(A) \neq 0$. – Geoff Oxberry Mar 14 '13 at 15:06
• Your comment about 'finite differencing used in an implicit ODE method' applies to me. Since the dimension $n$ of $A$ is much smaller than the dimension of my ODE system (this matrix arises merely in a mapping of a few variables), robustness is much more important at this stage than speed. In particular since in the development stage I will never know where encountered numerical errors arise if I don't make sure that individual components are safe. – highsciguy Mar 15 '13 at 11:36

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.

• Floating-point approximation (round off, cancellation, etc.) all count when it comes to stability. Even if you have a formula for the answer, you have to work out whether it can be calculated accurately in finite-precision arithmetic. – Bill Barth Mar 14 '13 at 1:33
• I do not see this answer as negative as others seem to see it. Of course the stability issue exists also for explicit results. But I believe that ctNGUYEN just wanted to say an approximate solution such as the expansion in a small quantity can actually be more precise than the full explicit result which, I think, is correct. In a sense I ask for explicit solutions which treat such difficult cases, such that the formula will always be stable. – highsciguy Mar 14 '13 at 11:30