function [w] = graddescent_vectorized(X, Y, w_init)
  
  w = w_init;
  epsilon = 0.000001;
  for iteration = 1 : 5000 %<- in practice, use some stopping criterion!
    %COMPUTE GRADIENT
    grad = 2*sum(repmat(w'*X-Y,size(X,1),1).*X,2);

    %TAKE STEP
    w = w - epsilon * grad;

    %OUTPUT THE CURRENT ERROR
    err = sum((Y - w'*X).^2)

  end;