NMinimize пожирает всю память ч / б ненужной символической работы - PullRequest
Новая
       25

NMinimize пожирает всю память ч / б ненужной символической работы

3 голосов
/ 07 августа 2011

Следующий код является наивным способом найти наименьшее число, у квадрата которого есть n делителей (минимум должен быть его логарифмом, а x_i - степенями в его первичной факторизации).Если я смотрю на случай n = 2000 и использую десять переменных вместо двадцати, это использует где-то около 600 МБ памяти.При значении n, на самом деле я пытаюсь найти ответ, мне нужно около 20 переменных, чтобы быть уверенным в том, что оно не пропустило реальное решение, и оно быстро использует всю доступную память, а затем сбрасывает своп.

Оказывается, что проблема вообще не связана с целочисленным программированием и т. Д. (Снятие ограничения на целые числа не помогает).

Мое лучшее предположение состоит в том, что проблема в том, что Mathematica тратит впустуювся эта память расширяется Product[2 Subscript[x, i] + 1, {i, 20}].Если я заменю продукт просто на Product[Subscript[x, i],{i,20}] и поменяю ограничения на >= 1 вместо 0, я получу результаты без хлопот и без ядра, использующего более 50 МБ памяти.(Хотя это сохраняет ограничение неравенства и не меняет задачу минимизации целевой функции, оно портит требование целостности - я получаю четные результаты, которые соответствуют половинным целым числам в актуальной задаче.)

У одного человека в StackOverflow была похожая проблема;в их случае у них была целевая функция, которая оценивалась символически с огромными затратами.Они смогли исправить это, заставив функцию принимать только числовой ввод, эффективно скрывая ее от стремления Mathematica «У меня есть молоток Expand [], и все выглядит как гвоздь».Но вы не можете скрыть ограничение за такой функцией (Mathematica будет жаловаться, что это недопустимое ограничение).

Есть мысли о том, как это исправить?

Редактировать: я знаю правильный ответ - после того, как мой код Mathematica не работал, я использовал AMPL и специальный решатель MINLP, чтобы получить его (тоже довольно быстро).Я просто хочу знать, как я могу надеяться когда-либо использовать встроенные в Mathematica функции нелинейной оптимизации в будущем, несмотря на сумасшедшие вещи, которые, как мне кажется, связаны с моими ограничениями, когда я вхожу в них только так, как я знаю.

Ответы [ 2 ]

7 голосов
/ 08 августа 2011

Может помешать этому условию ничего не делать, если входные данные не являются числовыми, как показано ниже. 

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

Думаю, я бы указал, что это можно установитькак целочисленная задача линейного программирования.Мы используем 0-1 переменные для всех возможных комбинаций простых чисел и степеней.

Мы можем ограничить число простых чисел, используя тот факт, что решение не может включать больше простых чисел, чем необходимый минимум, при условии, что все возведены в первую степень,Тогда цель минимальна, если они последовательно начинаются с 2.

Мы будем предполагать, что максимальный показатель не превышает 20. Вероятно, есть удобный способ показать это, но он еще не пришел в голову.

Цель и ограничения в этой формулировке следующие.Мы получаем полностью линейную задачу, беря логарифмы уравнения делителя сигмы. 

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 секунд.

--- end edit --

Даниэль Лихтблау

2 голосов
/ 07 августа 2011

Вы можете попробовать этот код вместо: 

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

, который возвращает 180180.

Сложение: 

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

Кажется, работает быстрее.

ДОПОЛНЕНИЕ 2:

Вот для этой ультра-улучшенной версии (но не доказано на 100%): 

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] дает 2768774904222066200260800 через ~ 4 секунды

Объяснение

Прежде всего необходимо доказать, что сумма показателей простых чисел в этом минимальном числе не более Log[2, n]. Я еще не нашел доказательства , но это может быть связано с соотношением последовательных простых чисел.

Flatten[IntegerPartitions /@ Range[Ceiling[Log[2, n]]], 1] дает вам все списки с Total <= Log[2, n], удобно отсортированные от большого к маленькому.

Product[Prime[k]^#[[k]], {k, 1, Length[#]}] & использовать их в качестве простых чисел для создания целых чисел.

Min@Select[..., DivisorSigma[0, #^2] > n &] выберите минимальное из них, соответствующее исходному условию.

...