%% Demonstrating good and bad ways to invert matrices
clear;
format long; 

% create some random data
K = 1000;   N = K+1;   b = ones(K,1);
x = randn(N,K); y = x*b + randn(N,1);
xy = x'*y;
xx = x'*x;

bhat = zeros(K,4);
% from slow to fast ...
tic; bhat(:,1) = (xx)^(-1)*xy; toc;
tic; bhat(:,2) = inv(xx)*xy; toc;
tic; bhat(:,3) = xx \ xy;    toc;
% the qr-way: slow but can be more accurate
tic;
[q r] = qr(x,0);
s = inv(r);
ixx = s*s';
bhat(:,4) = ixx*xy;;
toc;
sum(abs(xx*bhat-xy*ones(1,4))) % notice that \ gives smallest error
clear bhat;

%% Solving multiple equations 
M = 50; 
y = (x*b)*ones(1,M) + randn(N,M);
xy = x'*y;
tic; bhat{1} = inv(xx)*xy; toc;
tic; bhat{2} = xx*xy; toc;
% In a loop for illustrative purposes
tic;
ixx = inv(xx);
for m=1:size(y,2)
  % ... might be doing stuff that changes xy ... 
  bhat{1}(:,m) = ixx*xy(:,m); %#ok
end
toc;
tic; % for a general matrix, xx, \ uses LU decomp + backsub to solve, we
     % do that explicitly here
[L, U, p] = lu(xx,'vector');
lower.LT=true;
upper.UT=true;
for m=1:size(y,2);
  bhat{2}(:,m) = linsolve(U,linsolve(L,xy(p,m),lower),upper);
end
toc; 
clear lower;
tic; % for symmetrix matrix with positive diagonals (like we have), \ uses
     % cholesky decomposition and back substitution to solve, we do that
     % explicitly here
  R = chol(xx);
  lower.UT=true;
  lower.TRANSA=true;
  upper.UT=true;
  for m=1:size(y,2);
    bhat{3}(:,m) = linsolve(R,linsolve(R,xy(p,m),lower),upper);
  end
toc; 

%% Illustration of usefulness of QR 
% make inv and \ fail because of inaccuracy
d = 1.e-8;
x = [1 1; d 0; 0 d]
y = [1;2;3];
xx = x'*x;
xy = x'*y;
bbad1 = xx^(-1)*xy
bbad2 = inv(xx)*xy
bbad3 = xx \ xy
[q r] = qr(x,0);
s = inv(r);
ixx = s*s';
bbad4 = ixx*xy
%% Explanation
% The problem with using inv() and \ is that they require first computing
% x'*x, and computing x'*x involves adding a relatively big number, 1 to
% a tiny number, 1e-16.  Because of rounding, 1 + 1e-16 = 1, and x'*x
% becomes singular.  
xx
%% 
% Making d slightly bigger reduces the problem
d = 2.e-8;
x = [1 1; d 0; 0 d]
y = [1;2;3];
xx = x'*x
xy = x'*y;
bbad1 = xx^(-1)*xy
bbad2 = inv(xx)*xy
bbad3 = xx \ xy
[q r] = qr(x,0);
s = inv(r);
ixx = s*s';
bbad4 = ixx*xy

