Accelerazione della moltiplicazione e dell'esponenziazione di matrice vettoriale in Python, possibilmente chiamando C/C++

Question

Attualmente sto lavorando a un progetto di apprendimento automatico in cui - data una matrice di dati Z e un vettore rho - Devo calcolare il valore e la pendenza di il logistic loss function allo rho. Il calcolo comporta operazioni di moltiplicazione di matrice vettore-matrice e operazioni log/exp, con un trucco per evitare l'overflow numerico (descritto in questo previous post).

Attualmente sto facendo questo in Python usando NumPy come mostrato di seguito (come riferimento, questo codice viene eseguito in 0.2 s). Anche se questo funziona bene, vorrei accelerarlo poiché chiamo la funzione più volte nel mio codice (e rappresenta oltre il 90% del calcolo coinvolto nel mio progetto).

Sto cercando un modo per migliorare il runtime di questo codice senza parallelizzazione (vale a dire solo 1 CPU). Sono felice di utilizzare qualsiasi pacchetto disponibile pubblicamente in Python, o di chiamare C o C++ (poiché ho sentito che questo migliora i tempi di esecuzione di un ordine di grandezza). Anche la pre-elaborazione della matrice di dati Z andrebbe bene. Alcune cose che potrebbero essere sfruttati per meglio calcolo sono che il vettore rho è generalmente scarsa (con circa il 50% immissioni = 0) e ci sono solitamente molto più righe di colonne (nella maggior parte dei casi 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

Le librerie della famiglia BLAS sono già altamente ottimizzate per le migliori prestazioni. Quindi nessuno sforzo per collegarsi ad un codice C/C++ può darti dei benefici. Potresti comunque provare varie implementazioni BLAS, dal momento che ce ne sono alcune in giro, incluse alcune ottimizzate per alcune CPU.

L'altra cosa che mi viene in mente è quello di utilizzare una libreria come theano (o di tensorflow Google) che è in grado di rappresentare l'intero grafico computazionale (tutte le operazioni nella vostra funzione di cui sopra) e applicare le ottimizzazioni globali ad esso. Può quindi generare codice CPU da quel grafico tramite C++ (e capovolgendo un semplice switch anche codice GPU). Può anche calcolare automaticamente derivati ​​simbolici per te. Ho usato theano per problemi di machine learning ed è davvero una grande libreria, anche se non è la più facile da imparare.

(sto postando questo come una risposta perché è troppo lungo per un commento)

Edit:

realtà ho avuto un andare a questo a Teano, ma il risultato è in realtà circa 2x più lento sulla CPU, vedi sotto perché. Vi posto qui in ogni caso, forse è un punto di partenza per qualcun altro a fare qualcosa di meglio: (questo è solo il codice parziale, con il codice dal post originale)

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 è abbastanza ottimizzato. Il meglio che puoi fare è provare altre librerie con dati della stessa dimensione inizializzati in modo casuale (non inizializzati a 0) e fai il tuo benchmark.

Se si vuole provare, si può ovviamente provare BLAS. Dovresti provare anche a eigen, personalmente l'ho trovato più veloce su una delle mie applicazioni.