I'm trying my hand at regularized LR, simple with this formulas in matlab:

The cost function: 

J(theta) = 1/m*sum((-y_i)*log(h(x_i)-(1-y_i)*log(1-h(x_i))))+(lambda/2*m)*sum(theta_j)

The gradient: 

∂J(theta)/∂theta_0 = [(1/m)*(sum((h(x_i)-y_i)*x_j)] if j=0

∂j(theta)/∂theta_n = [(1/m)*(sum((h(x_i)-y_i)*x_j)]+(lambda/m)*(theta_j) if j>1

This is not matlab code is just the formula. 

So far I've done this: 

function [J, grad] = costFunctionReg(theta, X, y, lambda)

J = 0;
grad = zeros(size(theta));

temp_theta = [];

%cost function

%get the regularization term

for jj = 2:length(theta)

    temp_theta(jj) = theta(jj)^2;
end

theta_reg = lambda/(2*m)*sum(temp_theta);

temp_sum =[];

%for the sum in the cost function

for ii =1:m

   temp_sum(ii) = -y(ii)*log(sigmoid(theta'*X(ii,:)'))-(1-y(ii))*log(1-sigmoid(theta'*X(ii,:)'));

end

tempo = sum(temp_sum);

J = (1/m)*tempo+theta_reg;

%regulatization
%theta 0

reg_theta0 = 0;

for jj=1:m
 reg_theta0(jj) = (sigmoid(theta'*X(m,:)') -y(jj))*X(jj,1)
end    

reg_theta0 = (1/m)*sum(reg_theta0)

grad_temp(1) = reg_theta0

%for the rest of thetas

reg_theta  = [];
thetas_sum = 0;

for ii=2:size(theta)
    for kk =1:m
        reg_theta(kk) = (sigmoid(theta'*X(m,:)') - y(kk))*X(kk,ii)
    end
    thetas_sum(ii) = (1/m)*sum(reg_theta)+(lambda/m)*theta(ii)
    reg_theta = []
end

for i=1:size(theta)

    if i == 1
        grad(i) = grad_temp(i)
    else
        grad(i) = thetas_sum(i)
    end
end
end

And the cost function is giving correct results, but I have no idea why the gradient (one step) is not, the cost gives J = 0.6931 which is correct and the gradient grad = 0.3603 -0.1476 0.0320, which is not, the cost starts from 2 because the parameter theta(1) does not have to be regularized, any help? I guess there is something wrong with the code, but after 4 days I can't see it.Thanks 

Vectorized: 

function [J, grad] = costFunctionReg(theta, X, y, lambda)

hx = sigmoid(X * theta);
m = length(X);

J = (sum(-y' * log(hx) - (1 - y')*log(1 - hx)) / m) + lambda * sum(theta(2:end).^2) / (2*m);
grad =((hx - y)' * X / m)' + lambda .* theta .* [0; ones(length(theta)-1, 1)] ./ m ;

end