# FindMinimum gets confused by If[...]

I'm confused by why the following doesn't work. I'm trying to numerically minimize a function that has an If statement in its definition, as follows:

means = {3.0, 1.5, 2.1};
sigmashi = {0.5, 0.4, 0.3};
sigmaslo = {0.4, 0.6, 0.25};
coeffA = {1.2, 3.2, 2.0};
coeffB = {0.3, -0.4, 0.5};
chisq[x_] :=
Sum[If[coeffA[[i]] x + coeffB[[i]] >
means[[i]], (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
sigmashi[[i]]^2, (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
sigmaslo[[i]]^2], {i, 1, Length[means]}]
FindMinimum[chisq[x], {x, 1.0}]


This returns a long string of error messages:

Part::pspec: Part specification i is neither a machine-sized integer nor a list of machine-sized integers. >> Part::pspec: Part specification i is neither a machine-sized integer nor a list of machine-sized integers. >> Part::pspec: Part specification i is neither a machine-sized integer nor a list of machine-sized integers. >> General::stop: Further output of Part::pspec will be suppressed during this calculation. >> FindMinimum::nrnum: The function value ({0.3,-0.4,0.5}[[i]]+1. <<1>>-{3.,1.5,2.1}[[i]])^2/{0.4,0.6,0.25}[[i]]^2+(2 ({0.3,-0.4,0.5}[[i]]+<<1>>-{3.,1.5,2.1}[[i]])^2)/{0.5,0.4,0.3}[[i]]^2 is not a real number at {x} = {1.}. >>

Why doesn't FindMinimum like my function? Why does it even care about the "i"? Is there a way to force Mathematica to do the minimization purely numerically, and not care about how the function computes its result? (NMinimize returns the same errors.)

I can kluge my way around this with the following alternative:

myif[z_, x_, y_] := ((Tanh[100 z] + 1)/2) x + ((Tanh[-100 z] + 1)/2) y
chisqA[x_] :=
Sum[myif[coeffA[[i]] x + coeffB[[i]] -
means[[i]], (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
sigmashi[[i]]^2, (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
sigmaslo[[i]]^2], {i, 1, Length[means]}]


But it isn't pretty.

• try chisq[x_?NumericQ] !Mathematica graphics Commented Jan 7, 2014 at 1:31
• Thanks! Do you understand why it needs to be forced to use numerics in this case? Commented Jan 7, 2014 at 4:05
• Commented Jan 7, 2014 at 4:07

Using the initiative of bills and rewriting your function:

csq[x_] :=
Total[UnitStep[x #1 + #2 - #3] (#1 x + #2 - #3)^2/#4^2 +
UnitStep[-(x #1 + #2) + #3] (#1 x + #2 - #3)^2/#5^2 & @@@
Transpose[{coeffA, coeffB, means, sigmashi, sigmaslo}]]


then:

NMinimize[csq[x], x]


yields:

{21.6448, {x -> 0.798905}}

Or

FindMinimum[csq[x], {x, 1}]


yields:

{21.6448, {x -> 0.798905}}

Just for confirmation: plotting original function and min, redefined function and min, difference between two functions:

Code:

min = FindMinimum[csq[x], {x, 1}];
GraphicsRow[{Plot[chisq[x], {x, -3, 3}, PlotLabel -> "chisq(x)",
Epilog -> {Red, PointSize[0.04],
Point[{x /. min[[2]], min[[1]]}]}],
Plot[csq[x], {x, -3, 3},
Epilog -> {Red, PointSize[0.04], Point[{x /. min[[2]], min[[1]]}]},
PlotLabel -> "csq(x)"],
Plot[chisq[x] - csq[x], {x, -3, 3},
PlotLabel -> "chisq(x)-csq[x]"]}, ImageSize -> 700]


Let's see if we can figure out what's really causing the problem. Here's a simplified version of your problem:

f[x_] := If[x > 0, x^2, x^4];
Plot[f[x], {x, -1, 1}]
FindMinimum[f[x], {x, 0.1}]


{0., {x -> 0.}}


So there's no problem really If alone. Let's try with the indices into the lists:

a = {1, 2}; b = {3, 4};
g[x_] := Sum[If[a[[i]] x > 0, a[[i]] x^2, b[[i]] x^4], {i, 1, 2}];
FindMinimum[g[x], {x, 0.1}]
FindMinimum::nrnum: "The function value 0.02\ {1,2}[[i]] is not a real number at {x} = {0.1}"


So there is the problem. Here's a way to fix it, using the UnitStep function:

a = {1, 2}; b = {3, 4};
h[x_] := Sum[UnitStep[a[[i]] x]*a[[i]] x^2 + UnitStep[-a[[i]] x]*b[[i]] x^4, {i, 1, 2}];
FindMinimum[h[x], {x, 0.1}]
{2.23889*10^-29, {x -> -4.22896*10^-8}}


I'm sure you can figure out how to go back and apply UnitStep to your function.

• Thanks. Using UnitStep[] is basically an improved version of the kluge I was doing with Tanh[]. (For some reason I was worried that a non-smooth function would trip it up.) I think Nasser's ?NumericQ trick is somewhat nicer, although I still wonder why Mathematica doesn't get this right from the outset. Commented Jan 7, 2014 at 4:07
• @MattReece It's because Mma tries to use chisq'[x]` to find the minimum -- like every good calculus student ;) Commented Jan 7, 2014 at 4:19