% Lecture 11.3

clear;
close all;

m = 50; % constant set by problem
n = 12; % constant set by problem

t = linspace(0,1,m); % linearly spaced with 50 grid points
b = cos(4*t'); % function that we are trying to fit with polynomial
               % of degree n - 1;

	     
% creating Vandermonde matrix A associated
% with fitting a least square fitting on 
% this grid with a polynomial of degre n-1

A = vander(t); % creates a vandermonde matrix whose columns
               % are powers of t
A = fliplr(A); % flip matrix
A = A(:,1:n);  % fit by polynomial of degree n-1

% we now solve the least square coefficient vector via
% six methods - sort of...

x = zeros(n,6);

% (a) Formation and solution of the normal equations using \ 
  B = A'*A;                   % see ALG 11.1
  R = chol(B);                % form cholesky factorization
%  w1 = forwardsubs(R',A'*b); % my forward subs function
%  x(:,1) = backsubs(R,w);    % my backsubs function
  x(:,1) = R\(R'\(A'*b));     % Matlab's \ operator

  clear R B w;            % free memory

% (b) QR factorization computed by modified Gram-Schmidt
  [Q,R] = mgs(A);
  x(:,2) = backsubs(R,Q'*b);
  clear Q R;
  
% (c) QR factorization computed by householder triangulation
  [W,R] = house(A);
  Qb = formqb(W,b);
  x(:,3) = backsubs(R,Qb);
  clear W R;
  
% (d) use MATLAB's QR factorization in 
  [Q,R] = qr(A);
  x(:,4) = backsubs(R,Q'*b);
  clear Q R;

% (e) formation and solution using MATLAB's \
  x(:,5) = A\b; 

% (f) using MATLAB's SVD
  [U,S,V] = svd(A,0);   % note, A = USV', reduced svd
  c = U'*b;             % form U'b
  for i = 1:n       
    c(i) = c(i)/S(i,i); % divide through by singular values
  end 
  x(:,6) = V*c;         % form solution - see ALG 11.3
  clear U S V;