1
$\begingroup$

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}]

enter image description here

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?

$\endgroup$
  • $\begingroup$ 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}] $\endgroup$ – JimB Oct 29 at 2:51
  • $\begingroup$ @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. $\endgroup$ – Seyhmus Güngören Oct 29 at 15:26
  • $\begingroup$ @SeyhmusGüngören Last piece of code with While[] not so clear. What actually do you try to get? $\endgroup$ – Alex Trounev Oct 30 at 11:49
  • $\begingroup$ @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. $\endgroup$ – Seyhmus Güngören Oct 30 at 22:22
  • $\begingroup$ @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. $\endgroup$ – Alex Trounev Oct 31 at 8:32
2
$\begingroup$

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]] 

Figure 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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Seyhmus Güngören Oct 30 at 22:27
  • $\begingroup$ @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. $\endgroup$ – Alex Trounev Oct 30 at 23:26
  • $\begingroup$ 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? $\endgroup$ – Seyhmus Güngören Oct 31 at 0:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.