NMinimize吃掉不必要的符号工作的所有内存b/c [英] NMinimize eats all memory b/c of unnecessary symbolic work

问题描述

以下代码是一种幼稚的方法，可找到平方数为n的平方的最小数(最小应为对数，而x_i是其素因式分解的幂).如果我看n = 2000的情况，并使用10个变量而不是20个变量，那么这将使用大约600MB的内存.使用n的值，我实际上是在寻找答案.我需要大约20个变量，以确保不会丢失实际的解决方案，并且该解决方案很快会耗尽所有可用的内存，然后崩溃进行交换.

The following code is a naive way to find the least number whose square has n divisors (the minimum should be its log and the x_i the powers in its prime factorization). If I look at the case n=2000 and use ten variables instead of twenty, this uses somewhere around 600MB of memory. With the value of n I'm actually trying to find the answer for, I need around 20 variables to be sure of not missing the actual solution, and it quickly uses up all available memory and then thrashes swap.

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]

事实证明，问题根本与整数编程等无关(取消对整数的限制无济于事).

It turns out that the problem isn't related to integer programming etc at all (removing the restriction to the integers doesn't help).

我最好的猜测是问题在于Mathematica浪费了所有扩展Product[2 Subscript[x, i] + 1, {i, 20}]的内存.如果仅用Product[Subscript[x, i],{i,20}]替换产品，并将约束条件更改为>= 1而不是0，我将获得轻松的结果，并且内核不使用超过50MB的内存. (尽管它保留了不等式约束，并且没有改变使目标函数最小化的任务，但确实破坏了完整性要求-我得到了均匀的结果，对应于实际问题中的半整数.)

My best guess is that the problem is that Mathematica is wasting all that memory expanding Product[2 Subscript[x, i] + 1, {i, 20}]. If I replace the product with just Product[Subscript[x, i],{i,20}] and change the constraints to be >= 1 rather than 0 I get results without a hassle and without the kernel using more than 50MB of memory. (Though that preserves the inequality constraint and doesn't change the task of minimizing the objective function, it does mess up the integrality requirement- I get even results, which correspond to half-integers in the actual problem.)

StackOverflow上的一个人有类似的问题；在他们的案例中，他们有一个目标功能，该功能象征性地得到了巨大的代价.他们能够通过使该函数仅接受数字输入来对其进行补救，从而有效地将其隐藏在Mathematica的我有Expand []锤子，一切看上去像钉子"的趋势中.但是您不能将约束隐藏在这样的函数后面(Mathematica会抱怨这是无效的约束).

One person on StackOverflow had a similar problem; in their case, they had an objective function which was getting evaluated symbolically at a huge cost. They were able to remedy it by making the function only accept numeric input, effectively hiding it from Mathematica's "I have the Expand[] hammer, and everything looks like a nail" tendency. But you can't hide the constraint behind such a function (Mathematica will complain it's an invalid constraint).

关于如何解决此问题的任何想法?

Any thoughts on how to fix this?

我知道正确的答案-在我的Mathematica代码不起作用之后，我使用AMPL和专用的MINLP求解器来获取它(也很快).我只想知道我将来如何希望能够使用Mathematica的内置非线性优化功能，尽管当我以唯一的方式知道它们时，这似乎与我的约束有关. >

I know the correct answer- after my Mathematica code didn't work I used AMPL and a dedicated MINLP solver to get it (quite quickly too). I just want to know how I can ever hope to be able to use Mathematica's built-in nonlinear optimization features in the future despite the crazy things it seems to do with my constraints when I enter them in the only way I know how.

推荐答案

除非输入为数字，否则可以禁止该条件执行任何操作，如下所示.

Can inhibit that condition from doing anything unless inputs are numeric, as below.

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变量用于质数和幂的所有可能组合.

Thought I would point out that this can be set up as an integer linear programming problem. We use 0-1 variables for all possible combinations of primes and powers.

我们可以使用以下事实来限制素数的数量:该解决方案所涉及的素数不能超过所需的最小素数(假设所有素数都提高到第一乘方).如果目标从2开始连续，则目标是最小的.

We can limit the number of primes using the fact that the solution cannot involve more primes than the minimum needed assuming all are raised to the first power. The objective is then minimal if they are consecutive starting at 2.

我们将假定最大指数不超过20.可能有一种简便的方法可以显示这一点，但现在还没有想到.

We will assume the max exponent is no more than 20. There is probably a convenient way to show this but it has not come to mind as yet.

在此表述中，目标和制约因素如下.通过对除数sigma方程的对数得到一个完全线性的问题.

The objective and constraints, in this formulation, are as below. We get a fully linear problem by taking logs of the divisor sigma equation.

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

This approach handles n=10^10 is around 23 seconds.

-结束编辑---

丹尼尔·里奇布劳

这篇关于NMinimize吃掉不必要的符号工作的所有内存b/c的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！