优化 Numpy 操作

Posted 2023-03-13

【中文标题】优化 Numpy 操作

【英文标题】：Optimizing Numpy Operations

【发布时间】：2021-11-13 05:00:41

【问题描述】：

我正在尝试使用多项逻辑回归和梯度下降来训练多类分类器。具体来说，该模型将具有经过训练的权重矩阵 w，其形状为 (C, D)，其中 C 是类数，D 是每个输入的特征数。此外，我们将有一个维度为 (C,) 的偏置向量 b。我们有一个 (N, D) 输入矩阵 X，其中 N 是训练输入的数量，以及一个形状为 (N,) 的向量 y，其中 y 中的每个条目是一个从 0 到 C - 1 的数字，表示哪个类输入属于。我写了以下代码：

    for _ in range(max_iterations):
        z = np.apply_along_axis(lambda v: v - max(v), 1, X @ w.T + b)
        probs = np.exp(z)
        denom = np.sum(probs, axis=1)
        for i in range(C):
            for j in range(N):
                if i == y[j]:
                    w[i] -= (step_size / N) * ((probs[j][i] / denom[j]) - 1) * X[j]
                    b[i] -= (step_size / N) * ((probs[j][i] / denom[j]) - 1)
                else:
                    w[i] -= (step_size / N) * (probs[j][i] / denom[j]) * X[j]
                    b[i] -= (step_size / N) * (probs[j][i] / denom[j])

这会产生我想要的正确权重和偏差，但显然它没有利用 numpy 的操作来加快速度。因此，我尝试使用以下代码加快其中一些速度：

    for _ in range(max_iterations):
        z = np.apply_along_axis(lambda v: v - max(v), 1, X @ w.T + b)
        probs = np.exp(z)
        denom = np.sum(probs, axis=1)
        s = np.zeros((N, C))
        for i in range(N):
            s[i] = probs[i] / denom[i]
        for i in range(N):
            s[i][y[i]] += -1
        for c in range(C):
            grad_w = s.T[c] @ X
            w[c] += (step_size / N) * grad_w
            b[c] += (step_size / N) * sum(s.T[c])

我希望这会产生与上一部分相同的结果，同时速度更快......它设法更快，但结果不正确。

所以我有几个问题。首先，为什么我的第二段代码没有产生正确的结果，有什么办法可以解决它？其次，更重要的是，我将如何进一步优化它？这主要是为了让我学习如何利用numpy的向量化操作。

【问题讨论】：

不是说它处理更大的问题，而是np.apply_along_axis 不是一个速度工具。而且我怀疑它可以使用axis 和keepdims 参数（无迭代）替换为np.max。

【参考方案1】：

这可能有助于一些迭代。

从一个小的二维数组开始：

In [251]: probs = np.arange(12).reshape(3,4)
In [252]: denom = np.sum(probs, axis=1)
In [253]: denom
Out[253]: array([ 6, 22, 38])

要将 (3,4) 数组除以 (3,)，我们需要使后面的 (3,1)：

In [254]: probs/denom[:,None]
Out[254]: 
array([[0.        , 0.16666667, 0.33333333, 0.5       ],
       [0.18181818, 0.22727273, 0.27272727, 0.31818182],
       [0.21052632, 0.23684211, 0.26315789, 0.28947368]])

如果没有意义，请阅读并重新阅读 broadcasting 上的 numpy 文档。

获得所需 2d denom 的另一种方法是：

In [255]: denom = np.sum(probs, axis=1, keepdims=True)
In [256]: denom
Out[256]: 
array([[ 6],
       [22],
       [38]])
In [257]: probs/denom
Out[257]: 
array([[0.        , 0.16666667, 0.33333333, 0.5       ],
       [0.18181818, 0.22727273, 0.27272727, 0.31818182],
       [0.21052632, 0.23684211, 0.26315789, 0.28947368]])

与apply_along_axis 一起使用的max 减法同样适用。 apply... 不是速度工具，也不优于简单的迭代。

In [258]: np.max(probs, axis=1, keepdims=True)
Out[258]: 
array([[ 3],
       [ 7],
       [11]])
In [259]: probs - _
Out[259]: 
array([[-3, -2, -1,  0],
       [-3, -2, -1,  0],
       [-3, -2, -1,  0]])