% This script illustrates the stability of various
% least-squares algorithms
clear all
clc
m = 100; n = 15;
disp(sprintf('We use a degree %d polynomial to fit the function 1/(1+x^2)', n-1))
disp(sprintf('at %d equally spaced points in [-5,5].\n', m))
scrnsz = get(0, 'Screensize');
xi = linspace(-5,5,2*m+1);
eta = 1./(xi.^2+1);
figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
PolynomialLSQ(m,n-1);
clc
t = linspace(-5,5,m)';
A = [];
for i=1:n
    A = [A t.^(i-1)];
end
b = 1./(1+t.^2);
% normalized so that x(15) below should be 1
b = b/-.39345237595168944524e-7;
% Solve LSQ problem and compute condition number
x = A \ b;
y = A*x;
kappa = cond(A);
theta = asin(norm(b-y)/norm(b));
eta = norm(A)*norm(x)/norm(y);
disp(sprintf('cond(A) = %e', kappa))
disp(sprintf('A large condition number shows that the problem (matrix)'))
disp(sprintf('is ill-conditioned.\n'))
disp(sprintf('theta = %e', theta))
disp(sprintf('A small value of theta shows that it is possible'))
disp(sprintf('to fit the data well with a degree %d polynomial.\n', n-1))
disp(sprintf('eta = %e', eta))
disp(sprintf('eta can be anywhere between 1 and cond(A)\n'))
disp(sprintf('Condition of least squares problem'))
illcond = kappa + kappa^2*tan(theta)/eta;
disp(sprintf('kappa + kappa^2*tan(theta)/eta = %e', illcond))
loss = floor(log10(illcond));
disp(sprintf('We can expect to lose %d digits due to ill-conditioning', loss))
pause
clc
disp(sprintf('LSQ solution with degree %d polynomial.',n-1))
disp(sprintf('The problem is normalized so that the coefficient'))
disp(sprintf('computed below should be 1.\n'))
% Check Householder QR
[Q, R] = qr(A, 0);
x = R \ (Q'*b);
disp(sprintf('via Householder QR: c(%d) = %16.14f', n, x(15)))
disp(sprintf('                    error: %e', abs(1-x(15))))
loss_err = -floor(log10(abs(1-x(15))));
if loss_err < loss-2    % -2 as small tolerance
    disp(sprintf('Unstable algorithm\n'))
else
    disp(sprintf('Stable algorithm\n'))
end
pause
% Check modified Gram-Schmidt
[Q, R] = mgs(A);
x = R \ (Q'*b);
disp(sprintf('via modified Gram-Schmidt: c(%d) = %16.14f', n, x(15)))
disp(sprintf('                    error: %e', abs(1-x(15))))
loss_err = -floor(log10(abs(1-x(15))));
if loss_err < loss-2    % -2 as small tolerance
    disp(sprintf('Unstable algorithm\n'))
else
    disp(sprintf('Stable algorithm\n'))
end
pause
% Check stabilized MGS
[Q2, R2] = mgs([A b]);
Qb = R2(1:n,n+1);
R2 = R2(1:n,1:n);
x = R2 \ Qb;
disp(sprintf('via stabilized MGS: c(%d) = %16.14f', n, x(15)))
disp(sprintf('                    error: %e', abs(1-x(15))))
loss_err = -floor(log10(abs(1-x(15))));
if loss_err < loss-2    % -2 as small tolerance
    disp(sprintf('Unstable algorithm\n'))
else
    disp(sprintf('Stable algorithm\n'))
end
pause
% Check Normal equations
disp(sprintf('via normal equations:'))
x = (A'*A) \ (A'*b);
disp(sprintf('c(%d) = %16.14f', n, x(15)))
disp(sprintf('                    error: %e', abs(1-x(15))))
loss_err = -floor(log10(abs(1-x(15))));
if loss_err < loss-2    % -2 as small tolerance
    disp(sprintf('Unstable algorithm\n'))
else
    disp(sprintf('Stable algorithm\n'))
end
pause
% Check SVD
[U, S, V] = svd(A,0);
x = V*(S \ (U'*b));
disp(sprintf('via SVD: c(%d) = %16.14f', n, x(15)))
disp(sprintf('                    error: %e', abs(1-x(15))))
loss_err = -floor(log10(abs(1-x(15))));
if loss_err < loss-2    % -2 as small tolerance
    disp(sprintf('Unstable algorithm'))
else
    disp(sprintf('Stable algorithm'))
end
