suma - Acelerar la multiplicación y exponenciación de matrices-vectores en Python, posiblemente llamando a C/C++

multiplicar matrices python numpy

Question

Actualmente estoy trabajando en un proyecto de aprendizaje automático donde, dada una matriz de datos Z y un vector rho , tengo que calcular el valor y la pendiente de la función de pérdida logística en rho . El cálculo involucra operaciones básicas de multiplicación matriz-vector y log / exp, con un truco para evitar el desbordamiento numérico (descrito en esta publicación anterior ).

Actualmente estoy haciendo esto en Python usando NumPy como se muestra a continuación (como referencia, este código se ejecuta en 0.2s). Aunque esto funciona bien, me gustaría acelerarlo ya que llamo a la función varias veces en mi código (y representa más del 90% de la computación involucrada en mi proyecto).

Estoy buscando cualquier forma de mejorar el tiempo de ejecución de este código sin paralelización (es decir, solo 1 CPU). Me complace utilizar cualquier paquete disponible públicamente en Python, o llamar a C o C ++ (ya que he oído que esto mejora los tiempos de ejecución en un orden de magnitud). El preprocesamiento de la matriz de datos Z también estaría bien. Algunas cosas que podrían aprovecharse para un mejor cálculo son que el vector rho suele ser escaso (con alrededor del 50% de entradas = 0) y generalmente hay muchas más filas que columnas (en la mayoría de los casos n_cols <= 100 )

import time
import numpy as np
np.__config__.show() #make sure BLAS/LAPACK is being used
np.random.seed(seed = 0)
#initialize data matrix X and label vector Y
n_rows, n_cols = 1e6, 100
X = np.random.random(size=(n_rows, n_cols))
Y = np.random.randint(low=0, high=2, size=(n_rows, 1))
Y[Y==0] = -1
Z = X*Y # all operations are carried out on Z
def compute_logistic_loss_value_and_slope(rho, Z):
#compute the value and slope of the logistic loss function in a way that is numerically stable
#loss_value: (1 x 1) scalar = 1/n_rows * sum(log( 1 .+ exp(-Z*rho))
#loss_slope: (n_cols x 1) vector = 1/n_rows * sum(-Z*rho ./ (1+exp(-Z*rho))
#see also: https://stackoverflow.com/questions/20085768/
scores = Z.dot(rho)
pos_idx = scores > 0
exp_scores_pos = np.exp(-scores[pos_idx])
exp_scores_neg = np.exp(scores[~pos_idx])
#compute loss value
loss_value = np.empty_like(scores)
loss_value[pos_idx] = np.log(1.0 + exp_scores_pos)
loss_value[~pos_idx] = -scores[~pos_idx] + np.log(1.0 + exp_scores_neg)
loss_value = loss_value.mean()
#compute loss slope
phi_slope = np.empty_like(scores)
phi_slope[pos_idx] = 1.0 / (1.0 + exp_scores_pos)
phi_slope[~pos_idx] = exp_scores_neg / (1.0 + exp_scores_neg)
loss_slope = Z.T.dot(phi_slope - 1.0) / Z.shape[0]
return loss_value, loss_slope
#initialize a vector of integers where more than half of the entries = 0
rho_test = np.random.randint(low=-10, high=10, size=(n_cols, 1))
set_to_zero = np.random.choice(range(0,n_cols), size =(np.floor(n_cols/2), 1), replace=False)
rho_test[set_to_zero] = 0.0
start_time = time.time()
loss_value, loss_slope = compute_logistic_loss_value_and_slope(rho_test, Z)
print "total runtime = %1.5f seconds" % (time.time() - start_time)

Answer 1

Las bibliotecas de la familia BLAS ya están muy optimizadas para obtener el mejor rendimiento. Por lo tanto, es probable que ningún esfuerzo por vincularse a algún código C / C ++ le brinde ningún beneficio. Sin embargo, puedes probar varias implementaciones de BLAS, ya que hay bastantes de ellas, incluidas algunas especialmente adaptadas a algunas CPU.

Otra cosa que me viene a la mente es utilizar una biblioteca como theano (o el flujo de tensor de Google) que sea capaz de representar todo el gráfico computacional (todas las operaciones en su función anterior) y aplicarle optimizaciones globales. A continuación, puede generar código de CPU a partir de ese gráfico a través de C ++ (y volteando un simple interruptor también código de la GPU). También puede calcular automáticamente derivados simbólicos para usted. He utilizado theano para problemas de aprendizaje automático y es una gran biblioteca para eso, aunque no es la más fácil de aprender.

(Estoy publicando esto como respuesta porque es demasiado largo para un comentario)

Editar:

De hecho, tuve un intento en esto en theano, pero el resultado es en realidad 2 veces más lento en la CPU, ver a continuación por qué. Lo publicaré aquí de todos modos, tal vez sea un punto de partida para que alguien más haga algo mejor: (esto es solo un código parcial, completo con el código de la publicación original)

import theano
def make_graph(rho, Z):
scores = theano.tensor.dot(Z, rho)
# this is very inefficient... it calculates everything twice and
# then picks one of them depending on scores being positive or not.
# not sure how to express this in theano in a more efficient way
pos = theano.tensor.log(1 + theano.tensor.exp(-scores))
neg = theano.tensor.log(scores + theano.tensor.exp(scores))
loss_value = theano.tensor.switch(scores > 0, pos, neg)
loss_value = loss_value.mean()
# however computing the derivative is a real joy now:
loss_slope = theano.tensor.grad(loss_value, rho)
return loss_value, loss_slope
sym_rho = theano.tensor.col(''rho'')
sym_Z = theano.tensor.matrix(''Z'')
sym_loss_value, sym_loss_slope = make_graph(sym_rho, sym_Z)
compute_logistic_loss_value_and_slope = theano.function(
inputs=[sym_rho, sym_Z],
outputs=[sym_loss_value, sym_loss_slope]
)
# use function compute_logistic_loss_value_and_slope() as in original code

Answer 2

Numpy está bastante optimizado. Lo mejor que puede hacer es probar otras bibliotecas con datos del mismo tamaño inicializados al azar (no inicializados a 0) y hacer su propio punto de referencia.

Si quieres probar, puedes probar BLAS. También debería intentarlo, personalmente lo encontré más rápido en una de mis aplicaciones.