2
$\begingroup$

I have the following code, which is a bunch of linear constraints and tries to maximize a nonlinear objective using Maximize.

bounds=Table[b[i], {i, 128}];
boundCond=And@@Map[GreaterEqual[#,0]&, bounds];
upperCond=0.100897647440433502 + 0.0435742177069187164 b[2] + 0.110877290368080139 b[3] + 0.0627215802669525147 b[6] + 0.0303984060883522034 b[11] + 0.195413455367088318 b[12] + 0.00218466715887188911 b[13] + 0.178133592009544373 b[14] + 0.114736616611480713 b[16] + 0.0892148315906524658 b[19] + 0.0430798530578613281 b[20] + 0.127856075763702393 b[22] + 0.235492616891860962 b[24] + 0.100893214344978332 b[26] + 0.0550539642572402954 b[29] + 0.0691445842385292053 b[30] + 0.220859974622726440 b[33] + 0.0824572741985321045 b[34] + 0.111420944333076477 b[35] + 0.119025722146034241 b[36] + 0.0626275166869163513 b[37] + 0.0163705591112375259 b[38] + 0.163605749607086182 b[39] + 0.0159755386412143707 b[40] + 0.222821488976478577 b[41] + 0.0262936670333147049 b[42] + 0.176039412617683411 b[44] + 0.119316443800926209 b[46] + 0.00839051231741905212 b[48] + 0.0780717059969902039 b[49] + 0.0129324616864323616 b[51] + 0.160396441817283630 b[54] + 0.122262813150882721 b[55] + 0.225195974111557007 b[58] + 0.00371941062621772289 b[61] + 0.0782109126448631287 b[66] + 0.163961067795753479 b[67] + 0.0488019287586212158 b[68] + 0.0772711709141731262 b[70] + 0.0967618897557258606 b[72] + 0.0753138065338134766 b[74] + 0.138780370354652405 b[75] + 0.0779474973678588867 b[77] + 0.143815591931343079 b[78] + 0.216688260436058044 b[79] + 0.192488059401512146 b[82] + 0.0594606176018714905 b[83] + 0.0763093978166580200 b[84] + 0.197820603847503662 b[86] + 0.0981343537569046020 b[88] + 0.177558749914169312 b[89] + 0.0137362014502286911 b[91] + 0.0624699629843235016 b[93] + 0.172200247645378113 b[94] + 0.0640604123473167419 b[97] + 0.253324717283248901 b[99] + 0.0138462502509355545 b[102] + 0.00301329372450709343 b[103] + 0.102329082787036896 b[105] + 0.0432044677436351776 b[107] + 0.0967931896448135376 b[109] + 0.00924742128700017929 b[110] + 0.194597169756889343 b[112] + 0.0199854951351881027 b[115] + 0.00540251936763525009 b[117] + 0.0856299698352813721 b[118] + 0.00201270543038845062 b[119] + 0.0697343200445175171 b[120] + 0.0598938837647438049 b[121] + 0.115721419453620911 b[123] + 0.134053930640220642 b[126] + 0.0757555514574050903 b[128] <= 1.4
product:=Apply[Times, bounds];
Maximize[{product, upperCond && boundCond}, bounds];

Basically, bounds is a list of positive variables, and there is a linear constraint upperCond that involves some of the bounds variables, but not all. And I want to maximize the product of all bounds variables.

Now intuitively the maximal value of product should be unbounded because b[1], for example, could just be infinity as it does not appear in the linear constraint. But Mathematica gave me the following results:

Out[119]= {1.0149 10^-53, {b[1] -> 0.631844, b[2] -> 0.366275, b[3] -> 0.536577, b[4] -> 1.06186, b[5] -> 1.14916, b[6] -> 0.307205, b[7] -> 1.02957, b[8] -> 0.590336, b[9] -> 1.02625, b[10] -> 0.781978, 
b[11] -> 0.49093, b[12] -> 0.0258683, b[13] -> 0.585934, b[14] -> 0.0313637, b[15] -> 1.22059, b[16] -> 0.0865507, b[17] -> 0.744385, b[18] -> 0.652376, b[19] -> 0.191569, b[20] -> 0.401325, 
b[21] -> 0.591891, b[22] -> 0.0787587, b[23] -> 1.1079, b[24] -> 0.0211321, b[25] -> 0.673995, b[26] -> 0.270791, b[27] -> 1.09818, b[28] -> 1.15599, b[29] -> 0.847054, b[30] -> 0.286807, 
b[31] -> 1.12585, b[32] -> 1.12369, b[33] -> 0.0213811, b[34] -> 0.214411, b[35] -> 0.156409, b[36] -> 0.103377, b[37] -> 0.748765, b[38] -> 0.910576, b[39] -> 0.0460899, b[40] -> 0.943243, 
b[41] -> 0.0209113, b[42] -> 0.825903, b[43] -> 0.987693, b[44] -> 0.026583, b[45] -> 1.10995, b[46] -> 0.155892, b[47] -> 0.687417, b[48] -> 0.987659, b[49] -> 0.36043, b[50] -> 0.722281, 
b[51] -> 0.551261, b[52] -> 0.733565, b[53] -> 1.09543, b[54] -> 0.0416588, b[55] -> 0.596005, b[56] -> 0.599603, b[57] -> 1.16556, b[58] -> 0.0485591, b[59] -> 1.03535, b[60] -> 1.03462, 
b[61] -> 0.631659, b[62] -> 1.05244, b[63] -> 1.15654, b[64] -> 0.734218, b[65] -> 0.951136, b[66] -> 0.650616, b[67] -> 0.0792143, b[68] -> 0.354095, b[69] -> 0.559035, b[70] -> 0.295494, 
b[71] -> 0.627549, b[72] -> 0.153516, b[73] -> 0.578213, b[74] -> 0.375231, b[75] -> 0.0722742, b[76] -> 0.980545, b[77] -> 0.293393, b[78] -> 0.0657513, b[79] -> 0.028435, b[80] -> 1.06805, 
b[81] -> 1.03576, b[82] -> 0.0495011, b[83] -> 0.797685, b[84] -> 0.262131, b[85] -> 0.704682, b[86] -> 0.0488136, b[87] -> 0.730814, b[88] -> 0.129551, b[89] -> 0.0569268, b[90] -> 0.568371, 
b[91] -> 0.954501, b[92] -> 0.597521, b[93] -> 0.294341, b[94] -> 0.0337064, b[95] -> 1.1071, b[96] -> 0.97234, b[97] -> 0.254995, b[98] -> 0.71303, b[99] -> 0.0327845, b[100] -> 1.10424, 
b[101] -> 1.17713, b[102] -> 0.931491, b[103] -> 0.604491, b[104] -> 0.69352, b[105] -> 0.146641, b[106] -> 0.719575, b[107] -> 0.471206, b[108] -> 1.15883, b[109] -> 0.203697, b[110] -> 0.560459, 
b[111] -> 0.702077, b[112] -> 0.0634536, b[113] -> 1.11802, b[114] -> 1.08117, b[115] -> 0.583974, b[116] -> 0.987578, b[117] -> 0.655513, b[118] -> 0.214568, b[119] -> 0.652192, 
b[120] -> 0.426065, b[121] -> 0.307841, b[122] -> 0.996165, b[123] -> 0.130277, b[124] -> 1.02606, b[125] -> 1.17633, b[126] -> 0.0679813, b[127] -> 1.09422, b[128] -> 0.612329}}

Clearly Mathematica assigned a value to b[1], among other bounds variables that should have been unbounded, which leads to a "maximal" solution which should not exist.

But if I change the optimize objective from the product of all the bounds variables to the sum of all the bounds variables, it returned: NMaximize::ubnd: The problem is unbounded., which makes sense.

And also if I limit the number of bounds variables to say 5, it returned NMaximize::cvdiv: Failed to converge to a solution. The function may be unbounded., which again makes sense.

So seems like Maximize automatically uses NMaximize if constraints involve approximate real numbers, which clearly is the case here. So is the rounding error in NMaximize causing the problem when trying to find a numerical solution? And why is the rounding error not showing up in optimize for the sum of all the variables, and not showing up when the number of variables is small?

Improve this question
$\endgroup$
4
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) 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. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Aug 22, 2016 at 0:14
  • $\begingroup$ Perhaps, the computation does not have enough precision. $\endgroup$
    – bbgodfrey
    Aug 22, 2016 at 0:40
  • $\begingroup$ Actually all numbers' precisions are beyond the machine precision. So my understanding is that Mathematica would treat them as arbitrary‐precision numbers, and thus preserve the precision. I also tried to explicitly set the precision to be say 18, which is higher than machine precision, and ended up with the same result. $\endgroup$
    – queeten
    Aug 22, 2016 at 1:30
  • $\begingroup$ Not so. It treats them as precisely the accuracy and precision they have. See my answer below. It is not uncommon for inadequate precision to cause problems. $\endgroup$
    – bbgodfrey
    Aug 22, 2016 at 1:54

1 Answer 1

Reset to default
1
$\begingroup$

For a small number of b, say 23, Maximize returns the error message

Maximize::natt: The maximum is not attained at any point satisfying the given constraints.

as it should. However, for larger numbers of b it simply returns an incorrect answer, as noted in the question. This can be corrected by using Rationalize for upperCond and higher WorkingPrecision for Maximize. For instance, with 128 b, as in the question, use

upperCond = Rationalize[0.100897647440433502 + 0.0435742177069187164 b[2] + ..., 10^-18]

and

s = Maximize[{product, upperCond && boundCond}, bounds, WorkingPrecision -> 120]

which now also returns the error message above.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks for the answer! Seems like forcing Mathematica to work with precise numbers (through Rationalize) did work. But I am still a little confused. I thought the internal algorithms that the Wolfram Language uses to evaluate mathematical functions are set up to maintain as much precision as possible to the point that it ensures that all digits returned by the Wolfram Language are correct. So if that's the case, we should be able to trust the returned results at least to the level of precision that's contained in the results right? $\endgroup$
    – queeten
    Aug 22, 2016 at 21:49
  • $\begingroup$ Unfortunately, no. It maintains as much precision as possible within the constraint of the WorkingPrecision option, which defaults to machine precision (just less than 16 digits) To get more precision from a function, use a larger WorkingPrecision. However, Mathematica complains if you give it numbers with precision less than the WorkingPrecision you specify. That is why I used Rationalize. By the way, thanks for accepting my answer. Best wishes. $\endgroup$
    – bbgodfrey
    Aug 22, 2016 at 21:54
  • $\begingroup$ Your explanation makes sense, but I was still left wondering why Mathematica needed a higher precision to get a correct answer? That is, what exactly went wrong in the original calculation with lower precision? $\endgroup$
    – queeten
    Aug 22, 2016 at 22:37
  • $\begingroup$ It is difficult to know without understanding the details of the internal working of Maximize. Your problem seemed like a precision issue to me, so I explored that possibility. I am not sure how I would have proceeded, if my guess had been wrong. $\endgroup$
    – bbgodfrey
    Aug 22, 2016 at 23:02

Your Answer

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

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