How to find local minima of a multivariate function with singularity?

Question

I am working on a logistic regression problem which requires minimizing a cost function J[{theta0, theta1, theta2}, X, y] to find the optimal value for fitting parameter {theta0, theta1, theta2}.

X is a known M x 3 matrix of real numbers, and y an M-dimensional vector. M is the number of sample data points (approx. 300) in the training data set. The definition for J is

J[θ_, X_, y_] := -(1/Length[y]) (y.Log[h[θ, X]] + (1 - y).Log[1 - h[θ, X]])

where

h[θ_, X_] := Sigmoid[X.θ];
Sigmoid[z_] := 1/(1 + Exp[-z]);

When I try to minimize J with NMinimize and FindMinimum, both fail due to singularity:

NMinimize[J[{θ0, θ1, θ2}, X1, y1], {θ0, θ1, θ2}]
(* NMinimize::nnum: The function value Indeterminate is not a number at {θ0, θ1, θ2} = {0.673558,0.659492,0.0861047}. >> *)

and

FindMinimum[J[{θ0, θ1, θ2}, X1, y1], {θ0, 0}, {θ1, 0}, {θ2, 0}]
(* FindMinimum::nrnum: The function value Indeterminate is not a real number at {θ0, θ1, θ2} = {0.00303682,0.364698,0.342032}. >> *)

For comparison, in MATLAB, using fminunc, it suffers from NaN too

[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), [0;0;0], optimset('MaxIter', 200));
%Warning: Gradient must be provided for trust-region algorithm; using line-search algorithm instead. 
% In fminunc at 365 
%Error using roots (line 28)
%Input to ROOTS must not contain NaN or Inf.
%<snip>

however, setting the 'GradObj' option makes it work:

[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), [0;0;0], optimset('GradObj', 'on', 'MaxIter', 400));    
%Local minimum possible.
%fminunc stopped because the final change in function value relative to 
%its initial value is less than the default value of the function tolerance.    
%<stopping criteria details>

theta    
%theta =    
%  -24.9330
%    0.2044
%    0.1996

cost
% cost = 0.2035

Is there a way to tweak Mathematica to solve this? The notebook showing the details is here.

Edit (follow-up questions)

1. Is Indeterminate always an unachievable $+\infty$?

2. Is there a way to do automatic regularization (i.e. adding a regularization term $\lambda || \vec{\beta} ||^2 $ and optimizing $\lambda$) when fitting parameter vector $\vec{\beta}$ for a binary non-linear logistic regression by using some option/setting with LogitModelFit?

3. The best fit found with NMinimize is J_min as 0.203498 and {\[Beta]0 -> -25.1613, \[Beta]1 -> 0.206232, \[Beta]2 -> 0.201472} as compared to J_min as 0.203506 and [-24.932998 0.204408 0.199618] found with Matlab's fminunc. It seems Mathematica wins. Right?

4. Using LogitModelFit on the binary nonlinear logistic regression example, it found a decision boundary with a hole. Is there any way to control the complexity/topology of the decision boundary (thinking under-fitting and/or over-fitting)?

Appendix: complete code

The complete code and data is stored in Github. I intend to keep it for as long as Github allows.

Appendix: example data set

(as per comment)

34.62365962451697,78.0246928153624,0
30.28671076822607,43.89499752400101,0
35.84740876993872,72.90219802708364,0
60.18259938620976,86.30855209546826,1
79.0327360507101,75.3443764369103,1
45.08327747668339,56.3163717815305,0
61.10666453684766,96.51142588489624,1
75.02474556738889,46.55401354116538,1
76.09878670226257,87.42056971926803,1
84.43281996120035,43.53339331072109,1
95.86155507093572,38.22527805795094,0
75.01365838958247,30.60326323428011,0
82.30705337399482,76.48196330235604,1
69.36458875970939,97.71869196188608,1
39.53833914367223,76.03681085115882,0
53.9710521485623,89.20735013750205,1
69.07014406283025,52.74046973016765,1
67.94685547711617,46.67857410673128,0
70.66150955499435,92.92713789364831,1
76.97878372747498,47.57596364975532,1
67.37202754570876,42.83843832029179,0
89.67677575072079,65.79936592745237,1
50.534788289883,48.85581152764205,0
34.21206097786789,44.20952859866288,0
77.9240914545704,68.9723599933059,1
62.27101367004632,69.95445795447587,1
80.1901807509566,44.82162893218353,1
93.114388797442,38.80067033713209,0
61.83020602312595,50.25610789244621,0
38.78580379679423,64.99568095539578,0
61.379289447425,72.80788731317097,1
85.40451939411645,57.05198397627122,1
52.10797973193984,63.12762376881715,0
52.04540476831827,69.43286012045222,1
40.23689373545111,71.16774802184875,0
54.63510555424817,52.21388588061123,0
33.91550010906887,98.86943574220611,0
64.17698887494485,80.90806058670817,1
74.78925295941542,41.57341522824434,0
34.1836400264419,75.2377203360134,0
83.90239366249155,56.30804621605327,1
51.54772026906181,46.85629026349976,0
94.44336776917852,65.56892160559052,1
82.36875375713919,40.61825515970618,0
51.04775177128865,45.82270145776001,0
62.22267576120188,52.06099194836679,0
77.19303492601364,70.45820000180959,1
97.77159928000232,86.7278223300282,1
62.07306379667647,96.76882412413983,1
91.56497449807442,88.69629254546599,1
79.94481794066932,74.16311935043758,1
99.2725269292572,60.99903099844988,1
90.54671411399852,43.39060180650027,1
34.52451385320009,60.39634245837173,0
50.2864961189907,49.80453881323059,0
49.58667721632031,59.80895099453265,0
97.64563396007767,68.86157272420604,1
32.57720016809309,95.59854761387875,0
74.24869136721598,69.82457122657193,1
71.79646205863379,78.45356224515052,1
75.3956114656803,85.75993667331619,1
35.28611281526193,47.02051394723416,0
56.25381749711624,39.26147251058019,0
30.05882244669796,49.59297386723685,0
44.66826172480893,66.45008614558913,0
66.56089447242954,41.09209807936973,0
40.45755098375164,97.53518548909936,1
49.07256321908844,51.88321182073966,0
80.27957401466998,92.11606081344084,1
66.74671856944039,60.99139402740988,1
32.72283304060323,43.30717306430063,0
64.0393204150601,78.03168802018232,1
72.34649422579923,96.22759296761404,1
60.45788573918959,73.09499809758037,1
58.84095621726802,75.85844831279042,1
99.82785779692128,72.36925193383885,1
47.26426910848174,88.47586499559782,1
50.45815980285988,75.80985952982456,1
60.45555629271532,42.50840943572217,0
82.22666157785568,42.71987853716458,0
88.9138964166533,69.80378889835472,1
94.83450672430196,45.69430680250754,1
67.31925746917527,66.58935317747915,1
57.23870631569862,59.51428198012956,1
80.36675600171273,90.96014789746954,1
68.46852178591112,85.59430710452014,1
42.0754545384731,78.84478600148043,0
75.47770200533905,90.42453899753964,1
78.63542434898018,96.64742716885644,1
52.34800398794107,60.76950525602592,0
94.09433112516793,77.15910509073893,1
90.44855097096364,87.50879176484702,1
55.48216114069585,35.57070347228866,0
74.49269241843041,84.84513684930135,1
89.84580670720979,45.35828361091658,1
83.48916274498238,48.38028579728175,1
42.2617008099817,87.10385094025457,1
99.31500880510394,68.77540947206617,1
55.34001756003703,64.9319380069486,1
74.77589300092767,89.52981289513276,1

Answer 1

You could try something like the following.

bounded[Indeterminate] := $MaxMachineNumber;
bounded[x_?NumericQ] := x

This gives a very large number instead of Indeterminate so NMinimize keeps searching.

NMinimize[bounded[J[{θ0, θ1, θ2}, X1, y1]], {θ0, θ1, θ2}]    
(* {0.203498, {θ0 -> -25.1613, θ1 -> 0.206231, θ2 -> 0.201471}} *)

Incidentally, is there any reason you don't use LogitModelFit?

LogitModelFit[Transpose[Transpose[X1]~Join~{y1}],
    {θ0, θ1, θ2}, {θ0, θ1, θ2}, IncludeConstantBasis -> False]["BestFitParameters"]

(*{-25.1613, 0.206232, 0.201472}*)

Edit:

In response to the edits. The only way for Indeterminate to arise beyond the length of y being zero is for the sigmoid function to return exact zero or 1. This should never happen with real input.

Reduce[Sigmoid[z] == 0 || Sigmoid[z] == 1, z, Reals]

(*False *)

LogitModelFit cannot perform regularization as of version 8. You will have to do this using your NMinimize approach.