m = 200;
[w1,w2] = meshgrid(linspace(-2,2,m));

w_true = (rand(2,1)-.5)*3;
X = orth(rand(2))*diag((rand(2,1)+a)/(1+a))*orth(rand(2))';
y = X*w_true;
a = .5;

f = reshape(sum((repmat(y,1,m^2) - X*[w1(:) w2(:)]').^2,1),m,m);

figure(1);clf;
contour(w1,w2,f,20);
set(gca,'fontsize',20)
xlabel('$w_1$','Interpreter','Latex')
ylabel('$w_2$','Interpreter','Latex')
axis image;
title('$f({\bf w})$','Interpreter','Latex')
hold on;
scatter(w_true(1),w_true(2),50,'r','o','filled')

%%


figure(3);clf;hold on
contour(w1,w2,f,20);
set(gca,'fontsize',20)
xlabel('$w_1$','Interpreter','Latex')
ylabel('$w_2$','Interpreter','Latex')
scatter(w_true(1),w_true(2),50,'r','o','filled')

axis image

wk = [2;2];
scatter(wk(1),wk(2),50,1,'d','filled');
tau = 1/norm(X)^2;
for k = 1:6
   grad = X'*(X*wk-y);
   newiter = wk - tau*grad;
   plot([wk(1),newiter(1)],[wk(2),newiter(2)],'k','linewidth',2)
   wk = newiter;   
   scatter(wk(1),wk(2),50,1,'d','filled');
   title(['Gradient descent, iteration ' num2str(k)],'Interpreter','Latex')
   pause
end