1
$\begingroup$

I need to minimize numerically an expression with constraints. The expression also includes a numerical integral (called 'PNoninteract'). The minimization apparently stops prematurely (I know roughly what the answer should be and I see it is incorrect), while producing error messages regarding having reached complex numbers (where it should be restricted to Reals by the constraints). I receive messages such as:


I'm currently using `Method->"RandomSearch"` because my previous attempts (with a simpler version of this problem) only worked with this method and produced similar error messages with all other methods (including the other stochastic methods).
Any idea how to resolve this? Thanks!

My code:

    P = {29/61, 18/61, 10/61, 4/61};
    K = 10; 
    M = 1; 
    Eth = -6; 
    Lmod = 40; 
    RealL = 18; 
    V = Array[v, 4];
    U = Array[u, {4, 4}];
    Q = Array[q, 6]; 
    Q[[1]] = V[[1]]*U[[1, 1]]; 
    Q[[2]] = V[[1]]*(U[[1, 2]] + U[[1, 3]] + U[[1, 4]]) + 
       V[[2]]*U[[2, 1]] + V[[3]]*U[[3, 1]] + V[[4]]*U[[4, 1]]; 
    Q[[3]] = V[[2]]*U[[2, 2]];
    Q[[4]] = V[[2]]*(U[[2, 3]] + U[[2, 4]]) + V[[3]]*U[[3, 2]] + 
       V[[4]]*U[[4, 2]];
    Q[[5]] = V[[3]]*U[[3, 3]] + V[[4]]*U[[4, 4]];
    Q[[6]] = V[[3]]*U[[3, 4]] + V[[4]]*U[[4, 3]];
   

     EnergyVals = {-1, 0, -0.75, -0.25, 1, -1.25};
        PNoninteract = 
          Probability[{x, y, z, w, c, d} . EnergyVals >= 
            Eth, {x, y, z, w, c, d} \[Distributed] 
            MultinomialDistribution[RealL, Q]];
        allM = \!\(
        \*UnderoverscriptBox[\(∑\), \(a = 
          1\), \(4\)]\(V[\([\)\(a\)\(]\)]*\(
        \*UnderoverscriptBox[\(∑\), \(b = 1\), \(4\)]U[\([\)\(a, 
            b\)\(]\)]*\((K*
              Log[U[\([\)\(a, b\)\(]\)]/P[\([\)\(b\)\(]\)]] + 
             Log[V[\([\)\(a\)\(]\)]/P[\([\)\(a\)\(]\)]])\)\)\)\) - 
       Log[1 - PNoninteract^M];
        Energy = {{-1, 0, 0, 0}, {0, -0.75, -0.25, -0.25}, {0, -0.25, 
           1, -1.25}, {0, -0.25, -1.25, 
           1}}; 
        Emean = \!\(
        \*UnderoverscriptBox[\(∑\), \(a = 
            1\), \(4\)]\(V[\([\)\(a\)\(]\)] \(
        \*UnderoverscriptBox[\(∑\), \(b = 1\), \(4\)]U[\([\)\(a, 
              b\)\(]\)]\ Energy[\([\)\(a, 
              b\)\(]\)]\)\)\); 
        vars = Flatten[{U, V}];
    constraint1 = 
      Table[{0 <= vars[[i]] <= 1, vars[[i]] ∈ Reals}, {i, 1, 
        Length[vars]}];
    constraint2 = Total[V] == 1;
    constraint3 = Total[U[[1]]] == 1;
    constraint4 = Total[U[[2]]] == 1;
    constraint5 = Total[U[[3]]] == 1;
    constraint6 = Total[U[[4]]] == 1;
    constraint7 = Emean >= Eth/Lmod;
    cons = {constraint1, constraint2, constraint3, constraint4, 
        constraint5, constraint6, constraint7} // Flatten;
    SolK10M1 = NMinimize[{allM, cons}, vars, Method -> "RandomSearch"]
Improve this question
$\endgroup$
2
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Commented Feb 22, 2022 at 12:57
  • $\begingroup$ Is PNoninteract = Probability[{x, y, z, w, c, d} . EnergyVals >= Eth, {x, y, z, w, c, d} \[Distributed] MultinomialDistribution[RealL, Q]] OK? $\endgroup$
    – user64494
    Commented Feb 22, 2022 at 17:29

2 Answers 2

Reset to default
1
$\begingroup$

Your code with the command

SolK10M1 =  NMinimize[{allM, cons}, vars, Method -> "DifferentialEvolution"]

in the last line results in

{6.74324,{u[1,1]->0.13901,u[1,2]->0.31005,u[1,3]->0.405679,u[1,4]->0.145261,u[2,1]->0.284419,u[2,2]->0.416764,u[2,3]->0.128494,u[2,4]->0.170323,u[3,1]->0.136242,u[3,2]->0.675924,u[3,3]->0.121422,u[3,4]->0.066412,u[4,1]->0.602008,u[4,2]->0.124211,u[4,3]->0.120767,u[4,4]->0.153013,v[1]->0.281198,v[2]->0.0885391,v[3]->0.521515,v[4]->0.108748}}

and several warnings about complex numbers when calculating

NMinimize::nnum: The function value 8.327 -Log[1-PNoninteract] is not a number at {u[1,1],u[1,2],u[1,3],u[1,4],u[2,1],u[2,2],u[2,3],u[2,4],u[3,1],u[3,2],<<10>>} = {0.588221,0.101819,0.106911,0.20305,0.172127,0.241654,0.542321,0.0438981,0.145015,0.591377,<<10>>}.
NMinimize::nrnum: The function value 24.9444 -8.51506 I is not a real number at {u[1,1],u[1,2],u[1,3],u[1,4],u[2,1],u[2,2],u[2,3],u[2,4],u[3,1],u[3,2],<<10>>} = {0.195743,0.494297,0.106911,0.20305,0.357974,0.241654,0.356474,0.0438981,0.455513,0.280879,<<10>>}.

Warnings are not errors. Your constraints are rigid (==) so complex numbers appear.

Addition. Here is a partial verification.

constraint1[[1]] /. SolK10M1[[2]]

{True,True}

constraint2[[1]] /. SolK10M1[[2]]

1.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Indeed the solution obeys the constraints, but the answer is not the minimum. I can show that with different precisions I obtain different answers and can get lower values of the expression. So, rephrasing my question, I am not so much concerned about the error message, but that the answer is wrong and that optimization stopped too early perhaps. $\endgroup$
    – PalmTree
    Commented Feb 23, 2022 at 8:05
  • $\begingroup$ @PalmTree: Can you elaborate your "but the answer is not the minimum. I can show that with different precisions I obtain different answers and can get lower values of the expression", giving us those results? In the other case this is empty words $\endgroup$
    – user64494
    Commented Feb 23, 2022 at 11:39
  • $\begingroup$ The simplest demonstration is that changing the precision gives different solutions. with machine precision: minimum: 5.3911, v[1]->0.230147,v[2]->0.0337145,v[3]->0.258751,v[4]->0.477387, but with precision=30 the answer is v[1]->0.31653328919904558702484368865,v[2]->0.241291197405444009015328296861,v[3]->0.356553425478305113877827071994,v[4]->0.0856220879172052900820009424981} $\endgroup$
    – PalmTree
    Commented Feb 23, 2022 at 19:57
  • $\begingroup$ Could you please explain why the algorithm reaches complex-valued numbers here when it is constrained to real values? It is not obvious to me. $\endgroup$
    – PalmTree
    Commented Feb 23, 2022 at 20:15
  • $\begingroup$ Just to complete the previous example, with Method->"RandomSearch", WorkingPrecision->30 the minimum I obtain is 4.075. So in total we have 3 different solutions for the same problem (including the one you obtained which was 6.7) depending on the choice of algorithm and numerical precision. Does that convince that the answer the algorithm yields is not the minimum here? $\endgroup$
    – PalmTree
    Commented Feb 23, 2022 at 20:26
1
$\begingroup$

You can increase the number of search points for the random search method as below:

ParallelDo[
 Print[
  NMinimize[{allM, cons}, vars, 
   Method -> {"RandomSearch", "SearchPoints" -> i}]
  ],
 {i, 500, 3000, 500}
 ]

The solution to the above is :

{2.83582,{u[1,1]->0.474797,u[1,2]->0.295426,u[1,3]->0.164126,u[1,4]->0.0656503,u[2,1]->0.492646,u[2,2]->0.285607,u[2,3]->0.158391,u[2,4]->0.0633564,u[3,1]->0.416165,u[3,2]->0.240842,u[3,3]->0.278646,u[3,4]->0.0643475,u[4,1]->0.44781,u[4,2]->0.259156,u[4,3]->0.173101,u[4,4]->0.119934,v[1]->0.331787,v[2]->0.145927,v[3]->0.438083,v[4]->0.0842041}}

{2.83582,{u[1,1]->0.474798,u[1,2]->0.295426,u[1,3]->0.164126,u[1,4]->0.0656503,u[2,1]->0.492646,u[2,2]->0.285607,u[2,3]->0.158391,u[2,4]->0.0633563,u[3,1]->0.416165,u[3,2]->0.240842,u[3,3]->0.278646,u[3,4]->0.0643474,u[4,1]->0.44781,u[4,2]->0.259156,u[4,3]->0.173101,u[4,4]->0.119934,v[1]->0.331787,v[2]->0.145926,v[3]->0.438083,v[4]->0.0842042}}

{2.83582,{u[1,1]->0.474798,u[1,2]->0.295426,u[1,3]->0.164126,u[1,4]->0.0656503,u[2,1]->0.492646,u[2,2]->0.285607,u[2,3]->0.158391,u[2,4]->0.0633563,u[3,1]->0.416165,u[3,2]->0.240842,u[3,3]->0.278646,u[3,4]->0.0643474,u[4,1]->0.44781,u[4,2]->0.259156,u[4,3]->0.173101,u[4,4]->0.119934,v[1]->0.331787,v[2]->0.145926,v[3]->0.438083,v[4]->0.0842042}}

{2.83582,{u[1,1]->0.474798,u[1,2]->0.295426,u[1,3]->0.164126,u[1,4]->0.0656503,u[2,1]->0.492646,u[2,2]->0.285607,u[2,3]->0.158391,u[2,4]->0.0633563,u[3,1]->0.416165,u[3,2]->0.240842,u[3,3]->0.278646,u[3,4]->0.0643474,u[4,1]->0.44781,u[4,2]->0.259156,u[4,3]->0.173101,u[4,4]->0.119934,v[1]->0.331787,v[2]->0.145926,v[3]->0.438083,v[4]->0.0842042}}

As you see they are all the same and probably the minimum to your function

Improve this answer
$\endgroup$
3
  • $\begingroup$ Cannot reproduce it even with NMinimize[{allM, cons}, vars, Method -> {"RandomSearch", "SearchPoints" -> 500}] on my comp during more than two hours. My best is NMinimize[{allM, cons}, vars, Method -> {"RandomSearch", "SearchPoints" -> 150}] // AbsoluteTiming which results in {1363.39, {5.39125, {u[1, 1] -> 0.335922, u[1, 2] -> 0.23733, u[1, 3] -> 0.21026, ...,, u[4, 4] -> 0.192082, v[1] -> 0.230135, v[2] -> 0.0337154, v[3] -> 0.258754, v[4] -> 0.477396}}} and several warnings about complex numbers. $\endgroup$
    – user64494
    Commented Apr 26, 2022 at 9:40
  • $\begingroup$ What version of Mathematica are you using? If I remember correctly, this optimization took a lot of time. I will look through my documents and post what I find. $\endgroup$
    – Tom6639
    Commented Apr 27, 2022 at 11:09
  • $\begingroup$ @user64494 I found where I solved this problem. Although I did not time the solution, I remember it took quite a while. But the solutions I presented above are those I got in v12.2. I suppose you will need to let you Mathematica run further to get the minimum I obtained. Alternatively, you can check that the solution I post is a minimum by creating a grid of search points with Method -> {"RandomSearch", "InitialPoints" ->searchpoints}. Where searchpoints is a list of range for the independent variables. It will take a lesser time to solve this way. $\endgroup$
    – Tom6639
    Commented Apr 27, 2022 at 14:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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