在Python中加速矩阵向量乘法和指数运算，可能通过调用C/C++来实现

Question

我目前正在进行一个机器学习项目，在这个项目中，给定一个数据矩阵Z和一个向量rho，我必须计算logistic loss function在rho处的值和斜率。计算包括基本的矩阵向量乘法和log/exp运算，以及避免数值溢出的技巧（如本文previous post所述）。在

我目前使用的是Python中的NumPy，如下所示（作为参考，这段代码在0.2s中运行）。虽然这很好用，但是我想加快速度，因为我在代码中多次调用这个函数（它代表了我项目中90%以上的计算量）。在

<强>我正在寻找任何方法来提高这个代码的运行时间而不进行并行化（即只有1个CPU）。< /强>我很高兴使用Python中的任何公开可用的包，或者调用C或C++（因为我听说这会提高运行时间一个数量级）。预处理数据矩阵Z也可以。一些可以用来更好地计算的东西是向量rho通常是稀疏的（大约50%的条目=0），并且通常far行多于列（在大多数情况下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

Numpy非常优化。您所能做的最好的方法是尝试将相同大小的数据初始化为random（不初始化为0）的其他库，并执行自己的基准测试。在

如果你想试试，当然可以试试BLAS。你也应该尝试一下eigen，我个人在我的一个应用程序中发现它更快。在

Answer 2

BLAS家族的库已经进行了高度调整以获得最佳性能。因此，不必链接到某些C/C++代码可能会给你带来任何好处。但是，您可以尝试各种BLAS实现，因为有很多BLAS实现，包括一些专门针对某些cpu进行了调优的实现。在

我想到的另一件事是使用像theano（或Google的tensorflow）这样的库，它能够表示整个计算图（上面函数中的所有操作）并对其应用全局优化。然后它可以通过C++生成该图的CPU代码（并通过简单的GPU代码切换）。它还可以自动计算符号导数。我用过theano来解决机器学习问题，这是一个非常好的库，虽然不是最容易学习的。在

（我将此贴出来作为答复，因为它太长，无法发表评论）

编辑：

实际上，我在theano中尝试过这个，但结果是CPU速度慢了2倍，请看下面的原因。不管怎样，我还是把它贴在这里，也许这是别人做得更好的一个起点：（这只是部分代码，与原帖子中的代码一起完成）

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