在 Python 中加速矩阵向量乘法和求幂，可能通过调用 C/C++答案

【问题标题】：Speeding up matrix-vector multiplication and exponentiation in Python, possibly by calling C/C++在 Python 中加速矩阵向量乘法和求幂，可能通过调用 C/C++
【发布时间】：2016-05-11 13:44:43
【问题描述】：

我目前正在从事一个机器学习项目，其中 - 给定一个数据矩阵 Z 和一个向量 rho - 我必须在 rho 处计算 logistic loss function 的值和斜率。计算涉及基本的矩阵向量乘法和 log/exp 运算，并带有避免数值溢出的技巧（在此 previous post 中进行了描述）。

我目前正在使用 NumPy 在 Python 中执行此操作，如下所示（作为参考，此代码在 0.2 秒内运行）。虽然这很好用，但我想加快速度，因为我在代码中多次调用该函数（它代表了我项目中超过 90% 的计算）。

我正在寻找任何方法来改进此代码的运行时间，而无需并行化（即只有 1 个 CPU）。我很高兴在 Python 中使用任何公开可用的包，或者调用 C 或 C++（因为我听说这将运行时间提高了一个数量级）。预处理数据矩阵Z 也可以。可以利用一些东西来进行更好的计算是向量rho 通常是稀疏的（大约 50% 的条目 = 0）并且行通常比列多远（在大多数情况下 @ 987654328@)

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)

【问题讨论】：

为什么要排除超过 1 个 CPU？尽管 Python VM 本质上是单线程的，但在将数据复制到对线程更友好的数据结构之后，您可以从 C 扩展中调用 POSIX 线程。不使用多个 CPU 可能还有其他原因，但如果您转义到 C，则不受该限制的限制。
@rts 好问题。在这种情况下，我需要将其限制为 1 个 CPU，因为调用 compute_logistic_loss_function 的代码实际上是并行化的......所以调用函数时只有 1 个 CPU 可用。
对于大的n，运行时似乎由loss_slope = Z * (phi_slope - 1.0) 主导，它广播到与Z 相同的大小。由于您对行取平均值，因此您可以使用 Z.T.dot(phi_slope).T / Z.shape[0] 将其重写为点积，这在我的机器上提供了大约 4 倍的加速。
最昂贵的操作是矩阵乘积，因此您应该确保您的 numpy 版本与快速 BLAS 库（例如 OpenBLAS 或 MKL）链接。如果我是你，我会使用多线程 BLAS 来并行化矩阵乘积，并放弃调用 Python 代码的并行化。
*对不起，我在上面的评论中指的是Z.T.dot(phi_slope - 1.0).T / Z.shape[0]

标签： python c numpy machine-learning logistic-regression

【解决方案1】：

BLAS 系列的库已针对最佳性能进行了高度调整。因此，不努力链接到某些 C/C++ 代码可能会给您带来任何好处。但是，您可以尝试各种 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

【讨论】：

【解决方案2】：

Numpy 非常优化。您可以做的最好的事情是尝试其他库，将相同大小的数据初始化为随机（未初始化为 0）并进行自己的基准测试。

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

【讨论】：