% LSQGeom Script file demonstrates the geometry
%         of a linear least squares problem.
%
clc
close all
% 
% Table of Data for least squares fit
% 
x = [1:3]'; y = [2; 3; 8];
% 
% Get coefficients; plot curve and data points
% 
c = polyfit(x,y,1);
p = polyval(c,x);
plot(x,y,'ro',x,p,'b-')
xlabel('x'), ylabel('y')
title('Degree-one polynomial fitting three data points')
input('Press a key to continue');
% 
% Now re-conceive this problem as finding the nearest 
% point to y in the plane spanned by x and ones(size(x)).
% 
A = [x ones(size(x))]; b = y;
cqr = A \ b;
% 
% The coefficients are the same!
% 
c = c';
norm(c-cqr)
% 
% Here is the predicted value of y:
% 
yhat = A*c;
% 
% Observe that the residual y - A*c = yhat -y
% is perpendicular to the plane spanned by
% the columns of A:
A'*(yhat - y)
% 
% We'll plot this situation in 3-D
% 
figure
n = 8;
x1 = linspace(-2,5,n)';
x2 = linspace(-1,6,n)';
x3 = -ones(n,1)*x1' + 2*x2*ones(1,n);
% 
% Above is the equation of the plane
% 
%     x3 = - x1 + 2*x2
%     
% spanned by the vectors [1 1 1]' and [1 2 3]'.
% [Exercise!]
% 
mesh(x1,x2,x3)
hold on
% Plot the vectors that span the plane
plot3([0 1],[0 1],[0 1],'c-') 
plot3([0 1],[0 2],[0,3],'c-')
plot3(1,1,1,'b*',         ...
      x(1),x(2),x(3),'b*',...
      0,0,0,'bo')
% Plot the residual vector and its endpoints,
% which are y and the least squares approx.
plot3([y(1) c(1)*x(1)+c(2)], ...
      [y(2) c(1)*x(2)+c(2)], ...
	  [y(3) c(1)*x(3)+c(2)],'r-')
plot3(y(1),y(2),y(3),'r*', ...
      c(1)*x(1)+c(2),c(1)*x(2)+c(2),c(1)*x(3)+c(2),'r*')
xlabel('x1'), ylabel('x2'), zlabel('x3')
axis('equal')
title(['Residual vector is ' '{\perp}' ' span of columns of A'])
view(110,30)
text(2,3,8.5,'y')
text(1,2,3.5,'x')
text(1,1,1.5,'e')
text(1,4,7.5,'y*')