clear
format long
%well-conditioned example 4.14
%A=[1 -4;2 3;2 2]; y=[-3;15;9];
%ill-conditioned example 4.15
x=linspace(2,4,11)'; %11 points in [2,4]
y=1+x+x.^2+x.^3+x.^4+x.^5+x.^6+x.^7;
A=[x.^0 x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7]; %VanderMonde matrix
%page 233, ex 2, 4-by-3 matrix A
%s=1e-10; A=[1 1 1;s s/2 s/4;s/2 s/3 s/4;s/3 s/4 s/5]; y=A*ones(3,1);
%A=hilb(7); y=A*ones(7,1);
[m,n]=size(A);
condAA=cond(A'*A) %very large condition number
%fit a degree 7 ploynomial
c1=(A'*A)\(A'*y);
disp('normal equation method')

%
[Q,R]=clgs(A);
err_clgs=norm(A-Q*R);
loss_clgs=norm(Q'*Q-eye(n)); %loss of orthogonality
b=Q'*y;
c2=R\b;
disp('QR: Classical Gram-Schmidt algorithm')

[Q,R]=mgs(A);
err_mgs=norm(A-Q*R);
loss_mgs=norm(Q'*Q-eye(n));
b=Q'*y;
c3=R\b;

disp('QR: Modified Gram-Schmidt algorithm')

%QR method based on Householder algorithm (page 231)
[Q,R]=houseqr(A); %produce a full QR factorization
err_house=norm(A-Q*R);
loss_house=norm(Q(:,1:n)'*Q(:,1:n)-eye(n));
b=Q(:,1:n)'*y;
c4=R(1:n,1:n)\b(1:n);
disp('QR: Householder algorithm, best orthogonality')

%QR method of MATLAB
[Q,R]=qr(A,0); %produce a reduced QR factorization
err_qr=norm(A-Q*R);
loss_qr=norm(Q'*Q-eye(n));
b=Q'*y;
c5=R(1:n,1:n)\b(1:n);
disp('QR: from MATLAB')

[Q,R]=mgs([A y]); %use Augmented matrix for better accuracy
%norm(Q'*Q-eye(n+1))
c6=R(:,1:end-1)\R(:,end);
disp('QR: Modified Gram-Schmidt algorithm (Augmented)')

fprintf('===========Compare Loss of orthogonality:=============\n')
fprintf('\t\t clgs \t\t mgs \t\t householder \t qr(MATLAB) \n')
fprintf('||A-QR||=0?\t%1.2e \t%1.2e\t%1.2e \t%1.2e \n',err_clgs,err_mgs,err_house,err_qr)

fprintf('||Q''Q-I||=0?\t%1.2e \t%1.2e\t%1.2e \t%1.2e \n',loss_clgs,loss_mgs,loss_house,loss_qr)

fprintf('===========Compare Obtained solutions:===============\n')
fprintf('normal\t\t clgs \t\t mgs \t householder \t\t qr\t mgs(augmented) \n')
fprintf('%1.2e \t%1.2e \t%1.3f\t%1.8f\t%1.8f \t%1.8f\n',[c1 c2 c3 c4 c5 c6]')