n = 30;
rang = 5;
step = 1;
x = 5*(rand(n,1)-0.5);
y = (1+randn(1))*x+(randn(n,1));

A = [x ones(n,1)];

yy = y'*y;
Ay = A'*y
AA = A'*A

sol = AA\Ay
% set display range
a1 = (-rang:step:rang)+(sol(1));
a2 = (-rang:step:rang)+(sol(2));
k = size(a1,2);

F  = zeros(k);
dF1 = zeros(k);
dF2 = zeros(k);

for i = 1:k
  for j = 1:k
    a = [a1(i);a2(j)];
    % F(x) = x'A'Ax - 2x'A'y + y'y
    F(j,i)= a'*AA*a - 2*a'*Ay + yy;
    % dF(x) = 2A'Ax -2A'y
    gv = 2*AA*a - 2*Ay;
    dF1(i,j) = gv(1);
    dF2(i,j) = gv(2);
  end
end

% normalize
maxdf = max(max(abs(dF1)+abs(dF2)));
dF1 = dF1/maxdf;
dF2 = dF2/maxdf;

%contourc(a1, a2, F);
contour(a1, a2, F, 30);
%mesh(a1, a2, F);
hold on;
for i = 1:k
  for j = 1:k
    plot([a1(i) a1(i)+dF1(i,j)],[a2(j) a2(j)+dF2(i,j)],'k',...
        a1(i)+dF1(i,j), a2(j)+dF2(i,j),'rx');
  end
end
plot(sol(1), sol(2), 'ok', 'markersize', 30);
hold off;
