Error using symengine Singularity.

My code is showing Singularity problem. When int_theta is set to zero, there is no problem. But when it is set to values other than zero, I am getting the error for 'jc' stated below. The 'g0' gives me result perfectly but 'jc' shows error. Please can anyone show me where has it gone worng?

ERROR:

Error using symengine

Singularity.

Error in sym/double (line 692)

Xstr = mupadmex('symobj::double', S.s, 0);

Error in Logistic (line 51)

jc= double(subs(J,theta,int_theta));

syms a b c

data = load('exdata.txt'); % Load Data from a source
X = data(:, [1, 2]); 
y = data(:, 3);
[m, n] = size(X); % find X
X = [ones(m, 1) X]; % Add x0 as 1
int_theta = [-2;1;2];  % set initial theta values

theta = [a;b;c] ; % set theta as variable
g= 1./(1+ exp(-X*theta)); % set sigmoid function as eqn
g0= double(subs(g,theta,int_theta)); % find value of Sig funct for int_theta
J=  (-1/m).* sum(y.*log(g)+ (1-y).*log(1-g)); % cost funct
jc= double(subs(J,theta,int_theta)); % find value of cost funct for int_theta

j= gradient(J, theta); % find dJ/d(theta)
double(subs(j,theta,int_theta));

[Walter Roberson]

please attach the txt file

[RUHIN CHOWDHURY]

Here is the txt file.

Answer 1

[Walter Roberson]

Some of your terms are of the form

log(1 - 1/(exp(roughly -180)+1))

To the precision you are using, exp(roughly -180) + 1 is numerically the same as 1, so that becomes log(1 - 1/1) which is log(0) and that is generating the singularity, even though if you do the parts one-by-one you would get -inf instead of a singularity.

If you increase the digits used for calculation then the problem can be avoided.

syms a b c
data = load('exdata.txt'); % Load Data from a source
X = data(:, [1, 2]); 
y = data(:, 3);
[m, n] = size(X); % find X
X = [ones(m, 1) X]; % Add x0 as 1
int_theta = [-2;1;2]; % set initial theta values
theta = [a;b;c] ; % set theta as variable
g= 1./(1+ exp(-X*theta)); % set sigmoid function as eqn
olddigits = digits(100);
g0= double(subs(g,theta,int_theta)); % find value of Sig funct for int_theta
J= (-1/m).* sum(y.*log(g)+ (1-y).*log(1-g)); % cost funct
jc= double(subs(J,theta,int_theta)); % find value of cost funct for int_theta
j= gradient(J, theta); % find dJ/d(theta)
double(subs(j,theta,int_theta));
digits(olddigits);

[RUHIN CHOWDHURY]

Thank you so much.