clear all

x = 0.05:0.05:1;       % 20 evenly spaced points, need this to use vander
x_samples = x(4:4:20); % pick out the sample points
A = fliplr(vander(x)); % 20x20 vandermonde matrix
A = A(4:4:20,1:12);    % pick out 5x12 submatrix
y_samples = cos(4*x_samples)';

% solve regularized normal equations
mu = 0.001;
S = mu*diag(ones(12,1));
% S = diag([0.1,0.1,0.1,10,10,10,10,10,10,10,10,10]);
fprintf('cond(A^T*A) = %d\n',cond(A'*A));
fprintf('cond(A^T*A + S^T*S) = %d\n\n',cond(A'*A + S'*S));
b_regularized = (A'*A + S'*S) \ (A'*y_samples);

% print out ||r(b)||_2 and ||b||_2
r_regularized = y_samples - A*b_regularized;
fprintf('norm of regularized residual = %d\n',norm(r_regularized));
fprintf('norm of b_regularized = %d\n\n',norm(b_regularized));

% print out the coefficients obtained from the regularized problem
b_regularized

% plot the regularized best fit
xx = linspace(0,1);
figure(1)
% evaluate and plot the best fit at xx using polyval
p_regularized = polyval(flipud(b_regularized),xx);
plot(xx,p_regularized,'linewidth',2);
hold on
plot(x_samples,y_samples,'kx','markersize',14,'linewidth',2);
title(['regularized: ||r||_2 = ',num2str(norm(r_regularized)), ...
    ', ||b||_2 = ',num2str(norm(b_regularized))],'fontsize',16)

% % solve using backslash
% b_backslash = A \ y_samples
% 
% r_backslash = y_samples - A*b_backslash;
% fprintf('norm of backslash residual = %d\n',norm(r_backslash));
% fprintf('norm of b_backslash = %d\n\n',norm(b_backslash));
% 
% figure(2);
% p_backslash = polyval(flipud(b_backslash),xx);
% plot(xx,p_backslash,'linewidth',2);
% hold on
% plot(x_samples,y_samples,'kx','markersize',14,'linewidth',2);
% title(['backslash: ||r||_2 = ',num2str(norm(r_backslash)), ...
%     ', ||b||_2 = ',num2str(norm(b_backslash))],'fontsize',16)
