function [theta, theta_hist , J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESCENT(X, y, theta, alpha, num_iters) updates theta by 
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
    m = length(y); % number of training examples
    J_history = zeros(num_iters, 1);
    theta_hist = zeros(num_iters, length(theta));

    for iter = 1:num_iters
        theta = theta - alpha*1/m*X'*(X*theta-y);

        % Save the cost J and theta in every iteration    
        J_history(iter,1) = computeCost(X, y, theta);
        theta_hist(iter,:) = theta;

    end

end
