Numerically stable computation of the Characteristic Polynomial of a matrix for Cayley-Hamilton Theorem

I was wondering if there is any known way to compute the Charactaristic Polynomial P of a matrix A numerically stable in the sense of realizing the Cayley-Hamilton Theorem, i.e. that P(A)=0.

I've found and implemented a couple of ways, but when the dimension of the matrix is even a few hunderds, P(A) is not even close to the zero matrix.

'Numerically Reliable Computation of Characteristic Polynomials' by Pradeep Misra, Enrique S. Quintana, Paul M. Van Dooren claims that reduction to Frobenius canonical form, use of Hyman's method for computing the determinant of a Hessenberg matrix, Faddeev-LeVerrier recursion, finding the polynomial by first computing the eigenvalues of the matrix are not proved to be numerically stable, and presents an algorithm which is numerically stable.

I implemented the algorithm on MATLAB:

function p = characteristic(A)

n = size(A,1);
[AA,BB,~,~] = hess(A, eye(n,n));
F = zeros(n,n); G = zeros(n,n);

F(:, 2:n) = AA(:,1:n-1); F(1,1) = -1;
fn = AA(:,n);

G(:, 2:n) = BB(:,1:n-1);
gn = BB(:,n);

res = zeros(n, n+1);
for i=0:n
if i == 0
res(:,n+1) = F\(-fn);
elseif i == 1
res(:,n) = F\(G*res(:,n+1) + gn);
else
res(:,n+1-i) = F\(G*res(:,n+1-i+1));
end
end

p = res(1,:);
for i=1:n-1
p = p*AA(i+1,i);
end
p = p*(-1)^(n-1);


which does work well for small dimenstions (P(A) is close to the zero matrix), but for larger dimenstions, it's better than MATLAB's 'poly' but still very far from the zero matrix.

>>A = rand(10,10); p = characteristic(A); polyvalm(p,A); ans(1,1)
ans =
-2.1472e-10
>>A = rand(100,100); p = characteristic(A); polyvalm(p,A); ans(1,1)
ans =
9.8824e+152

>>A = rand(20,20); p = characteristic(A); polyvalm(p,A); ans(1,1)
ans =
-581.7226
>>A = rand(20,20); p = poly(A); polyvalm(p,A); ans(1,1)
ans =
3.9789e+03


So I still didn't find an algorithm to calculate the characteristic polynomial in a way that P(A) will be close to the zeros matrix.

I'll much appreciate any answers, Niv

Generally speaking, $P(A)$ would not be close to the zero matrix: if you compute the polynomial $P$ to a relative accuracy of $\epsilon$ (which is $10^{-16}$ in double-precision), then $P(A)$ would be something on the order of $\epsilon \|A^n\|$ ($n$—the size of the matrix), which would be huge.

So you can only ask that $P(A)$ be as small as the typical perturbation due to roundoff error. That would be as good as an exact zero in floating-point arithmetic.

I don't think it matters much how you compute the characteristic polynomial itself, so the exact algorithm probably doesn't matter.

This is not specific to matrices either: even for a real roots of $p(x)=0$, you can't really ask that $p(x)$ be exactly zero in floating-point arithmetic. So to determine if $x$ is a root of $p$, you compare $p(x)$ not to zero, but to some error bound derived from the coefficients of $p$ and machine epsilon. There is a clear, advanced-level discussion of these kinds of issues in Higham's Accuracy and Stability of Numerical Algorithms (see Section 5.1 about Horner's method and its error bounds).

In addition, the roots of polynomials are notoriously sensitive to the coefficients of those polynomials expressed in the monomial basis.

Since $P(A)=0$ is an unreachable goal, you have to consider what you're trying to achieve, because whatever it is, it can probably be done without evaluating the coefficients of the characteristic polynomial explicitly or requiring that $P(A)=0$. If you absolutely must have $P(A)=0$, I'd recommend simply using arbitrary-precision floating-point arithmetic with number of digits proportional to the size of the matrix.

• "In addition, the roots of polynomials are notoriously sensitive to the coefficients of those polynomials expressed in the monomial basis." - The Perfidious Polynomial is a classic of the genre. – J. M. Dec 29 '16 at 19:29
• @Niv, even with that paper's claim, I would still not trust a method for producing coefficients in the monomial basis. Something fun: try to look at the condition number of F and F\G (MATLAB has rcond() IIRC), which is involved in all your backslash calls. – J. M. Dec 29 '16 at 19:34
• Thank you so much for this usefull answer @Kirill. P(A) should not be huge for example when A is orthogonal (A^n will be orthogonal too), but P(A) was still horrible in these cases. I do need the polynomials coefficients explicitly, because I want to solve matrix polynomials with large degrees (much larger than the dimension), and I want to compute a polynomial division with remainders. By dividing the polynomial by the characteristic polynomial, I'll be able to calculate a matrix polynomial with a much lower degree and achieve the same result. – Niv Hoffman Dec 29 '16 at 19:36
• As I see it now, @Niv, you've asked an XY problem. Consider asking a new question regarding your large degree matrix polynomials; I suspect that a clever linearization would be better than your current pursuit. – J. M. Dec 29 '16 at 20:27
• @NivHoffman Also, even for an orthogonal $A$, the coefficients of $P$ would get quite large, and $P(A)$ will be proportional to $\epsilon_{\mathrm{mach}}\|P\|\|A^n\|$, which should explain that. – Kirill Dec 30 '16 at 18:04