python - neural - Diferencia en el rendimiento entre numpy y matlab
que es mejor matlab o c++ (3)
Estoy computando el algoritmo backpropagation para un autoencoder escaso. Lo he implementado en python usando numpy y en matlab . El código es casi el mismo, pero el rendimiento es muy diferente. El tiempo que tarda el matlab en completar la tarea es 0.252454 segundos mientras que numpy 0.973672151566, es decir, casi cuatro veces más. Llamaré a este código varias veces más tarde en un problema de minimización, por lo que esta diferencia lleva a varios minutos de retraso entre las implementaciones. ¿Es esto un comportamiento normal? ¿Cómo podría mejorar el rendimiento en numpy?
Implementación Numpy:
Sparse.rho es un parámetro de ajuste, sparse.nodes es el número de nodos en la capa oculta (25), sparse.input (64) el número de nodos en la capa de entrada, theta1 y theta2 son las matrices de peso para el primero y segunda capa respectivamente con dimensiones 25x64 y 64x25, m es igual a 10000, rhoest tiene una dimensión de (25), x tiene una dimensión de 10000x64, a3 10000x64 y a2 10000x25.
UPDATE : He introducido cambios en el código siguiendo algunas de las ideas de las respuestas. El rendimiento ahora es numpy: 0.65 vs matlab: 0.25.
partial_j1 = np.zeros(sparse.theta1.shape)
partial_j2 = np.zeros(sparse.theta2.shape)
partial_b1 = np.zeros(sparse.b1.shape)
partial_b2 = np.zeros(sparse.b2.shape)
t = time.time()
delta3t = (-(x-a3)*a3*(1-a3)).T
for i in range(m):
delta3 = delta3t[:,i:(i+1)]
sum1 = np.dot(sparse.theta2.T,delta3)
delta2 = ( sum1 + sum2 ) * a2[i:(i+1),:].T* (1 - a2[i:(i+1),:].T)
partial_j1 += np.dot(delta2, a1[i:(i+1),:])
partial_j2 += np.dot(delta3, a2[i:(i+1),:])
partial_b1 += delta2
partial_b2 += delta3
print "Backprop time:", time.time() -t
Implementación de Matlab:
tic
for i = 1:m
delta3 = -(data(i,:)-a3(i,:)).*a3(i,:).*(1 - a3(i,:));
delta3 = delta3.'';
sum1 = W2.''*delta3;
sum2 = beta*(-sparsityParam./rhoest + (1 - sparsityParam) ./ (1.0 - rhoest) );
delta2 = ( sum1 + sum2 ) .* a2(i,:).'' .* (1 - a2(i,:).'');
W1grad = W1grad + delta2* a1(i,:);
W2grad = W2grad + delta3* a2(i,:);
b1grad = b1grad + delta2;
b2grad = b2grad + delta3;
end
toc
Comenzaría con operaciones en el lugar para evitar asignar nuevas matrices cada vez:
partial_j1 += np.dot(delta2, a1[i,:].reshape(1,a1.shape[1]))
partial_j2 += np.dot(delta3, a2[i,:].reshape(1,a2.shape[1]))
partial_b1 += delta2
partial_b2 += delta3
Puedes reemplazar esta expresión:
a1[i,:].reshape(1,a1.shape[1])
de una manera más simple y rápida (gracias a Bi Rico ):
a1[i:i+1]
Además, esta línea:
sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest))
parece ser el mismo en cada ciclo, no necesita recalcularlo.
Y, una optimización probablemente menor, puede reemplazar todas las ocurrencias de x[i,:] con x[i] .
Finalmente, si puede permitirse asignar los m veces más memoria, puede seguir la sugestión unutbu y vectorizar el ciclo:
for m in range(m):
delta3 = -(x[i]-a3[i])*a3[i]* (1 - a3[i])
con:
delta3 = -(x-a3)*a3*(1-a3)
Y siempre puedes usar Numba y ganar velocidad significativamente sin vectorizar (y sin usar más memoria).
La diferencia en el rendimiento entre numpy y matlab siempre me ha frustrado. A menudo, al final se reducen a las bibliotecas subyacentes de lapack. Por lo que sé, matlab usa el atlas lapack completo como predeterminado, mientras que Numpy usa una luz lapack. Matlab reconoce que las personas no se preocupan por el espacio y el volumen, mientras que Numpy cree que las personas sí. Pregunta similar con una buena respuesta.
Sería un error decir "Matlab siempre es más rápido que NumPy" o viceversa. A menudo, su rendimiento es comparable. Al usar NumPy, para obtener un buen rendimiento, debe tener en cuenta que la velocidad de NumPy proviene de llamar a funciones subyacentes escritas en C / C ++ / Fortran. Funciona bien cuando aplica esas funciones a arreglos completos. En general, obtiene un peor rendimiento cuando llama a esa función NumPy en matrices o escalares más pequeños en un bucle de Python.
¿Qué pasa con un ciclo de Python que preguntas? Cada iteración a través del ciclo de Python es una llamada al next método. Cada uso de indexación de [] es una llamada a un método __getitem__ . Cada += es una llamada a __iadd__ . Cada búsqueda de atributos punteados (como en np.dot similar) implica llamadas a funciones. Esas llamadas a funciones se suman a un obstáculo significativo a la velocidad. Estos ganchos le dan a Python un poder expresivo: la indexación de cadenas significa algo diferente a la indexación de dicts, por ejemplo. Misma sintaxis, diferentes significados. La magia se logra dando a los objetos diferentes métodos __getitem__ .
Pero ese poder expresivo tiene un costo en velocidad. Entonces, cuando no necesita toda esa expresividad dinámica, para obtener un mejor rendimiento, trate de limitarse a las llamadas a la función NumPy en arreglos completos.
Por lo tanto, elimine el for-loop; use ecuaciones "vectorizadas" cuando sea posible. Por ejemplo, en lugar de
for i in range(m):
delta3 = -(x[i,:]-a3[i,:])*a3[i,:]* (1 - a3[i,:])
puedes calcular delta3 para cada uno de una vez:
delta3 = -(x-a3)*a3*(1-a3)
Mientras que en el for-loop delta3 es un vector, usar la ecuación vectorizada delta3 es una matriz.
Algunos de los cálculos en for-loop no dependen de i y, por lo tanto, deben levantarse fuera del ciclo. Por ejemplo, sum2 parece una constante:
sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho) / (1.0 - rhoest) )
Aquí hay un ejemplo ejecutable con una implementación alternativa ( alt ) de su código ( orig ).
Mi punto de referencia de tiempo muestra una mejora de velocidad de 6.8 veces :
In [52]: %timeit orig()
1 loops, best of 3: 495 ms per loop
In [53]: %timeit alt()
10 loops, best of 3: 72.6 ms per loop
import numpy as np
class Bunch(object):
""" http://code.activestate.com/recipes/52308 """
def __init__(self, **kwds):
self.__dict__.update(kwds)
m, n, p = 10 ** 4, 64, 25
sparse = Bunch(
theta1=np.random.random((p, n)),
theta2=np.random.random((n, p)),
b1=np.random.random((p, 1)),
b2=np.random.random((n, 1)),
)
x = np.random.random((m, n))
a3 = np.random.random((m, n))
a2 = np.random.random((m, p))
a1 = np.random.random((m, n))
sum2 = np.random.random((p, ))
sum2 = sum2[:, np.newaxis]
def orig():
partial_j1 = np.zeros(sparse.theta1.shape)
partial_j2 = np.zeros(sparse.theta2.shape)
partial_b1 = np.zeros(sparse.b1.shape)
partial_b2 = np.zeros(sparse.b2.shape)
delta3t = (-(x - a3) * a3 * (1 - a3)).T
for i in range(m):
delta3 = delta3t[:, i:(i + 1)]
sum1 = np.dot(sparse.theta2.T, delta3)
delta2 = (sum1 + sum2) * a2[i:(i + 1), :].T * (1 - a2[i:(i + 1), :].T)
partial_j1 += np.dot(delta2, a1[i:(i + 1), :])
partial_j2 += np.dot(delta3, a2[i:(i + 1), :])
partial_b1 += delta2
partial_b2 += delta3
# delta3: (64, 1)
# sum1: (25, 1)
# delta2: (25, 1)
# a1[i:(i+1),:]: (1, 64)
# partial_j1: (25, 64)
# partial_j2: (64, 25)
# partial_b1: (25, 1)
# partial_b2: (64, 1)
# a2[i:(i+1),:]: (1, 25)
return partial_j1, partial_j2, partial_b1, partial_b2
def alt():
delta3 = (-(x - a3) * a3 * (1 - a3)).T
sum1 = np.dot(sparse.theta2.T, delta3)
delta2 = (sum1 + sum2) * a2.T * (1 - a2.T)
# delta3: (64, 10000)
# sum1: (25, 10000)
# delta2: (25, 10000)
# a1: (10000, 64)
# a2: (10000, 25)
partial_j1 = np.dot(delta2, a1)
partial_j2 = np.dot(delta3, a2)
partial_b1 = delta2.sum(axis=1)
partial_b2 = delta3.sum(axis=1)
return partial_j1, partial_j2, partial_b1, partial_b2
answer = orig()
result = alt()
for a, r in zip(answer, result):
try:
assert np.allclose(np.squeeze(a), r)
except AssertionError:
print(a.shape)
print(r.shape)
raise
Consejo: Note que dejé en los comentarios la forma de todas las matrices intermedias. Conocer la forma de las matrices me ayudó a entender lo que estaba haciendo su código. La forma de las matrices puede guiarlo hacia las funciones de NumPy correctas que debe usar. O al menos, prestar atención a las formas puede ayudarlo a saber si una operación es sensata. Por ejemplo, cuando calcula
np.dot(A, B)
y A.shape = (n, m) y B.shape = (m, p) , entonces np.dot(A, B) será una matriz de forma (n, p) .
Puede ayudar a construir las matrices en orden C_CONTIGUA (al menos, si usa np.dot ). Puede haber tanto como 3 veces la velocidad al hacerlo:
A continuación, x es lo mismo que xf excepto que x es C_CONTIGUA y xf es F_CONTIGUA - y la misma relación para y y yf .
import numpy as np
m, n, p = 10 ** 4, 64, 25
x = np.random.random((n, m))
xf = np.asarray(x, order=''F'')
y = np.random.random((m, n))
yf = np.asarray(y, order=''F'')
assert np.allclose(x, xf)
assert np.allclose(y, yf)
assert np.allclose(np.dot(x, y), np.dot(xf, y))
assert np.allclose(np.dot(x, y), np.dot(xf, yf))
%timeit puntos de referencia %timeit muestran la diferencia de velocidad:
In [50]: %timeit np.dot(x, y)
100 loops, best of 3: 12.9 ms per loop
In [51]: %timeit np.dot(xf, y)
10 loops, best of 3: 27.7 ms per loop
In [56]: %timeit np.dot(x, yf)
10 loops, best of 3: 21.8 ms per loop
In [53]: %timeit np.dot(xf, yf)
10 loops, best of 3: 33.3 ms per loop
En cuanto a la evaluación comparativa en Python:
Puede ser engañoso utilizar la diferencia en pares de llamadas time.time() para comparar la velocidad del código en Python. Debe repetir la medición muchas veces. Es mejor desactivar el recolector de basura automático. También es importante medir grandes intervalos de tiempo (por ejemplo, al menos 10 segundos de repeticiones) para evitar errores debido a una resolución deficiente en el temporizador de reloj y para reducir la importancia de la time.time llamadas time.time . En lugar de escribir todo ese código usted mismo, Python le proporciona el módulo de tiempo . Básicamente, estoy usando eso para sincronizar las partes del código, excepto que lo llamo a través de un terminal IPython para mayor comodidad.
No estoy seguro de si esto está afectando sus puntos de referencia, pero tenga en cuenta que podría marcar la diferencia. En la pregunta a la que me time.time , según time.time dos fragmentos de código diferían en un factor de 1.7x, mientras que los benchmarks que usaban timeit mostraban que los fragmentos de código se ejecutaban en cantidades de tiempo esencialmente idénticas.