Существует ли быстрая операция с продуктом для PackedArrays? - PullRequest
Новая
       116

Существует ли быстрая операция с продуктом для PackedArrays?

8 голосов
/ 14 марта 2011

В Mathematica вектор (или прямоугольный массив), содержащий все целые числа или числа с плавающей запятой, может храниться в упакованном массиве.Эти объекты занимают меньше памяти, а некоторые операции с ними выполняются намного быстрее.

RandomReal по возможности создает упакованный массив.Упакованный массив может быть распакован с помощью функции Developer FromPackedArray

. Рассмотрим эти временные параметры 

lst = RandomReal[1, 5000000];

Total[lst] // Timing
Plus @@ lst // Timing

lst = Developer`FromPackedArray[lst];

Total[lst] // Timing
Plus @@ lst // Timing

Out[1]= {0.016, 2.50056*10^6}

Out[2]= {0.859, 2.50056*10^6}

Out[3]= {0.625, 2.50056*10^6}

Out[4]= {0.64, 2.50056*10^6}

Следовательно, в случае упакованного массива Total во много раз быстреечем Plus @@, но примерно так же для неупакованного массива.Обратите внимание, что Plus @@ на самом деле немного медленнее в упакованном массиве.

Теперь рассмотрим 

lst = RandomReal[100, 5000000];
Times @@ lst // Timing

lst = Developer`FromPackedArray[lst];
Times @@ lst // Timing

Out[1]= {0.875, 5.8324791357*10^7828854}

Out[1]= {0.625, 5.8324791357*10^7828854}

Наконец, мой вопрос: есть ли быстрый метод в Mathematica для списка продуктовупакованный массив, аналогичный Total?

Я подозреваю, что это может быть невозможно из-за того, что числовые ошибки сочетаются с умножением.Кроме того, функция должна иметь возможность возвращать не машинные значения с плавающей точкой, чтобы быть полезной.

Ответы [ 3 ]

9 голосов
/ 14 марта 2011

Я также задавался вопросом, был ли мультипликативный эквивалент Total.

A действительно не так уж и плохое решение 

In[1]:= lst=RandomReal[2,5000000];
        Times@@lst//Timing
        Exp[Total[Log[lst]]]//Timing
Out[2]= {2.54,4.370467929041*10^-666614}
Out[3]= {0.47,4.370467940*10^-666614}

Пока числа положительные и не слишком большие или маленькие, ошибки округления не будутТоже плохо.Предположение относительно того, что может происходить во время оценки, заключается в следующем: (1) Если числа являются положительными числами, операция Log может быть быстро применена к упакованному массиву.(2) Числа могут быть быстро добавлены с помощью метода упакованного массива Total.(3) Тогда это только последний шаг, когда возникает необходимость в плавании не машинного размера.

См. этот ответ SO для решения, которое работает как для положительных, так и для отрицательных чисел.

Давайте быстро проверим, что это решение работает с числами с плавающей точкой, которые дают ответ не машинного размера.Сравните с Эндрю (намного быстрее) compiledListProduct: 

In[10]:= compiledListProduct = 
          Compile[{{l, _Real, 1}}, 
           Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot], 
           CompilationTarget -> "C"]

In[11]:= lst=RandomReal[{0.05,.10},15000000];
         Times@@lst//Timing
         Exp[Total[Log[lst]]]//Timing
         compiledListProduct[lst]//Timing
Out[12]= {7.49,2.49105025389*10^-16998863}
Out[13]= {0.5,2.4910349*10^-16998863}
Out[14]= {0.07,0.}

Если вы выберете более крупные (>1) реалы, то compiledListProduct выдаст предупреждение CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. и потребуется некоторое время, чтобы датьрезультат ...

Любопытно, что и Sum, и Product могут принимать произвольные списки.Sum отлично работает 

In[4]:= lst=RandomReal[2,5000000];
         Sum[i,{i,lst}]//Timing
         Total[lst]//Timing
Out[5]= {0.58,5.00039*10^6}
Out[6]= {0.02,5.00039*10^6}

, но для длинных PackedArray с, таких как тестовые примеры здесь, Product завершается неудачно, поскольку автоматически скомпилированный код (в версии 8.0) не улавливает неполадки / переполнения должным образом: 

In[7]:= lst=RandomReal[2,5000000];
         Product[i,{i,lst}]//Timing
         Times@@lst//Timing
Out[8]= {0.,Compile`AutoVar12!}
Out[9]= {2.52,1.781498881673*10^-666005}

Обход, предоставляемый службой технической поддержки WRI, состоит в том, чтобы отключить компиляцию продукта с помощью SetSystemOptions["CompileOptions" -> {"ProductCompileLength" -> Infinity}].Другой вариант - использовать lst=Developer`FromPackedArray[lst].

4 голосов
/ 14 марта 2011

Во-первых, чтобы избежать путаницы, взгляните на пример, все результаты которого представлены в виде чисел точности аппаратного обеспечения, которые должны быть меньше 

In[1]:= $MaxMachineNumber

Out[1]= 1.79769*10^308

Ваш общий пример уже имел это хорошее (и быстрое) свойство. Вот пример на вашем примере Times с использованием машинных номеров: 

In[2]:= lst = RandomReal[{0.99, 1.01}, 5000000];
Times @@ lst // Timing

Out[3]= {1.435, 1.38851*10^-38}

Теперь мы можем использовать Compile , чтобы сделать скомпилированную функцию для эффективного выполнения этой операции: 

In[4]:= listproduct = 
 Compile[{{l, _Real, 1}}, 
  Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot]]

Out[4]= CompiledFunction[{l},Module[{tot=1.},Do[tot*=x,{x,l}];tot],-CompiledCode-]

Это намного быстрее: 

In[5]:= listproduct[lst] // Timing

Out[5]= {0.141, 1.38851*10^-38}

Предполагая, что у вас есть компилятор C и Mathematica 8, вы также можете автоматически скомпилировать весь код на C. Временная DLL создается и связывается обратно с Mathematica во время выполнения. 

In[6]:= compiledlistproduct = 
 Compile[{{l, _Real, 1}}, 
  Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot], 
  CompilationTarget -> "C"]

Out[6]= CompiledFunction[{l},Module[{tot=1.},Do[tot*=x,{x,l}];tot],-CompiledCode-]

Это дает производительность, не сильно отличающуюся от той, которая была бы у встроенной функции Mathematica: 

In[7]:= compiledlistproduct[lst] // Timing

Out[7]= {0.015, 1.38851*10^-38}

Обратите внимание, что если ваш продукт действительно выйдет за пределы $ 1021 * MaxMachineNumber (или $ MinMachineNumber ), то вам лучше придерживаться Apply[Times, list]. Тот же комментарий относится и к Total, если ваши результаты могут стать такими большими: 

In[11]:= lst = RandomReal[10^305, 5000000];
Plus @@ lst // Timing

Out[12]= {1.435, 2.499873364498981*10^311}

In[13]:= lst = RandomReal[10^305, 5000000];
Total[lst] // Timing

Out[14]= {1.576, 2.500061580905602*10^311}
3 голосов
/ 15 марта 2011

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

Функция 

f = (-1)^(-1 /. Rule @@@ Tally@Sign@# /. -1 -> 0) * Exp@Total@Log@Abs@# &;

Тестирование 

lst = RandomReal[{-50, 50}, 5000000];

Times @@ lst // Timing
f@lst // Timing

lst = Developer`FromPackedArray[lst];

Times @@ lst // Timing
f@lst // Timing

{0.844, -4.42943661963*10^6323240}

{0.062, -4.4294366*10^6323240}

{0.64, -4.42943661963*10^6323240}

{0.203, -4.4294366*10^6323240}
...