Przyspieszenie matrycę wektor mnożenie i potęgowanie w Pythonie, ewentualnie poprzez wywołanie C/C++

Question

Obecnie pracuję nad projektem uczenia maszynowego gdzie - danej macierzy danych Z oraz wektor rho - Muszę obliczyć wartość i nachylenie logistic loss function pod adresem rho. Obliczenia obejmują podstawowe operacje mnożenia macierzy-wektora i log/exp, z trikiem, aby uniknąć przepełnienia liczbowego (opisanego w tym previous post).

Obecnie robię to w Pythonie za pomocą NumPy, jak pokazano poniżej (jako odniesienie, ten kod działa w 0.2s). Chociaż to działa dobrze, chciałbym go przyspieszyć, ponieważ wiele razy nazywam tę funkcję w moim kodzie (i reprezentuje ona ponad 90% obliczeń związanych z moim projektem).

Szukam sposobu na poprawę środowiska wykonawczego tego kodu bez równoległego (tj. Tylko 1 procesora). Cieszę się, używając dowolnego publicznie dostępnego pakietu w Pythonie, lub wywołując C lub C++ (ponieważ słyszałem, że poprawia to runtimes o rząd wielkości). Wstępne przetwarzanie matrycy danych Z również byłoby OK. Niektóre rzeczy, które mogłyby być wykorzystane do lepszego obliczeń są takie, że wektor rho jest zwykle rzadki (o około 50% zgłoszeń = 0) i są zwykle daleko więcej wierszy niż kolumn (w większości przypadków 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

Biblioteki rodziny BLAS są już wysoko dostrojone, aby uzyskać najlepszą wydajność. Zatem żadne wysiłki w celu połączenia z jakimś kodem C/C++ nie przyniosą żadnych korzyści. Możesz jednak wypróbować różne implementacje BLAS, ponieważ jest ich sporo, w tym niektóre specjalnie przystosowane do niektórych procesorów.

Inną rzeczą, która przychodzi mi do głowy, jest użycie biblioteki takiej jak theano (lub Google tensorflow), która jest w stanie przedstawić cały wykres obliczeniowy (wszystkie operacje w twojej funkcji powyżej) i zastosować do niego globalne optymalizacje. Następnie może wygenerować kod CPU z tego wykresu za pośrednictwem C++ (i przerzucając prosty przełącznik również kod GPU). Może również automatycznie obliczać dla ciebie symboliczne pochodne. Używałem theano do problemów z uczeniem maszynowym i jest to świetna biblioteka do tego, chociaż nie jest najłatwiejsza do nauczenia się.

(jestem delegowania to jako odpowiedź, ponieważ jest zbyt długi dla komentarzu)

Edit:

I rzeczywiście miał iść w tym w Theano, ale wynik jest faktycznie około 2x wolniej na CPU, patrz poniżej, dlaczego. Ja to pisać tutaj i tak, może to stanowić punkt wyjścia dla kogoś do zrobienia czegoś lepiej: (jest to tylko częściowy kod, wraz z kodem z oryginalnego postu)

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 jest dość zoptymalizowane. Najlepsze, co możesz zrobić, to wypróbować inne biblioteki z danymi o tym samym rozmiarze zainicjowanym losowo (nie zainicjowanym na 0) i zrobić własny test porównawczy.

Jeśli chcesz spróbować, możesz oczywiście wypróbować BLAS. Powinieneś także spróbować uzyskać eigen, osobiście znalazłem to szybciej w jednej z moich aplikacji.