close all; clc; clear; 

% Problem EPage 45
% 
A = hilb(5)
AInv = inv(A);
ATAInv = inv(A' * A);

AInv
ATAInv

for n=3:6
  A = hilb(n);
  AInv_approx = inv(A);
  AInv_exact = invhilb(n);
  
  % Compute various approximations of (A^T A)^{-1}: 
  %
  R = AInv_exact * AInv_exact'; % the exact value of this inverse
  P = inv( A' * A ); % multiply A together first and then take the numerical inverse 
  Q = AInv_approx * AInv_approx'; % take the numerical inverse of A and then multiply together 
  
  disp( sprintf( 'n=%5d, ||P-R||= %10.6e, ||Q-R||= %10.6e', n, norm( P-R ), norm( Q-R ) ) );
end