Pregunta ¿Hay una operación de producto rápida para PackedArrays?


En Mathematica, un vector (o matriz rectangular) que contiene todos los enteros o flotantes del tamaño de la máquina puede almacenarse en una matriz empaquetada. Estos objetos consumen menos memoria y algunas operaciones son mucho más rápidas.

RandomReal produce una matriz empaquetada cuando es posible. Una matriz empaquetada se puede desempaquetar con Developer función FromPackedArray

Considere estos tiempos

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}

Por lo tanto, en el caso de una matriz empaquetada, Total es mucho más rápido que Plus @@ pero más o menos lo mismo para una matriz no empaquetada. Tenga en cuenta que Plus @@ es en realidad un poco más lento en la matriz empaquetada.

Ahora considera

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}

Finalmente, mi pregunta: ¿existe un método rápido en Mathematica para el producto de la lista de una matriz empaquetada, análoga a Total?

Sospecho que esto puede no ser posible debido a la forma en que los errores numéricos se combinan con la multiplicación. Además, la función deberá ser capaz de devolver flotadores que no sean de máquina para que sea útil.


8
2018-03-14 04:12


origen


Respuestas:


También me he preguntado si había un multiplicativo equivalente a Total.

UN De Verdad no tan mala solución es

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}

Siempre que los números sean positivos y no sean demasiado grandes o pequeños, entonces los errores de redondeo no son tan malos. Una suposición acerca de lo que podría estar sucediendo durante la evaluación es que: (1) siempre que los números sean flotantes positivos, Log la operación se puede aplicar rápidamente a la matriz empaquetada. (2) Los números se pueden agregar rápidamente usando Totalestá lleno de método de matriz. (3) Entonces es solo el último paso donde se necesita un flotador sin tamaño de máquina.

Ver esta respuesta SO para una solución que funcione para flotadores positivos y negativos.

Comprobemos rápidamente que esta solución funciona con flotadores que ofrecen una respuesta que no es de tamaño de máquina. Comparar con Andrew's (mucho más rápido) 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.}

Si elige más grande (>1) reales, luego compiledListProduct rendirá la advertencia CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. y tomará un tiempo para dar un resultado ...


Una currículum es que ambos Sum y Product puede tomar listas arbitrarias Sum funciona bien

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}

pero por mucho tiempo PackedArrays, como los ejemplos de prueba aquí, Product falla ya que el código compilado automáticamente (en la versión 8.0) no captura los subdesbordamientos / desbordamientos correctamente:

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}

El trabajo de reparación proporcionado por el útil soporte técnico de WRI consiste en desactivar la compilación del producto utilizando SetSystemOptions["CompileOptions" -> {"ProductCompileLength" -> Infinity}]. Otra opción es usar lst=Developer`FromPackedArray[lst].


9
2018-03-14 05:05



Primero, para evitar confusiones, eche un vistazo a un ejemplo cuyos resultados son representables como números de precisión de la máquina de hardware, que deben ser todos menos de

In[1]:= $MaxMachineNumber

Out[1]= 1.79769*10^308

Su ejemplo Total ya tenía esta propiedad agradable (y rápida). Aquí hay una variante en su ejemplo de Times usando números de máquina:

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

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

Ahora podemos usar Compilar para hacer una función compilada para realizar esta operación de manera eficiente:

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

Es mucho más rápido:

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

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

Suponiendo que tiene un compilador de C y Mathematica 8, también puede compilar automáticamente todo el camino hasta el código C. Se crea una DLL temporal y se vincula nuevamente a Mathematica en tiempo de ejecución.

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

Esto le da un rendimiento no muy diferente al que tendría una función incorporada de Mathematica:

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

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

Tenga en cuenta que si su producto realmente va más allá de $MaxMachineNumber (o $ MinMachineNumber), entonces es mejor quedarse con Apply[Times, list]. El mismo comentario se aplica a Total, si tus resultados pueden ser tan grandes:

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}

4
2018-03-14 05:30



El método de Simon es rápido, pero falla en valores negativos. Combinándolo con su respuesta a mi otra pregunta, aquí hay una solución rápida que maneja los negativos. Gracias, Simon.

Función

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

Pruebas

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}

3
2018-03-15 06:53