# Non-Linear optimization with NMaximize fails to find the parameters

I am trying to solve an optimization problem with several parameters and here is my code:

m = 4;

\[Kappa] = 10;

smin = 0;

smax = 80;

step = (smax - smin)/m;

thtmp = Table[Range[smin, smax, step], {k, 1, \[Kappa]}];

th = thtmp + 0.01;

i = 0;

threshold = 0.001;

\[Omega] = {10, 10, 10, 10, 10, 10, 10, 10, 10, 10};

\[Gamma] = 10.^(Range[-3, 2, 5/(\[Kappa] - 1)]/10); // Rationalize

f0[k_, y_] := PDF[ChiSquareDistribution[\[Omega][[k]]], y]

f00[k_, z_] := Piecewise[{{Gamma[\[Omega][[k]]/2, 0, z/2]/Gamma[\[Omega][[k]]/2], z > 0}}]

f1[k_, y_] := Piecewise[{{y^(-1 + \[Omega][[k]]/2)/((2*(\[Gamma][[k]] + 1))^(\[Omega][[k]]/2)*
E^(y/(2*(\[Gamma][[k]] + 1))))/Gamma[\[Omega][[k]]/2], y > 0}}]

f11[k_, z_] := Piecewise[{{Gamma[\[Omega][[k]]/2, 0, z/(2 (\[Gamma][[k]] + 1))]/
Gamma[\[Omega][[k]]/2], z > 0}}]

p0[k_, idx_] := N[f00[k, thtmp[[k]][[idx + 1]]] - f00[k, thtmp[[k]][[idx]]], 30]

p1[k_, idx_] := N[f11[k, thtmp[[k]][[idx + 1]]] - f11[k, thtmp[[k]][[idx]]], 30]

p00[k_, idx_, \[Lambda]_] := N[f00[k, \[Lambda]] - f00[k, thtmp[[k]][[idx]]], 30]

p01[k_, idx_, \[Lambda]_] := N[f00[k, thtmp[[k]][[idx + 2]]] - f00[k, \[Lambda]], 30]

p10[k_, idx_, \[Lambda]_] := N[f11[k, \[Lambda]] - f11[k, thtmp[[k]][[idx]]], 30]

p11[k_, idx_, \[Lambda]_] := N[f11[k, thtmp[[k]][[idx + 2]]] - f11[k, \[Lambda]], 30]

utmp1[k_, idx_, \[Lambda]_] := Log[p10[k, idx, \[Lambda]]/p00[k, idx, \[Lambda]]];

utmp2[k_, idx_, \[Lambda]_] := Log[p11[k, idx, \[Lambda]]/p01[k, idx, \[Lambda]]];

utmp[k_, idx_] := Log[p1[k, idx]/p0[k, idx]]

utable[k_, idx_] := Join[Table[utmp[k, n], {n, 1, idx - 1}], Table[utmp[k, n], {n, idx + 2, m}]];

p0table[k_, idx_] := Join[Table[p0[k, n], {n, 1, idx - 1}], Table[p0[k, n], {n, idx + 2, m}]];

p1table[k_, idx_] := Join[Table[p1[k, n], {n, 1, idx - 1}], Table[p1[k, n], {n, idx + 2, m}]];

c02[k_, idx_] := (utable[k, idx]^2).p0table[k, idx];

c12[k_, idx_] := (utable[k, idx]^2).p1table[k, idx];

c01[k_, idx_] := utable[k, idx].p0table[k, idx];

c11[k_, idx_] := utable[k, idx].p1table[k, idx];

c0[k_, idx_] := Total[p0table[k, idx]];

c1[k_, idx_] := Total[p1table[k, idx]];

Subscript[gm, 0][k_] := Sum[p0[k, idx]*utmp[k, idx], {idx, 1, m}];

Subscript[gm, 1][k_] := Sum[p1[k, idx]*utmp[k, idx], {idx, 1, m}];

Subscript[mm, 0][k_, idx_, \[Lambda]_] := c01[k, idx] + (p00[k, idx, \[Lambda]]*utmp1[k, idx, \[Lambda]] + p01[k, idx, \[Lambda]]*utmp2[k, idx, \[Lambda]]);

Subscript[mm, 1][k_, idx_, \[Lambda]_] := c11[k, idx] + (p10[k, idx, \[Lambda]]*utmp1[k, idx, \[Lambda]] + p11[k, idx, \[Lambda]]*utmp2[k, idx, \[Lambda]]);

Subscript[c\[Mu], 0][k_, idx_] := Sum[Subscript[gm, 0][n], {n, 1, \[Kappa]}] - (p0[k, idx]*utmp[k, idx] + p0[k, idx + 1]*utmp[k, idx + 1]);

Subscript[c\[Mu], 1][k_, idx_] := Sum[Subscript[gm, 1][n], {n, 1, \[Kappa]}] - (p1[k, idx]*utmp[k, idx] + p1[k, idx + 1]*utmp[k, idx + 1]);

Subscript[s, 0][k_] := Sum[Sum[p0[j, idx]*(utmp[j, idx] - Subscript[gm, 0][j])^2, {idx, 1,
m}], {j, 1, \[Kappa]}] - Sum[p0[k, idx]*(utmp[k, idx] - Subscript[gm, 0][k])^2, {idx, 1, m}];

Subscript[s, 1][k_] := Sum[Sum[p1[j, idx]*(utmp[j, idx] - Subscript[gm, 1][j])^2, {idx, 1,
m}], {j, 1, \[Kappa]}] - Sum[p1[k, idx]*(utmp[k, idx] - Subscript[gm, 1][k])^2, {idx, 1, m}];

Subscript[c\[Sigma], 0][k_, idx_, \[Lambda]_] := c02[k, idx] - 2*c01[k, idx]*Subscript[mm, 0][k, idx, \[Lambda]] + Subscript[mm, 0][k, idx, \[Lambda]]^2*c0[k, idx];

Subscript[c\[Sigma], 1][k_, idx_, \[Lambda]_] := c12[k, idx] - 2*c11[k, idx]*Subscript[mm, 1][k, idx, \[Lambda]] + Subscript[mm, 1][k, idx, \[Lambda]]^2*c1[k, idx];

Subscript[\[Mu], 0][k_, idx_, \[Lambda]_] := utmp1[k, idx, \[Lambda]]*p00[k, idx, \[Lambda]] +utmp2[k, idx, \[Lambda]]*p01[k, idx, \[Lambda]] + Subscript[c\[Mu], 0][k, idx]

Subscript[\[Mu], 1][k_, idx_, \[Lambda]_] := utmp1[k, idx, \[Lambda]]*p10[k, idx, \[Lambda]] +utmp2[k, idx, \[Lambda]]*p11[k, idx, \[Lambda]] + Subscript[c\[Mu], 1][k, idx]

Subscript[\[Sigma], 0][k_, idx_, \[Lambda]_] := Sqrt[(p00[k, idx, \[Lambda]]*(utmp1[k, idx, \[Lambda]] - Subscript[mm, 0][k, idx, \[Lambda]])^2 + p01[k, idx, \[Lambda]]*(utmp2[k, idx, \[Lambda]] - Subscript[mm, 0][k, idx, \[Lambda]])^2) + Subscript[c\[Sigma], 0][k, idx, \[Lambda]] + Subscript[s, 0][k]]

Subscript[\[Sigma], 1][k_, idx_, \[Lambda]_] := Sqrt[(p10[k, idx, \[Lambda]]*(utmp1[k, idx, \[Lambda]] - Subscript[mm, 1][k, idx, \[Lambda]])^2 + p11[k, idx, \[Lambda]]*(utmp2[k, idx, \[Lambda]] - Subscript[mm, 1][k, idx, \[Lambda]])^2) + Subscript[c\[Sigma], 1][k, idx, \[Lambda]] + Subscript[s, 1][k]]

a[k_, idx_, \[Lambda]_] := Total[Join[
Table[Min[Table[utmp[jj, n], {n, 1, m}]], {jj, 1, k - 1}],
Table[Min[Table[utmp[jj, n], {n, 1, m}]], {jj, k + 1, \[Kappa]}]]] + Min[Join[Table[utmp[k, n], {n, 1, idx - 1}], Table[utmp[k, n], {n, idx + 2, m}]], utmp1[k, idx, \[Lambda]],
utmp2[k, idx, \[Lambda]]];

b[k_, idx_, \[Lambda]_] := Total[Join[
Table[Max[Table[utmp[jj, n], {n, 1, m}]], {jj, 1, k - 1}],
Table[Max[Table[utmp[jj, n], {n, 1, m}]], {jj,
k + 1, \[Kappa]}]]] + Max[Join[Table[utmp[k, n], {n, 1, idx - 1}],
Table[utmp[k, n], {n, idx + 2, m}]], utmp1[k, idx, \[Lambda]],
utmp2[k, idx, \[Lambda]]];

result[\[Kappa]x_, a0_, b0_, a1_, b1_] := FindRoot[PDF[SkewNormalDistribution[a0, b0, \[Kappa]x], t] - PDF[SkewNormalDistribution[a1, b1, \[Kappa]x], t], {t, 0}]

rr[k_, idx_, \[Lambda]_, \[Kappa]x_] := -Erf[(Subscript[\[Mu], 0][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]][[1, 2]]))/(Sqrt[2] Subscript[\[Sigma], 0][k, idx, \[Lambda]])] + Erf[(Subscript[\[Mu], 1][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]][[1, 2]]))/(Sqrt[2] Subscript[\[Sigma], 1][k, idx, \[Lambda]])] - 4 OwenT[(Subscript[\[Mu], 0][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]][[1, 2]]))/
Subscript[\[Sigma], 0][k, idx, \[Lambda]], \[Kappa]x] + 4 OwenT[(Subscript[\[Mu], 1][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]][[1, 2]]))/
Subscript[\[Sigma], 1][k, idx, \[Lambda]], \[Kappa]x]


I am able to evaluate as well as plot the objective function 'rr' as follows:

 rr[1, 1, 10, 1]

1.98504

Plot[rr[1, 1, \[Lambda], 1], {\[Lambda], 0, 40}]


But I am unable to get any result for the following optimization problem here:

 NMaximize[{rr[1, 1, \[Lambda], \[Kappa]x], {thtmp[[1]][[1]] <= \[Lambda] <=
thtmp[[1]][[1 + 2]], -10 <= \[Kappa]x <= 10}}, {\[Lambda], \[Kappa]x}, WorkingPrecision -> 30]


My ultimate aim is actually to solve this bigger problem but without the previous one I cannot succeed here:

 While[Total[Total[Abs[thtmp - th]]] > threshold, {i = i + 1; Print[Total[Total[Abs[thtmp - th]]], thtmp];, th = thtmp;, results =
Table[NMaximize[{rr[k, idx, \[Lambda], \[Kappa]x],
thtmp[[k]][[idx]] <= \[Lambda] <=
thtmp[[k]][[idx + 2]], -10 <= \[Kappa]x <=
10}, {\[Lambda], \[Kappa]x}, WorkingPrecision -> 30], {k,
1, \[Kappa]}, {idx, 1, m - 1}];, Table[thtmp[[k]][[2 ;; m]] = \[Lambda] /.
results[[k]][[All, 2]], {k, 1, \[Kappa]}]}]


I already set-up a certain precision to the problem but it still tells me that I have problems with it. Moreover, it tells me that I got results which are not numbers, when the objective funtion is evaluated. This is also not clear to me.

How should one deal with this problem?

• What do you mean by "I am unable to get any result" ? Are there error messages? Is the command just returned? Does it just keep executing until you stop it? The "bigger problem" you describe does mention that the results are not numbers. Do you get the same error/warning for the initial problem? Maybe a contour plot would help get better starting values: ContourPlot[rr[1, 1, \[Lambda], \[Kappa]x], {\[Lambda], 0, 40}, {\[Kappa]x, -10, 10}]
– JimB
Oct 29, 2020 at 2:51
• @JimB I am getting some bullshit as a result from Mathematica. Please run the code and you will see. In the bigger problem there are also warnings. In the smaller one, simply no outcome. I don’t interrupt the program. It finds something but it’s just a list of nonsense expressions. Oct 29, 2020 at 15:26
• @SeyhmusGüngören Last piece of code with While[] not so clear. What actually do you try to get? Oct 30, 2020 at 11:49
• @AlexTrounev I am iterating the results until convergence. Take [Kappa]x=2 and run the while loop as given in the question. You will see that it will find the best thresholds after several iterations. [Kappa]x=2 mean Gaussian density. I am trying the get even better thresholds considering generalized Guassians with parameter [Kappa]x. So there is the second optimization over [Kappa]x. Considering the below written answer, I was able to run the code. But it doesnt work. I think it is because of FindRoot. There are more than one root and it is possible that FindRoot is finding the wrong one. Oct 30, 2020 at 22:22
• @SeyhmusGüngören It means there is no unique solution of this problem. In this case we can add some constraints to get right solution. Oct 31, 2020 at 8:32

We can improve code in two step. First, we use Compile[] for result and, second, Module for rr:

result = Compile[{{\[Kappa]x, _Real}, {a0, _Real}, {b0, _Real}, {a1, \
_Real}, {b1, _Real}},
t /. FindRoot[
PDF[SkewNormalDistribution[a0, b0, \[Kappa]x], t] -
PDF[SkewNormalDistribution[a1, b1, \[Kappa]x], t], {t, 1/10}]];
rr[k_, idx_, \[Lambda]_, \[Kappa]x_] := Module[{},

rrr = -Erf[(Subscript[\[Mu], 0][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]]))/(Sqrt[
2] Subscript[\[Sigma], 0][k, idx, \[Lambda]])] +
Erf[(Subscript[\[Mu], 1][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]]))/(Sqrt[
2] Subscript[\[Sigma], 1][k, idx, \[Lambda]])] -
4 OwenT[(Subscript[\[Mu], 0][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]]))/
Subscript[\[Sigma], 0][k, idx, \[Lambda]], \[Kappa]x] +
4 OwenT[(Subscript[\[Mu], 1][k,
idx, \[Lambda]] - (result[\[Kappa]x,
Subscript[\[Mu], 0][k, idx, \[Lambda]],
Subscript[\[Sigma], 0][k, idx, \[Lambda]],
Subscript[\[Mu], 1][k, idx, \[Lambda]],
Subscript[\[Sigma], 1][k, idx, \[Lambda]]]))/
Subscript[\[Sigma], 1][k, idx, \[Lambda]], \[Kappa]x]; rrr];


Now we can plot rr[1,1,x,y] using code

lst = Table[{x, y, rr[1,1,x, y]}, {x, .1, 40, 1}, {y, -10, 10, .5}];

ListPlot3D[Flatten[lst, 1]]


Finally we compute

NMaximize[{rr[1, 1, x, y],
Element[{x, y}, Rectangle[{0, -10}, {40, 10}]]}, {x, y}] // Quiet

Out[]= {2., {x -> 20.5113, y -> 9.61341}}


We can recalculate this result with

FindMaximum[{rr[1, 1, x, y],
Element[{x, y}, Rectangle[{0, -10}, {40, 10}]]}, {{x, 20.5}, {y,
9.6}}] // Quiet

Out[]= {2., {x -> 20.39, y -> 6.70582}}


Therefore we get same max value 2 but in different points. As it shown in Figure 1 function rr looks like constant=2 in these points.

• Thank you very much for the answer. It is working well, although still with many warnings. I have one more issue I think. The bigger problem with 'While'. The results are not converging unfortunately. For example if I choose [Kappa]x=2 or [Kappa]x=4 and run the while loop over one dimensional optimization. Then the algorithm is converging. I suspect that FindRoot is finding the wrong root for other values of [Kappa]x if the optimization runs also on [Kappa]x. Is the result of FindRoot currently always real number? Oct 30, 2020 at 22:27
• @SeyhmusGüngören The problem with FindRoot[], Compile and NMaximize[] is not so clear. It is not affected result, but looks like a bug. I think that it depends on some function definition including Real and Complex. Oct 30, 2020 at 23:26
• rr[a0_, b0_, a1_, b1_, [Kappa]x_, t_] := PDF[SkewNormalDistribution[a0, b0, [Kappa]x], t] - PDF[SkewNormalDistribution[a1, b1, [Kappa]x], t] and assume FindRoot[rr[-5, 1, -1, 1, .2, t], {t, 1/10}] then we find the wrong answer and if we use FindRoot[rr[-5, 1, -1, 1, .2, t], {t, -2}] then we get the correct answer... so sometimes we get wrong answers.. is it possible to prevent wrong answers? Oct 31, 2020 at 0:20