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