Why do I get two different inverse matrices for the same N * N matrix as N increases in Matlab?

Question

As an experiment I just wanted to see the computational times between Cayley-Hamilton theory and the MATLAB inv() function. I knew that C-H would be slower on a CPU due to the amount of matrix products, however I wasn't expecting them to give me different answers as N increases.

For square matrices less than around 30 * 30, the inverses are about the same. But after this point they begin to differ from each other quite drastically. By the time N = 100 they share no similarities at all.

Is this a numerical computation issue, or is there something else happening here? Also which can I trust? I'm assuming inv()is highly optimised and trust worthy, but it would be nice to have some input from others.

Here is the code:

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

Thank you to anyone who takes the time to answer.

Answer 1

The Cayley-Hamilton method is a very unstable method for computing inverses because it involves raising matrices to high powers.

Consider a matrix that can be diagonalized into A=inv(P)DP where D is a diagonal matrix. When raised to the 100th power, this becomes A^100 = inv(P) D^100 P. Any difference in size between the diagonal entries in D will be blown up by this operation. For example, consider the difference between 2^100 and 0.5^100.

It is actually easy to see this within your Matlab program. Print out A * A_inv and A * A_inv_2. The first is very close to the identity, while the second contains nonsense:

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