如何改进 Python 代码中的循环？

Question

我正在将这段代码从 Matlab 翻译成 Python。代码功能很好，但在 python 中速度很慢。在 Matlab 中，代码运行不到一分钟，在 python 中则需要 30 分钟！！！有python模式经验的人可以帮帮我吗？

# P({ai})
somai = 0
for i in range(1, n):
    somaj = 0
    for j in range(1, n):
        exponencial = math.exp(-((a[i] - a[j]) * (a[i] - a[j])) / dev_a2 - ((b[i] - b[j]) * (b[i] - b[j])) / dev_b2)
        somaj = somaj + exponencial

    somai = somai + somaj

Answer 1

与 MATLAB 一样，我建议您使用 vectorize your code。通过 for 循环进行迭代可能比 MATLAB 和 numpy 的低级实现要慢得多。

---

您的操作(a[i] - a[j])*(a[i] - a[j]) 是所有N 数据点的成对平方欧几里得距离。您可以使用 scipy 的 pdist 和 squareform 函数计算成对距离矩阵 -- pdist、squareform。

然后计算成对距离矩阵A 和B 之间的差，并求和指数衰减。所以你可以得到一个矢量化的代码，比如：

import numpy as np
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform

# Example data
N = 1000
a = np.random.rand(N,1)
b = np.random.rand(N,1)
dev_a2 = np.random.rand()
dev_b2 = np.random.rand()

# `a` is an [N,1] matrix (i.e. column vector)
A = pdist(a, 'sqeuclidean')
# Change to pairwise distance matrix
A = squareform(A)
# Divide all elements by same divisor
A = A / dev_a2

# Then do the same for `b`'s
# `b` is an [N,1] matrix (i.e. column vector)
B = pdist(b, 'sqeuclidean')
B = squareform(B)
B = B / dev_b2

# Calculate exponential decay
expo = np.exp(-(A-B))

# Sum all elements
total = np.sum(expo)

---

这里是迭代方法和这个向量化代码之间的快速时序比较。

N: 1000 | Iter Output: 2729989.851117 | Vect Output: 2732194.924364
Iter time: 6.759 secs | Vect time: 0.031 secs

N: 5000 | Iter Output: 24855530.997400 | Vect Output: 24864471.007726
Iter time: 171.795 secs | Vect time: 0.784 secs

请注意，最终结果并不完全相同。我不确定这是为什么，这可能是我的舍入错误或数学错误，但我会留给你.

Answer 2

TLDR

使用 numpy

为什么是 Numpy？

默认情况下，Python 很慢。 python 的强大功能之一是它可以很好地与 C 配合使用，并且拥有大量的库。可以帮助您听到的是numpy。 Numpy 主要是用 C 语言实现的，如果使用得当，它的速度非常快。诀窍是用这样一种方式来表达代码，让你在 numpy 和 python 之外正确地执行。

代码和结果

import math
import numpy as np

n = 1000
np_a = np.random.rand(n)
a = list(np_a)
np_b = np.random.rand(n)
b = list(np_b)
dev_a2, dev_b2 = (1, 1)

def old():
    somai = 0.0
    for i in range(0, n):
        somaj = 0.0
        for j in range(0, n):
            tmp_1 = -((a[i] - a[j]) * (a[i] - a[j])) / dev_a2
            tmp_2 = -((b[i] - b[j]) * (b[i] - b[j])) / dev_b2
            exponencial = math.exp(tmp_1 + tmp_2)
            somaj += exponencial
        somai += somaj
    return somai


def new():
    tmp_1 = -np.square(np.subtract.outer(np_a, np_a)) / dev_a2
    tmp_2 = -np.square(np.subtract.outer(np_b, np_b)) / dev_a2
    exponential = np.exp(tmp_1 + tmp_2)
    somai = np.sum(exponential)
    return somai

旧 = 1.76 s ± 48.3 ms 每个循环（平均值 ± 标准偏差，7 次运行，每个循环 1 个）
新 = 24.6 ms ± 66.1 µs 每个循环（平均值 ± 标准偏差。7 次运行，每次 10 个循环） 这大约是 70 倍的改进
old 产生 740919.6020840995
new 产生 740919.602084099

说明

您会注意到，为了清楚起见，我将您的代码拆分为 tmp_1 和 tmp_2。

np.random.rand(n)：这将创建一个长度为 n 的数组，其中包含从 0 到 1（不包括 1）的随机浮点数 (documented here)。

np.subtract.outer(a, b)：Numpy 为所有运算符提供了模块，允许你用它们做各种事情。假设你有np_a = [1, 2, 3]，np.subtract.outer(np_a, np_a) 会产生

array([[ 0, -1, -2],
       [ 1,  0, -1],
       [ 2,  1,  0]])

Here's a stackoverflow link 如果您想更深入地了解这一点。 （“外”这个词也来自“外积”，如线性代数）

np.square: 简单地将矩阵中的每个元素平方。

/：在 numpy 中，当您在标量和矩阵之间进行算术运算符时，它会执行适当的操作并将该运算应用于矩阵中的每个元素。

np.exp: 喜欢np.square

np.sum：将每个元素相加并返回一个标量。