Eigenvalues of Small Matrices

I'm writing a small numerical library for 2x2, 3x3, and 4x4 matrices (real, unsymmetric).

A lot of numerical analysis texts highly recommend against computing the roots of the characteristic polynomial and recommend using the double-shifted QR algorithm. However, the size of the matrices makes me question whether or not simply computing the characteristic polynomial and finding the roots can be sufficient for this. I found this answer on StackExchange that supports this, but I know that errors in the polynomial coefficients can produce inaccurate zeros of the polynomial (and thus different eigenvalues). On the other hand, the equation is only quartic at best and we have analytical formulas for the polynomial roots so we shouldn't get too far off.

What are the pros and cons to using the characteristic polynomial for getting eigenvalues this specifically for this case?

The first thing to note is that the correspondence between finding roots of a polynomial (any polynomial) and finding the eigenvalues of an arbitrary matrix is really direct, and it's a rich subject, see Pseudozeros of polynomials and pseudospectra of companion matrices by Toh and Trefethen and the references there.

Basically, the 2×2 case is trivial and the standard formula, $$x_1 = \frac{-b-\mathrm{sign}(b)\sqrt{Δ}}{2a}, \qquad x_2 = c/(ax_1), \qquad Δ = \det \begin{pmatrix} b & 2a\\ 2c & b \end{pmatrix}$$ is numerically stable and accurate, so long as the determinant $Δ$ is evaluated accurately—the direct formula will be inaccurate when $b^2≈4ac$, but there is an accurate formula due to Kahan that uses FMAs (https://hal.inria.fr/ensl-00649347v1/document).

There is no such direct equivalent even for cubic polynomials. (Edit: see the link to Kahan's method below in CADJunkie's comment. This might well be wrong.) The direct formula is not always numerically stable (contrary to your assumption, I think), and it can't be made numerically stable the same way as the quadratic formula by inserting the right signs somewhere. You could try evaluating it with extra precision, e.g. with double-native arithmetic. But the approaches that do work directly on the polynomial are pretty complicated. For example (https://doi.org/10.1145/2699468, which also works on quartic polynomials) you could use Newton's method with a good precomputed first guess, but it gets really rather complicated, and the speedup is not even all that large.

The explicit formulas for degree-4 polynomials are likewise not always numerically stable. The polynomials that are the hardest tend to have uncommon roots (small, or close to each other, differing in magnitude by a lot), but even testing your code on a few billion purely random polynomials can usually reveal numerical errors.

One curious thing related to this is that Jenkins-Traub, which is a good common way of finding polynomial roots, is actually an eigenvalue algorithm (inverse iteration) in disguise.

I would say the explicitness of the formulas is, in a way, misleading: it tricks you into thinking that because the formula has a closed form that means it's somehow cheaper/easier. I would really recommend that you actually test/benchmark this on some test data. It doesn't have to be true: determining the roots of degree-$≥3$ polynomials is within a small integer factor of difficulty of the full problem of determining the eigenvalues of a small matrix, and the standard library routines for eigenvalues are much more robust and well-tested. So by reducing the small eigenvalue problem to a low-degree polynomial roots problem, you aren't necessarily simplifying it.

What are the pros and cons to using the characteristic polynomial for getting eigenvalues this specifically for this case?

I think the main con is that this assumption that you make:

On the other hand, the equation is only quartic at best and we have analytical formulas for the polynomial roots so we shouldn't get too far off.

that because a formula is in closed form, and analytical, that means it's easy/cheap/accurate, is not necessarily true. It could be true on specific data you might have, but as far as I know it's not true in general.

P.S. The whole distinction between in-closed-form and not-in-closed-form gets really finicky with computer arithmetic: you might think that $\cos (\cdots)$ in a cubic formula is a closed form, but as far as computer arithmetic is concerned, that's just another approximate rational function, it is maybe faster, but not fundamentally different, from the approximate rational function that defines the result of an eigenvalue algorithm.

• Thank you for the great answer. The bit about eigenvalues and poly roots being connected was new to me. I see that closed-form solutions aren't necessarily better when dealing with computer arithmetic. I was planning on using a quadratic and cubic poly root solver Kahan has written about (people.eecs.berkeley.edu/~wkahan/Math128/Cubic.pdf) and converting the quartic to a cubic using Descartes factorization for my solutions. Would you recommend implementing the QR algorithm or my stated approach? Sep 25 '17 at 14:55
• @CADJunkie That's cool, thank you, I didn't know Kahan wrote that, I'll read that later. It's difficult to recommend like this, rather the better way of settling these questions is to implement the ideas, then test and benchmark them and go with those benchmarks, that's much more definitive. Performance is very difficult to predict like this, before seeing results. But at least make sure to compare with all the standard library eigenvalue solvers. Sep 25 '17 at 16:04
• @CADJunkie Based on what Kahan wrote, the cubic case might be better done directly too. Assuming that, I think what I wrote is okay for quartics and higher. Sep 25 '17 at 18:18

Using a QR algorithm is the better way. I think it is best to use an algorithm best-suited for the task at hand.

In fact, even if you were trying to compute the roots of a polynomial, without intending to use them as the eigenvalues of a matrix, it is often recommended to create the Companion Matrix for that polynomial, and solve for the eigenvalues of that matrix. (i.e., do the opposite process to what you are considering.) The process is extremely robust, but not very computationally efficient. An algorithm specific to the task of finding polynomial roots (e.g., Jenkins-Traub, Method of Laguerre, etc.) would be more efficient. And even if you use one of these methods, there may still be some cases for which forming the Companion Matrix and computing its eigenvalues still gives better results.

Also, as Kirill indicated, there is no way to take advantage of the closed form solutions for a 3rd or 4th degree polynomial. I looked into this some years ago before writing a translation of the Jenkins-Traub algorithm. For numerical results, it is still best to write an algorithm from the ground up as a discrete solver and ignore the closed form solutions.