% Computer problem 3.5 on p.153 Heath
% Least squares problem.
% K. Ming Leung, 04/26/03

clear all; hold off; xPlot=0:0.01:1.02;
x=[1.02 0.95 0.87 0.77 0.67 0.56 0.44 0.30 0.16 0.01]';
y=[0.39 0.32 0.27 0.22 0.18 0.15 0.13 0.12 0.13 0.15]';
plot(x,y,'ok'); axis([0 1.2 0.1 0.4]); hold on; % plot data
m=max(size(x)); n=5;
A=[y.*y x.*y x y ones(m,1)], cond(A), % matrix A
b=x.*x;
AT=A';
ATA=AT*A;
pNormal=ATA\AT*b,
pMatlab=A\b,
[Q,R]=qr(A); RS=R(1:n,1:n); pqr=backsubstitution(RS,Q'*b,n),
ellipsePlot(pMatlab,xPlot,'m-');

% Introduce random noise
L=-0.015; H=-L;
xp=x+L+rand(m,1)*(H-L); yp=y+L+rand(m,1)*(H-L);
plot(xp,yp,'xk');
Ap=[yp.*yp xp.*yp xp yp ones(m,1)]; cond(Ap),
bp=xp.*xp;
pp=Ap\bp,
ellipsePlot(pp,xPlot,'r-');

nn=n-1;
[U, S, V]=svd(A),
psvd=zeros(n,1);
for i=1:nn, psvd=psvd+(U(:,i)'*b/S(i,i))*V(:,i); end
psvd, ellipsePlot(psvd,xPlot,'b--');

[Up, Sp, Vp]=svd(Ap);
ppsvd=zeros(n,1);
for i=1:nn, ppsvd=ppsvd+(Up(:,i)'*bp/Sp(i,i))*Vp(:,i); end
ppsvd, ellipsePlot(ppsvd,xPlot,'g--');