% Stability of Least Squares algorithms
format long
disp('___---------------- Test 1: polynomial fit of a function')
disp('___m=100;')
disp('___n=15;')
disp('___t=(0:m-1)''/(m-1);')
disp('___A=[];')
disp('___for k=1:n')
disp('___    A = [A t.^(k-1)];  % Vandermont matrix')
disp('___end')
disp('___b = exp(sin(4*t));')
disp('___b = b/2006.787453080206;')

m=100;
n=15;
t=(0:m-1)'/(m-1);
A=[];
for k=1:n
    A = [A t.^(k-1)];  % Vandermont matrix
end
b = exp(sin(4*t));
b = b/2006.787453080206;
wait
disp(' ....')
disp('___Now, solve the system: ')
disp('___x1=A\b; y1=A*x1; ')
x1=A\b; y1=A*x1;
waitt('Done.')
disp('___Compute the condition number of A: kappa(A)=sigma_1/sigma_n')
kapa = cond(A)
disp('___Compute a measure of the closeness of the fit:')
disp('___theta=asin(norm(b-y1)/norm(b) )')
theta=asin(norm(b-y1)/norm(b) )
wait
disp('___eta=norm(A)*norm(x1)/norm(y1)')
eta=norm(A)*norm(x1)/norm(y1)
disp('___x1(15)')
x1(15),wait
disp('___Next we test the QR factorization ')
disp('___with Householder triangularization:'),wait
disp('___[Q,R]=qr(A,0);')
disp('___x2=R\(Q''*b); ')
disp('___x2(15)')
[Q,R]=qr(A,0);
x2=R\(Q'*b);
x2(15)
wait
disp('___Third test: QR factorization ')
disp('___with modif. Gram-Schmidt:'),wait
disp('___[QM,RM]=MGS_qr(A);')
disp('___x3=RM*(QM''*b);  %OBS! AR=Q')
disp('___x3(15)')
[QM,RM]=MGS_qr(A);
x3=RM*(QM'*b); 
x3(15),wait
disp('___One may wonder how good the factorization is: ')
disp('___norm(A*RM-QM)')
norm(A*RM-QM),wait
disp(' ')
disp('___Fourth test: the normal equation '),wait
disp('___x4=(A''*A)\(A''*b);')
disp('___x4(15)')
AA=A'*A;
x4=AA\(A'*b);
x4(15),wait

disp(' ...')
disp('___Fifth test: SVD '),wait
disp('___[U,S,V]=svd(A,0);  % reduced SVD')
disp('___x5=V*(S\(U''*b));')
disp('___x5(15)')
[U,S,V]=svd(A,0);  % reduced SVD
x5=V*(S\(U'*b));
x5(15),wait