Почему я получаю две разные обратные матрицы для одной и той же матрицы N * N при увеличении N в Matlab?

Question

В качестве эксперимента я просто хотел увидеть время вычислений между теорией Кэли-Гамильтона и MATLAB inv() функция. Я знал, что CH будет медленнее на CPU из-за количества матричных продуктов, однако я не ожидал, что они дадут мне разные ответы при увеличении N.

Для квадратных матриц меньше 30 * 30 обратные значения примерно одинаковы. Но после этого они начинают резко отличаться друг от друга. К тому времени, когда N = 100, они вообще не имеют ничего общего.

Это проблема числовых вычислений, или здесь что-то еще происходит? И чему я могу доверять? Я предполагаю, что inv()высоко оптимизирован и заслуживает доверия, но было бы неплохо получить некоторую информацию от других.

Вот код:

n = 100;
A = randn(n);

% MATLAB inv()

tic;
initime = cputime;
time1   = clock;
A_inv = inv(A);
fintime = cputime;
elapsed = toc;
time2   = clock;
fprintf('TIC TOC: %g\n', elapsed);
fprintf('CPUTIME: %g\n', fintime - initime);
fprintf('CLOCK:   %g\n', etime(time2, time1));

% Cayley-Hamilton inversion

tic;
initime = cputime;
time1   = clock;
p_coeff = poly(A);
A_inv_2 = 0;
for ii = 1:n-1
    A_inv_2 = A^(ii)*p_coeff(end-1-ii) + A_inv_2;
end
A_inv_2 = 1/-p_coeff(end) * (A_inv_2 + eye(n)*p_coeff(end-1));    

fintime = cputime;
elapsed = toc;
time2   = clock;
fprintf('TIC TOC: %g\n', elapsed);
fprintf('CPUTIME: %g\n', fintime - initime);
fprintf('CLOCK:   %g\n', etime(time2, time1));

Спасибо всем, кто нашел время, чтобы ответить.

Answer 1

Метод Кэли-Гамильтона является очень нестабильным методом для вычисления инверсий, потому что он предполагает поднятие матриц до больших степеней.

Рассмотрим матрицу, которая может быть диагонали в A=inv(P)DP где D - диагональная матрица. Когда поднят до 100-й степени, это становится A^100 = inv(P) D^100 P, Любая разница в размере между диагональными элементами в D будет взорвана этой операцией. Например, рассмотрим разницу между 2^100 и 0,5^100.

Это действительно легко увидеть в вашей программе Matlab. Распечатайте A * A_inv и A * A_inv_2. Первый очень близок к идентичности, а второй содержит чепуху:

A*A_inv_2
ans = 1.0e10 *
  0.2278  0.3500 -0.2564 ...