为何我的神经网络中Numpy版本比纯Python版本运行更慢？

我正在学习神经网络，分别用纯Python和Numpy实现了线性回归模型，原本预期Numpy版本会更快，但实际运行下来它慢了很多。运行环境是M1 MacBook Pro，Python 3.9.13，Numpy 1.24.1。

核心原因分析 #

1. 数据集规模太小，Numpy的开销抵消了优势 #

你的数据集只有30个样本，Numpy的向量化操作依赖C级别的底层计算加速，但小数据量下，Python列表的循环、求和等操作的开销反而比Numpy的数组内存管理、维度校验等额外开销更小。纯Python的列表推导和sum操作在小样本场景下的效率反而更高。

2. M1芯片的Numpy适配问题 #

早期版本的Numpy在Apple Silicon架构上可能未完全优化，比如未充分利用ARM的NEON指令集，或者使用了通用编译版本而非针对M1定制的版本。你使用的Numpy 1.24.1虽不算老旧，但仍可能存在部分适配短板，导致向量化操作未发挥出应有的性能。

3. 计算逻辑的细节差异 #

• 纯Python的sigmoid函数用了近似值e=2.718，而Numpy版本使用了精确的np.exp(z)，高精度计算带来了额外开销。
• 每次迭代中，Numpy的数组操作都要处理类型、维度等信息，100万次循环累计下来，这部分开销会被显著放大。

优化建议 #

• 增大数据集规模：当样本量达到数千甚至数万时，Numpy的向量化优势会完全体现，纯Python的循环会因解释器开销急剧变慢，而Numpy的C级计算会大幅领先。
• 升级或优化Numpy版本：尝试升级到Numpy 1.26+的新版本，或通过Homebrew安装针对Apple Silicon优化的Numpy，充分利用M1的硬件特性。
• 减少循环内的重复计算：在train方法的循环中，可合并重复的梯度、损失计算逻辑，避免每次迭代都创建新数组，复用中间结果。
• 优化sigmoid实现：改用1/(1+np.exp(-z)) - 0.5的写法，减少不必要的除法操作；或者使用scipy.special.expit实现更高效的sigmoid计算。

纯Python实现代码 #

class PureNeuralNetwork():

    def sigmoid(self, z):
        e = 2.718
        z = min(100, z)
        return (e ** z) / (1 + (e ** z)) - 0.5

    def predict(self, w, X, b):
        return w * X + b

    def loss(self, w, b, X, Y): # the average squared loss
        losses = [(self.predict(w, X[i], b) - Y[i]) ** 2 for i in range(len(X))]  
        return sum(losses) / len(X)

    def gradient_w(self, w, b, X, Y):
        inner = [2 * (w * X[i] + b - Y[i]) * X[i] for i in range(len(X))]
        grad = sum(inner) / len(X)
        return self.sigmoid(grad)

    def gradient_b(self, w, b, X, Y):
        inner = [2 * (w * X[i] + b - Y[i]) for i in range(len(X))]
        grad = sum(inner) / len(X)
        return self.sigmoid(grad)

    def accuracy(self, w, b, X, Y):
        accuracies = [1 - (abs(self.predict(w, X[i], b) - Y[i]) / Y[i]) for i in range(len(X))]
        return sum(accuracies) / len(accuracies) # average mean

    def report(self, i, w, b, X, Y, wgradient, bgradient):
        print("\nITERATION ", i)
        print("W        = ", w)
        print("B        = ", b)
        print("ERROR    = ", self.loss(w, b, X, Y))
        print("GRADIENT_W = ", wgradient)
        print("GRADIENT_B = ", bgradient)
        print("ACCURACY = ", self.accuracy(w, b, X, Y))

    def train(self, w, b, X, Y, iterations, lr):
        for i in range(iterations):
            wgradient = self.gradient_w(w, b, X, Y) * lr
            bgradient = self.gradient_b(w, b, X, Y) * lr

            w -= wgradient
            b -= bgradient

            if (i % (iterations / 10 ) == 0):
                self.report(i, w, b, X, Y, wgradient, bgradient)
    
        return w, b

    
# COLLECT DATA

f = open("pizza_comma.txt", "r")
all_nums = []
temp = []
for i in f.read():
    if i.isnumeric():
        temp.append(i)

    if i == "," or i == "\n":
        temp_string = ''.join(temp)
        all_nums.append(int(temp_string))
        temp = []
X = []
Y = []
for i in range(len(all_nums)):
    if i % 2 == 0:
        X.append(all_nums[i])
    else:
        Y.append(all_nums[i])

# CODE

NN = PureNeuralNetwork()

iterations = 1000000
w = 1
b = 1
lr = 0.01

w, b = NN.train(w, b, X, Y, iterations, lr)

print("\nFinal accuracy = ", NN.accuracy(w, b, X, Y))

纯Python版使用的数据（pizza_comma.txt） #

13,33
2,16
14,32
23,51
13,27
1,16
18,34
10,17
26,29
3,15
3,15
21,32
7,22
22,37
2,13
27,44
6,16
10,21
18,37
15,30
9,26
26,34
8,23
15,39
10,27
21,37
5,17
6,18
13,25
13,23

Numpy实现代码 #

import numpy as np

class NpNeuralNetwork():

    def sigmoid(self, z):
        z = np.minimum(100, z)
        return np.exp(z) / (1 + np.exp(z)) - 0.5

    def predict(self, w, X, b):
        return w * X + b

    def loss(self, w, b, X, Y):
        squared_error = np.power((self.predict(w, X, b) - Y), 2)
        return np.mean(squared_error)

    def gradient_w(self, w, b, X, Y):
        inner = 2 * (w * X + b - Y) * X
        mean = np.mean(inner)
        return self.sigmoid(mean)

    def gradient_b(self, w, b, X, Y):
        inner = 2 * (w * X + b - Y)
        mean = np.mean(inner)
        return self.sigmoid(mean)

    def accuracy(self, w, b, X, Y):
        accuracies = 1 - (np.absolute(self.predict(w, X, b) - Y) / Y)
        return np.mean(accuracies)

    def report(self, i, w, b, X, Y, wgradient, bgradient):
        print("\nITERATION ", i)
        print("W        = ", w)
        print("B        = ", b)
        print("ERROR    = ", self.loss(w, b, X, Y))
        print("GRADIENT_W = ", wgradient)
        print("GRADIENT_B = ", bgradient)
        print("ACCURACY = ", self.accuracy(w, b, X, Y))

    def train(self, X, Y, iterations, lr, w=1, b=1):
        for i in range(iterations):
            wgradient = self.gradient_w(w, b, X, Y) * lr
            bgradient = self.gradient_b(w, b, X, Y) * lr

            w -= wgradient
            b -= bgradient

            if (i % (iterations / 10 ) == 0):
                self.report(i, w, b, X, Y, wgradient, bgradient)

        
        return w, b



data = np.loadtxt("pizza_space.txt").T
X = data[0]
Y = data[1]

NN = NpNeuralNetwork()

w, b = NN.train(X, Y, 1000000, 0.01)

print("\nFinal accuracy = ", NN.accuracy(w, b, X, Y))

Numpy版使用的数据（pizza_space.txt） #

13  33
2  16
14  32
23  51
13  27
1  16
18  34
10  17
26  29
3  15
3  15
21  32
7  22
22  37
2  13
27  44
6  16
10  21
18  37
15  30
9  26
26  34
8  23
15  39
10  27
21  37
5  17
6  18
13  25
13  23

内容的提问来源于stack exchange，提问作者shrivelbottom