clear all
colr='grkxgrkwgrbkwgrbkw';
disp('+++++++++++++++++++++++++++++++++++++++++++++++++++')
disp('++++++++++welcome to the least squares demo++++++++')
disp('+++++++++++++++++++++++++++++++++++++++++++++++++++')

%in this demo, we try to experiment with least squares
%first we generate a matrix with more row than columns

disp('A = [0.1389    0.2722    0.4451')
disp('     0.2028    0.1988    0.9318')
disp('     0.1987    0.0153    0.4660')
disp('     0.6038    0.7468    0.4186];')


A = [0.1389    0.2722    0.4451
    0.2028    0.1988    0.9318
    0.1987    0.0153    0.4660
    0.6038    0.7468    0.4186];

null=input(' ');

%now we select some vector `x0', and select the right hand side `b0' in 
%a way that x0 is the solution to Ax=b0. remember that in general there is
%no exact solution to the system Ax=b

disp('x0=[1:3]''')
x0=[1:3]'

null=input(' ');

disp('b0=A*x0')
b0=A*x0

null=input(' ');

%now, we change a bit the vector b0 by adding a `noise' vector.
%the `noise' a simply a vector with small entries

disp('noise=[0.0846; 0.0525; 0.0203; 0.0672]')

noise=[0.0846
       0.0525
       0.0203
       0.0672]

null=input(' ');

disp('b=b1+noise')
b=b0+noise
% right now, we should not expect to have an exact solution to Ax=b
%let's compute now the least squares solution

null=input(' ');

disp('xleast=A\b')
xleast=A\b

null=input(' ');
%let's now compare the error in the least square solution,
%to that of the original x0 (remember that x0 was the exact
%solution to a `near-by' problem

disp('norm(A*xleast-b,2), norm(A*x0-b,2)')
norm(A*xleast-b,2), norm(A*x0-b,2)

null=input(' ');
%now, let's pick a random x and try its performance
disp('xrand=rand(3,1)')
xrand=rand(3,1)

null=input(' ');
disp('norm(A*xrand-b,2)')
norm(A*xrand-b,2)

null=input(' ');
%the least squares solution is the best performer in the 2-norm
%let's see how well it fdoes in other norms.
disp('norm(A*xleast-b,1), norm(A*x0-b,1), norm(A*xrand-b,1)')
norm(A*xleast-b,1), norm(A*x0-b,1), norm(A*xrand-b,1)

null=input(' ');
disp('norm(A*xleast-b,inf), norm(A*x0-b,inf), norm(A*xrand-b,inf)')
norm(A*xleast-b,inf), norm(A*x0-b,inf), norm(A*xrand-b,inf)

null=input(' ');
%but, if I work hard enough, I may be able to find, in the other norm,
%a better performer than the least solution

disp('xtry=[-.01 -.01 .01]''+xleast')
xtry=[-.01 -.01 .01]'+xleast

null=input(' ');
disp('norm(A*xtry-b,inf) norm(A*xleast-b,inf)')
norm(A*xtry-b,inf),  norm(A*xleast-b,inf)


disp('+++++++++++++++++++++++++++++++++++++++++++++++++++')
disp('++++++++++END OF PART I OF THE DEMO++++++++++++++++')
disp('+++++++++++++++++++++++++++++++++++++++++++++++++++')
