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

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

-
can you add some example data? the code you gave doesn't really work without the data (it gives a different error and can be confusing). It's a bad idea to rely on externally-linked files as they'll eventually disappear, and then this question is impossible to understand. –  acl Oct 27 '12 at 23:05
@acl: Added data set. –  Problemaniac Oct 29 '12 at 1:16
You should probably reformulate your last follow up question as a new question. –  Andy Ross Oct 29 '12 at 2:31
this is a problem from Coursera, right? –  magma Nov 5 '12 at 1:51

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. - Could you explain a little that why when Sigmoid returns exactly 0 or 1, there will then be a Indeterminate? Is it that y_i - Sigmoid[...] is somehow used as a denominator in the computation? – Problemaniac Oct 29 '12 at 1:50 Because your cost function contains Log[h[...]] and Log[1-h[...]] If h[...] gives 0 or 1 one of the logs would return Indeterminate inside NMinimize since -Infinity is not a real number. – Andy Ross Oct 29 '12 at 2:26 I'm still not exactly sure why converting Indeterminate to $MaxMachineNumber is safe that no $-\infty$ will not be mistakenly converted. –  Problemaniac Oct 29 '12 at 3:14
Because over the machine reals that NMinmize will try there is no instance that gives -Infinity. –  Andy Ross Oct 29 '12 at 4:09