% ----------------------------------------
%   Solution of least squares  problem  min_x || Ax - y ||_2
%   using the method of normal equations.
%   Matrix  A is constructed as a Vandermonde matrix.
%   Program performs fitting to the function y = 5x+7
% --------------------------------------------------------------

d=15;  % degree of the polynomial
m=15;%number of discretization points or rows in the matrix A

x=zeros(1,m);
y=zeros(1,m);
b=zeros(1,m);

x = linspace(0,10.0,m); 

A=[];
brus = 0.1;
brus = brus*60.0*randn(size(x));

for i=1:1:m
	
  %  exact function which should be approximated
b(i) = 5*x(i) + 7;

% genererar observationer  med random brus

y(i) =  b(i) + brus(i);

end

% construction of a Vandermonde matrix

for i=1:1:m
  for j=1:1:d+1
    A(i,j)=power(x(i),j-1);
  end
end

% computing the right hand side in the method of normal equations
c=A'*y';

% computing matrix in the left hand side in the method of normal equations
C=A'*A;

l=zeros(d+1);

% solution of the normal equation using Cholesky decomposition

for j=1:1:d+1
  s1=0;
  for k=1:1:j-1
    s1=s1+l(j,k)*l(j,k);
  end
  l(j,j)=(C(j,j)-s1)^(1/2);
  for i=j+1:1:d+1
    s2=0;
    for k=1:1:j-1
      s2=s2+l(i,k)*l(j,k);
    end
    l(i,j)=(C(i,j)-s2)/l(j,j);
  end
end
for i=1:1:d+1
  for k=1:1:i-1
    c(i)=c(i)-c(k)*l(i,k);
  end
  c(i)=c(i)/l(i,i);
end
for i=d+1:-1:1
  for k=d+1:-1:i+1
    c(i)=c(i)-c(k)*l(k,i);
  end
  c(i)=c(i)/l(i,i);
end




% computed coefficients

c

% compute approximation to this exact polynomial with comp. coefficients c
approx = A*c;

figure(1)
plot(x,approx,'*- r', 'linewidth',1)
hold on
% plot of noisy  data
plot(x,y,'o b', 'linewidth',2)

hold off

str_xlabel = ['poly.degree d=', num2str(d)];

legend(str_xlabel,'noisy data');

xlabel('x')
  title('Polynomial fitting of random data via Cholesky  ')

% computation of the relative error  in coefficients to be recovered as
%  norm(approx. value - true value) / norm(true value)

% computation of the relative error  in function to be recovered as
%  norm(approx. value - true value) / norm(true value)


relerror_func = norm(y'- approx)/norm(y')
