NMinimize吃掉所有不必要的符号工作的记忆

Question

下面的代码是一种天真的方法来查找其正方形具有n个除数的最小数字(最小值应为其对数,x_i为其素数因子分解中的幂).如果我看一下n = 2000的情况并且使用十个变量而不是二十个变量,那么这将使用大约600MB的内存.使用n的值我实际上试图找到答案,我需要大约20个变量以确保不会错过实际的解决方案,并且它会快速耗尽所有可用内存然后打开交换.

n=8*10^6;
a = Table[N[Log[Prime[i]]], {i, 20}];
b = Table[Subscript[x, i], {i, 20}];
cond = Fold[And, Product[2 Subscript[x, i] + 1, {i, 20}] > n,
   Table[Subscript[x, i] >= 0, {i, 20}]] && b \[Element] Integers;
NMinimize[{a.b, cond}, b, MaxIterations -> 1000]

事实证明,问题与整数编程等无关(取消对整数的限制没有帮助).

我最好的猜测是问题是Mathematica浪费了所有内存扩展Product[2 Subscript[x, i] + 1, {i, 20}].如果我用just替换产品Product[Subscript[x, i],{i,20}]并更改约束>= 1而不是0我没有麻烦得到结果,并且没有内核使用超过50MB的内存.(虽然这保留了不等式约束并且没有改变最小化目标函数的任务,但它确实搞乱了完整性要求 - 我得到了甚至结果,这对应于实际问题中的半整数.)

StackOverflow上的一个人有类似的问题; 在他们的情况下,他们有一个客观的功能,正在以巨大的代价象征性地进行评估.他们能够通过使函数只接受数字输入来有效地修复它,从而有效地将它隐藏在Mathematica的"我有扩展[]锤子之后,一切看起来像钉子"的倾向.但是你不能隐藏这种函数背后的约束(Mathematica会抱怨它是一个无效的约束).

有关如何解决此问题的任何想法？

编辑:我知道正确的答案 - 在我的Mathematica代码不起作用之后我使用了AMPL和一个专用的MINLP解算器来获得它(很快).我只是想知道我怎么能希望能够在将来使用Mathematica的内置非线性优化功能,尽管当我以我知道的唯一方式进入它们时,它似乎与我的约束有关.

Answer 1

除非输入是数字,否则可以禁止该条件做任何事情,如下所示.

n = 8*10^6;
nvars = 20;
a = Table[N[Log[Prime[i]]], {i, nvars}];
b = Table[Subscript[x, i], {i, nvars}];
c1[xx : {_?NumericQ ..}] := Times @@ (2 xx + 1);
c2 = Map[# >= 0 &, b];
c3 = b \[Element] Integers;
cond = Join[{c1[b] > n}, c2, {c3}];

In[29]:= Timing[NMinimize[{a.b, cond}, b, MaxIterations -> 400]]

Out[29]= {43.82100000000008, {36.77416664719056, {Subscript[x, 1] -> 
    3, Subscript[x, 2] -> 3, Subscript[x, 3] -> 2, 
   Subscript[x, 4] -> 2, Subscript[x, 5] -> 1, Subscript[x, 6] -> 1, 
   Subscript[x, 7] -> 1, Subscript[x, 8] -> 1, Subscript[x, 9] -> 1, 
   Subscript[x, 10] -> 1, Subscript[x, 11] -> 1, 
   Subscript[x, 12] -> 1, Subscript[x, 13] -> 0, 
   Subscript[x, 14] -> 0, Subscript[x, 15] -> 0, 
   Subscript[x, 16] -> 0, Subscript[x, 17] -> 0, 
   Subscript[x, 18] -> 0, Subscript[x, 19] -> 0, 
   Subscript[x, 20] -> 0}}}

- -编辑 - -

我想我会指出这可以设置为整数线性规划问题.我们使用0-1变量来表示素数和幂的所有可能组合.

我们可以使用以下事实来限制质数的数量:假设所有质数都提高到第一个幂,解决方案不会涉及比最小需求更多的素数.如果从2开始连续,则目标是最小的.

我们假设最大指数不超过20.可能有一种方便的方式来显示它,但它还没有想到.

在这个表述中,目标和约束如下.我们通过取除除数sigma方程的对数得到一个完全线性的问题.

n = 8*10^6;
nprimes = Ceiling[Log[2, n]];
maxexpon = 20;
vars = Array[x, {maxexpon, nprimes}];
fvars = Flatten[vars];
c1 = Map[0 <= # <= 1 &, fvars];
c2 = {Element[fvars, Integers]};
c3 = Thread[Total[vars] <= 1];
c4 = {Total[N[Log[2*Range[maxexpon] + 1]].vars] >= N@Log[n]};
constraints = Join[c1, c2, c3, c4];
obj = Range[maxexpon].vars.N[Log[Prime[Range[nprimes]]]];

Timing[{min, vals} = NMinimize[{obj, constraints}, fvars];]

Out[118]= {5.521999999999991, Null}

Pick[fvars, fvars /. vals, 1] /. x[j_, k_] :> {Prime[k], j}

Out[119]= {{11, 1}, {13, 1}, {17, 1}, {19, 1}, {23, 1}, {29, 1}, {31, 
  1}, {37, 1}, {5, 2}, {7, 2}, {2, 3}, {3, 3}}

这种方法处理n = 10 ^ 10约为23秒.

---结束编辑---

Daniel Lichtblau