计算科学 - 凯莱-汉密尔顿定理矩阵特征多项式的数值稳定计算 - 吾爱随笔录

我想知道是否有任何已知的方法来计算矩阵 A 的特征多项式 P 在实现 Cayley-Hamilton 定理的意义上数值稳定，即 P(A)=0。

我已经找到并实现了几种方法，但是当矩阵的维数甚至是几百倍时，P(A) 甚至不接近零矩阵。

Pradeep Misra、Enrique S. Quintana、Paul M. Van Dooren 的“特征多项式的数值可靠计算”声称，简化为 Frobenius 规范形式，使用 Hyman 的方法计算 Hessenberg 矩阵的行列式，Faddeev-LeVerrier 递归，找到多项式通过先计算矩阵的特征值并不能证明是数值稳定的，并提出了一种数值稳定的算法。

我在 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);

这对于小维度（P(A) 接近零矩阵）效果很好，但对于更大的维度，它比 MATLAB 的“多边形”要好，但离零矩阵仍然很远。

>>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

所以我仍然没有找到一种算法来计算特征多项式，使 P(A) 接近零矩阵。

我会非常感谢任何答案，尼夫