optimization - NMinimize 吃掉所有不必要的符号工作的内存 b/c

Question

下面的代码是一种简单的方法来找到平方有 n 个除数的最小数(最小值应该是它的对数，x_i 是它的素数分解的幂)。如果我查看 n=2000 的情况并使用 10 个变量而不是 20 个变量，这将使用大约 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}]。如果我只用 Product[Subscript[x, i],{i,20}] 替换产品并将约束更改为 >= 1 而不是 0 我可以轻松获得结果，并且内核使用超过 50MB 的内存。 (虽然这保留了不等式约束并且不会改变最小化目标函数的任务，但它确实弄乱了完整性要求 - 我得到了偶数结果，这对应于实际问题中的半整数。)

One person on StackOverflow有类似的问题；在他们的例子中，他们有一个目标函数，它以巨大的代价被象征性地评估。他们能够通过使函数只接受数字输入来解决这个问题，从而有效地将其隐藏在 Mathematica 的“我有 Expand[] 锤子，一切看起来都像钉子”的趋势中。但是你不能隐藏这样一个函数背后的约束(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 秒。

---结束编辑---

丹尼尔·利希特布劳