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

В качестве эксперимента я просто хотел увидеть время вычислений между теорией Кэли-Гамильтона и функцией MATLAB inv(). Я знал, что C-H будет медленнее на процессоре из-за количества матричных произведений, однако я не ожидал, что они будут давать мне разные ответы по мере увеличения 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));

Спасибо всем, кто найдет время ответить.


matlab algorithm matrix inversion
person Bradley    schedule 18.04.2017    source источник
comment
Если я правильно помню, вычислить характеристический полином численно устойчивым способом сложно. Поэтому я подозреваю, что виновата poly(). Однако это всего лишь предположение. Я голосую в надежде, что кто-то более знающий, чем я, ответит на этот вопрос.   -  person Ash    schedule 18.04.2017
comment
Ставлю на числовые вопросы в качестве объяснения. Я помню книгу Noble и Daniel's Applied Linear Algebra, всегда настаивавшую на том, что вы не должны реализовывать свои собственные обратные маршруты, потому что они обязательно будут страдать от числовых проблем. Чтобы избежать этих проблем, нужно сделать много оптимизаций, и inv, несомненно, имеет эти оптимизации. Насчет того, кому можно доверять: inv, конечно. Сравните imagesc(A*A_inv) и imagesc(A*A_inv_2) в своем примере.   -  person Luis Mendo    schedule 18.04.2017

Ответы (1)


arrow_upward
1
arrow_downward

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

Рассмотрим матрицу, которую можно преобразовать в 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 ...
person Peter de Rivaz    schedule 18.04.2017