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

Question

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

Answer 1

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.