Cómo acelerar múltiples productos internos en python
performance algorithm (5)
Tengo un código simple que hace lo siguiente.
Se itera sobre toda la longitud posible n enumera F con + -1 entradas. Para cada uno, itera sobre toda la longitud posible 2n listas S con + -1 entradas, donde la primera mitad de $ S $ es simplemente una copia de la segunda mitad. El código calcula el producto interno de F con cada sub-lista de S de longitud n . Para cada F, S cuenta los productos internos que son cero hasta el primer producto interno que no es cero.
Aquí está el código.
#!/usr/bin/python
from __future__ import division
import itertools
import operator
import math
n=14
m=n+1
def innerproduct(A, B):
assert (len(A) == len(B))
s = 0
for k in xrange(0,n):
s+=A[k]*B[k]
return s
leadingzerocounts = [0]*m
for S in itertools.product([-1,1], repeat = n):
S1 = S + S
for F in itertools.product([-1,1], repeat = n):
i = 0
while (i<m):
ip = innerproduct(F, S1[i:i+n])
if (ip == 0):
leadingzerocounts[i] +=1
i+=1
else:
break
print leadingzerocounts
La salida correcta para n=14 es
[56229888, 23557248, 9903104, 4160640, 1758240, 755392, 344800, 172320, 101312, 75776, 65696, 61216, 59200, 59200, 59200]
Usando pypy, esto toma 1 min 18 segundos para n = 14. Desafortunadamente, realmente me gustaría ejecutarlo para 16,18,20,22,24,26. No me importa usar numba o cython, pero me gustaría estar cerca de python si es posible.
Cualquier ayuda para acelerar esto es muy apreciada.
Mantendré un registro de las soluciones más rápidas aquí. (Por favor, avíseme si pierdo una respuesta actualizada).
- n = 22 a 9m35.081s por Eisenstat (C)
- n = 18 a 1m16.344s por Eisenstat (pypy)
- n = 18 a 2m54.998s por Tupteq (pypy)
- n = 14 en 26s por Neil (numpy)
- n - 14 a 11m59.192s por kslote1 (pypy)
En mi opinión, una buena manera de obtener un aumento de rendimiento es utilizar los elementos incorporados de python.
Primer mapa de uso para calcular el producto de las entradas:
>>> a =[1,2,3]
>>> b = [4,5,6]
>>>map(lambda x,y : x*y, a , b)
[4, 10, 18]
Luego use reduce para calcular las sumas:
>>> reduce(lambda v,w: v+w, map(lambda x,y :x*y, a, b))
32
Entonces tu función se convierte en
def innerproduct(A, B):
assert (len(A) == len(B))
return reduce(lambda v,w: v+w, map(lambda x,y :x*y, A, B))
A continuación, podemos eliminar todos esos "para bucles" y reemplazarlos con generadores y capturar el StopIteration.
#!/usr/bin/python
from __future__ import division
import itertools
import operator
import math
n=14
m=n+1
def innerproduct(A, B):
assert (len(A) == len(B))
return reduce(lambda v,w: v+w, map(lambda x,y :x*y, A, B))
leadingzerocounts = [0]*m
S_gen = itertools.product([-1,1], repeat = n)
try:
while(True):
S = S_gen.next()
S1 = S + S
F_gen = itertools.product([-1,1], repeat = n)
try:
while(True):
F = F_gen.next()
for i in xrange(m):
ip = innerproduct(F, S1[i:i+n])
if (ip == 0):
leadingzerocounts[i] +=1
i+=1
else:
break
except StopIteration:
pass
except StopIteration as e:
print e
print leadingzerocounts
Observé una aceleración para n menor, pero mi jalopy carecía del poder de caballo para calcular mi versión ni el código original para n = 14. Una forma de acelerar esto sería memorizar la línea:
F_gen = itertools.product([-1,1], repeat = n)
Este nuevo código se acelera en otro orden de magnitud al aprovechar la simetría cíclica del problema. Esta versión de Python enumera collares con el algoritmo de Duval; La versión C usa fuerza bruta. Ambos incorporan las aceleraciones que se describen a continuación. En mi máquina, la versión C resuelve n = 20 en 100 segundos! Un cálculo en la parte posterior del sobre sugiere que, si lo dejara correr durante una semana en un solo núcleo, podría hacer n = 26 y, como se indica a continuación, es susceptible al paralelismo.
import itertools
def necklaces_with_multiplicity(n):
assert isinstance(n, int)
assert n > 0
w = [1] * n
i = 1
while True:
if n % i == 0:
s = sum(w)
if s > 0:
yield (tuple(w), i * 2)
elif s == 0:
yield (tuple(w), i)
i = n - 1
while w[i] == -1:
if i == 0:
return
i -= 1
w[i] = -1
i += 1
for j in range(n - i):
w[i + j] = w[j]
def leading_zero_counts(n):
assert isinstance(n, int)
assert n > 0
assert n % 2 == 0
counts = [0] * n
necklaces = list(necklaces_with_multiplicity(n))
for combo in itertools.combinations(range(n - 1), n // 2):
for v, multiplicity in necklaces:
w = list(v)
for j in combo:
w[j] *= -1
for i in range(n):
counts[i] += multiplicity * 2
product = 0
for j in range(n):
product += v[j - (i + 1)] * w[j]
if product != 0:
break
return counts
if __name__ == ''__main__'':
print(leading_zero_counts(12))
Versión C:
#include <stdio.h>
enum {
N = 14
};
struct Necklace {
unsigned int v;
int multiplicity;
};
static struct Necklace g_necklace[1 << (N - 1)];
static int g_necklace_count;
static void initialize_necklace(void) {
g_necklace_count = 0;
for (unsigned int v = 0; v < (1U << (N - 1)); v++) {
int multiplicity;
unsigned int w = v;
for (multiplicity = 2; multiplicity < 2 * N; multiplicity += 2) {
w = ((w & 1) << (N - 1)) | (w >> 1);
unsigned int x = w ^ ((1U << N) - 1);
if (w < v || x < v) goto nope;
if (w == v || x == v) break;
}
g_necklace[g_necklace_count].v = v;
g_necklace[g_necklace_count].multiplicity = multiplicity;
g_necklace_count++;
nope:
;
}
}
int main(void) {
initialize_necklace();
long long leading_zero_count[N + 1];
for (int i = 0; i < N + 1; i++) leading_zero_count[i] = 0;
for (unsigned int v_xor_w = 0; v_xor_w < (1U << (N - 1)); v_xor_w++) {
if (__builtin_popcount(v_xor_w) != N / 2) continue;
for (int k = 0; k < g_necklace_count; k++) {
unsigned int v = g_necklace[k].v;
unsigned int w = v ^ v_xor_w;
for (int i = 0; i < N + 1; i++) {
leading_zero_count[i] += g_necklace[k].multiplicity;
w = ((w & 1) << (N - 1)) | (w >> 1);
if (__builtin_popcount(v ^ w) != N / 2) break;
}
}
}
for (int i = 0; i < N + 1; i++) {
printf(" %lld", 2 * leading_zero_count[i]);
}
putchar(''/n'');
return 0;
}
Puede obtener un poco de aceleración explotando la simetría de signos (4x) e iterando solo sobre los vectores que pasan la primera prueba interna del producto (asintóticamente, O (sqrt (n)) x).
import itertools
n = 10
m = n + 1
def innerproduct(A, B):
s = 0
for k in range(n):
s += A[k] * B[k]
return s
leadingzerocounts = [0] * m
for S in itertools.product([-1, 1], repeat=n - 1):
S1 = S + (1,)
S1S1 = S1 * 2
for C in itertools.combinations(range(n - 1), n // 2):
F = list(S1)
for i in C:
F[i] *= -1
leadingzerocounts[0] += 4
for i in range(1, m):
if innerproduct(F, S1S1[i:i + n]):
break
leadingzerocounts[i] += 4
print(leadingzerocounts)
Versión C, para tener una idea de la cantidad de rendimiento que estamos perdiendo con PyPy (16 para PyPy es aproximadamente 18 para C):
#include <stdio.h>
enum {
HALFN = 9,
N = 2 * HALFN
};
int main(void) {
long long lzc[N + 1];
for (int i = 0; i < N + 1; i++) lzc[i] = 0;
unsigned int xor = 1 << (N - 1);
while (xor-- > 0) {
if (__builtin_popcount(xor) != HALFN) continue;
unsigned int s = 1 << (N - 1);
while (s-- > 0) {
lzc[0]++;
unsigned int f = xor ^ s;
for (int i = 1; i < N + 1; i++) {
f = ((f & 1) << (N - 1)) | (f >> 1);
if (__builtin_popcount(f ^ s) != HALFN) break;
lzc[i]++;
}
}
}
for (int i = 0; i < N + 1; i++) printf(" %lld", 4 * lzc[i]);
putchar(''/n'');
return 0;
}
Este algoritmo es vergonzosamente paralelo porque simplemente se está acumulando sobre todos los valores de xor . Con la versión C, un cálculo del fondo del sobre sugiere que unas pocas miles de horas de tiempo de CPU serían suficientes para calcular n = 26 , lo que equivale a unos doscientos dólares a las tasas actuales en EC2. Sin duda, hay algunas optimizaciones por realizar (p. Ej., Vectorización), pero para un caso único como este, no estoy seguro de cuánto más valga la pena el esfuerzo de los programadores.
He tenido la oportunidad de transferir esto a los arreglos NumPy y he tomado prestado de esta pregunta: los productos de itertools se aceleran
Esto es lo que tengo (puede haber más aceleraciones aquí):
def find_leading_zeros(n):
if n % 2:
return numpy.zeros(n)
m = n+1
leading_zero_counts = numpy.zeros(m)
product_list = [-1, 1]
repeat = n
s = (numpy.array(product_list)[numpy.rollaxis(numpy.indices((len(product_list),) * repeat),
0, repeat + 1).reshape(-1, repeat)]).astype(''int8'')
i = 0
size = s.shape[0] / 2
products = numpy.zeros((size, size), dtype=bool)
while i < m:
products += (numpy.tensordot(s[0:size, 0:size],
numpy.roll(s, i, axis=1)[0:size, 0:size],
axes=(-1,-1))).astype(''bool'')
leading_zero_counts[i] = (products.size - numpy.sum(products)) * 4
i += 1
return leading_zero_counts
Corriendo para n = 14 obtengo:
>>> find_leading_zeros(14)
array([ 56229888., 23557248., 9903104., 4160640., 1758240.,
755392., 344800., 172320., 101312., 75776.,
65696., 61216., 59200., 59200., 59200.])
Así que todo se ve bien. En cuanto a la velocidad:
>>> timeit.timeit("find_leading_zeros_old(10)", number=10)
28.775046825408936
>>> timeit.timeit("find_leading_zeros(10)", number=10)
2.236745834350586
Mira lo que piensas.
EDITAR:
La versión original usó 2074MB de memoria para N = 14, así que numpy.roll matriz concatenada y usé numpy.roll en numpy.roll lugar. También al cambiar los tipos de datos para usar una matriz booleana, la memoria baja a 277MB para n = 14.
En el tiempo la edición es un poco más rápida de nuevo:
>>> timeit.timeit("find_leading_zeros(10)", number=10)
1.3816070556640625
EDIT2:
Ok, así que agregando la simetría como lo señaló David, reduzco esto nuevamente. Ahora usa 213MB. El tiempo de comparación en comparación con ediciones anteriores:
>>> timeit.timeit("find_leading_zeros(10)", number=10)
0.35357093811035156
Ahora puedo hacer el caso n = 14 en 14 segundos en mi Mac book, lo cual no es malo para "python puro", creo.
Intenté acelerar esto y fallé mal :( Pero estoy enviando el código, de alguna manera es más rápido, pero no lo suficientemente rápido para valores como n=24 .
Mis suposiciones
Sus listas constan de valores, por lo que decidí usar números en lugar de listas: cada bit representa uno de los valores posibles: si se establece un bit, entonces esto significa 1 , si se pone a cero significa -1 . El único resultado posible de la multiplicación {-1, 1} es 1 o -1 , así que usé XOR bits en lugar de la multiplicación. También noté que hay una simetría, así que solo necesitas revisar el subconjunto (una cuarta parte) de las listas posibles y multiplicar el resultado por 4 (David explicó esto en su respuesta).
Finalmente puse los resultados de las posibles operaciones en tablas para eliminar la necesidad de cálculos. Se necesita mucha memoria, pero ¿a quién le importa (para n=24 fue alrededor de 150 MB)?
Y luego, @David Eisenstat respondió la pregunta :) Entonces, tomé su código y lo modifiqué a base de bits. Es aproximadamente 2-3 veces más rápido (para n=16 tomó ~ 30 segundos, comparado con ~ 90 de la solución de David), pero creo que todavía no es suficiente para obtener resultados para n=26 o menos.
import itertools
n = 16
m = n + 1
mask = (2 ** n) - 1
# Create table of sum results (replaces innerproduct())
tab = []
for a in range(2 ** n):
s = 0
for k in range(n):
s += -1 if a & 1 else 1
a >>= 1
tab.append(s)
# Create combination bit masks for combinations
comb = []
for C in itertools.combinations(range(n - 1), n // 2):
xor = 0
for i in C:
xor |= (1 << i)
comb.append(xor)
leadingzerocounts = [0] * m
for S in xrange(2 ** (n-1)):
S1 = S + (1 << (n-1))
S1S1 = S1 + (S1 << n)
for xor in comb:
F = S1 ^ xor
leadingzerocounts[0] += 4
for i in range(1, m):
if tab[F ^ ((S1S1 >> i) & mask)]:
break
leadingzerocounts[i] += 4
print(leadingzerocounts)
Conclusiones
Pensé que había inventado algo brillante y esperaba que todo este lío con los bits proporcionara un gran impulso de velocidad, pero el impulso fue decepcionantemente pequeño :(
Creo que la razón es la forma en que Python utiliza los operadores: llama a la función para cada operación aritmética (o lógica), incluso si se pudiera realizar mediante un comando de ensamblador único (esperaba que pypy pudiera simplificar las operaciones a ese nivel, pero no lo hizo). ''t). Por lo tanto, probablemente si se usara C (o ASM) con esta solución operativa de bits, funcionaría muy bien (quizás podría llegar a n=24 ).
Una simple aceleración de un factor de n es cambiar este código:
def innerproduct(A, B):
assert (len(A) == len(B))
for j in xrange(len(A)):
s = 0
for k in xrange(0,n):
s+=A[k]*B[k]
return s
a
def innerproduct(A, B):
assert (len(A) == len(B))
s = 0
for k in xrange(0,n):
s+=A[k]*B[k]
return s
(No sé por qué tiene el bucle sobre j, pero solo hace el mismo cálculo cada vez, por lo que es innecesario).