Acelerando a multiplicação e exponenciação de vetores matriciais em Python, possivelmente chamando C / C ++

Atualmente, estou trabalhando em um projeto de aprendizado de máquina em que - dada uma matriz de dadosZ e um vetorrho - Eu tenho que calcular o valor e a inclinação dofunção de perda logística àsrho. O cálculo envolve operações básicas de multiplicação de vetores matriciais e operações log / exp, com um truque para evitar o excesso numérico (descrito nestepostagem anterior)

Atualmente, estou fazendo isso no Python usando o NumPy, como mostrado abaixo (como referência, esse código é executado em 0.2s). Embora isso funcione bem, eu gostaria de acelerá-lo, já que chamo a função várias vezes no meu código (e representa mais de 90% da computação envolvida no meu projeto).

Estou procurando alguma maneira de melhorar o tempo de execução desse código sem paralelização (ou seja, apenas 1 CPU). Fico feliz em usar qualquer pacote disponível publicamente em Python ou em chamar C ou C ++ (desde que ouvi dizer que isso melhora os tempos de execução por uma ordem de magnitude). Pré-processando a Matriz de DadosZ também seria bom. Algumas coisas que poderiam ser exploradas para uma melhor computação são que o vetorrho geralmente é escasso (com cerca de 50% das entradas = 0) e geralmente existemlonge mais linhas do que colunas (na maioria dos casosn_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)

questionAnswers(2)

yourAnswerToTheQuestion