NMinimize eats all memory b/c of unnecessary symbolic work

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

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

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 makin开发者_开发百科g 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?

Edit: 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}}}

---edit---

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.

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.

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.

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}}

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

---end edit ---

Daniel Lichtblau

You can try this code instead:

Catch[Do[If[DivisorSigma[0, k^2] > 2000, Throw[k]], {k, 1000000}]]

which returns 180180.

ADDITION:

Catch[Do[If[Times@@(2 FactorInteger[k][[All, 2]] + 1) > 2000, Throw[k]], {k, 1000000}]]

Seems to work faster.

ADDITION 2:

Behold for this ultra-improved version (but not 100% proved):

MinSquareWithDivisors[n_] :=
 Min@Select[
   Product[Prime[k]^#[[k]], {k, 1, Length[#]}] & /@ 
    Flatten[IntegerPartitions /@ Range[Ceiling[Log[2, n]]], 1], 
   DivisorSigma[0, #^2] > n &]

MinSquareWithDivisors[2000000000] gives 2768774904222066200260800 in ~4 seconds

Explanation:

First of all one needs to prove that the sum of the prime's exponents in this minimum number is at most Log[2, n]. I haven't found a proof yet, but it might be related to the ratio between successive primes.

Flatten[IntegerPartitions /@ Range[Ceiling[Log[2, n]]], 1] gives you all the lists with Total <= Log[2, n], conveniently sorted from large to small.

Product[Prime[k]^#[[k]], {k, 1, Length[#]}] & use them as prime's exponents to create integers.

Min@Select[..., DivisorSigma[0, #^2] > n &] choose the minimal of them that agrees with the original condition.