function ShowLSfit()
% Polynomial LS fits of 
%    y = .05*sin(1000*t) +.5*cos(pi*t) -.4*sin(10*t)

% Note that if A = QR then A(:,1:k) = Q(:,1:k)R(1:k,1:k)
% This means that we can "live" off of single QR factorization..
t = linspace(0,1,101)';
At = [ones(101,1) t  t.^2 t.^3 t.^4  t.^5  t.^6];
[Q,R] = qr(At,0);
fvec = f(t);
ftilde = Q'*fvec;
% Note that At(:,1:k+1) = Q(:,1:k+1)R(1:k+1,1:k+1) is the thin QR
% factorization that we need to determine p_k. Moreover,
% if we go on to solve R(1:k+1,1:k+1)*alfa = ftilde(1:k+1), then
%   p_k(z) = alfa(1) + alfa(2)z + alfa(3)z^2 +...+alfa(k+1)z^k

% Some precomputation in anticipation of the required plotting...
z = linspace(0,1,500)';
Az = [ones(500,1) z z.^2 z.^3 z.^4  z.^5  z.^6];
fz = f(z);
close all
for k=1:6
    % Generate and display the best polynomial fit of the
    % data assuming that the allowable degree is k..
    figure
    alfa = R(1:k+1,1:k+1)\ftilde(1:k+1);
    plot(z,fz)
    hold on
    plot(z,Az(:,1:k+1)*alfa,'r','linewidth',3)
    title(sprintf('Degree = %1d',k),'fontsize',14)
    shg
    pause(2)
end

function y = f(t)
y = .05*sin(1000*t) +.5*cos(pi*t) -.4*sin(10*t);

% -2 for doing six QR factorizations
% -1 for uninformative plots
% -1 for least squares on the larger mesh
% -1 for doing QR on A'*A instead of A