python - 为什么这个脚本在 Python 中这么慢？

Question

我编写了一个脚本，用于在 MATLAB 中训练一维 Kohonen 网络，它非常有用。然后，我尝试将它翻译成 Python 2.7，这是一种我很陌生的语言，而且脚本需要很长时间才能运行。

我将解释我正在做什么，看看这里是否有人可以对此事有所了解。我在矩阵 y 中有一个给定的数据集，我想用它训练不同的 SOM。 SOM 是一维的(一条线)，其神经元数量各不相同。我首先训练一个大小为 N=2 的 SOM，最后训练一个 N=NMax，总共有 NMax-2+1 个 SOM .对于每个 SOM，我想在训练结束后存储权重，然后再进入下一个 SOM。

在 MATLAB 中，对于 NMax = 5 和 iterMax = 50，需要 9.74 秒。在 Python 中，54.04 秒。这种差异是巨大的，而实际数据集、SOM 数量和迭代次数甚至更大，因此 Python 代码需要永远结束。

我当前的代码如下:

import numpy as np
import time
y = np.random.rand(2500,3) # Create random dataset to test
def A(d,s): # Neighborhood function
    return np.exp(-d**2 / (2*s**2))
sigma_0 = float(5) # Initial standard deviation for A
eta_0 = float(1) # Initial learning rate
iterMax = 250 # Maximum number of iterations
NMax = 10 # Maximum number of neurons
w = range(NMax - 1) # Initialize the size of the weight matrix (it will store NMax-2+1 sets of weights, each of varying size depending on the value of N)
#%% KOHONEN 1D
t = time.time() # Start time
for N in np.arange(2,NMax + 1): # Size of the network
    w[N - 2] = np.random.uniform(0,1,(N,np.size(y,axis=1))) - 0.5 # Initial weights
    iterCount = 1; # Iteration counter
    while iterCount < iterMax:
        # Mix the datapoints to choose them in random order
        mixInputs = y[np.random.permutation(np.size(y,axis = 0)),:]
        # Decrease the value of the variance and the learning rate
        sigma = sigma_0 - (sigma_0/(iterMax + 1)) * iterCount
        eta = eta_0 - (eta_0/(iterMax + 1)) * iterCount
        for kk in range(np.size(mixInputs,axis = 0)): # Picking one datapoint at a time
            selectedInput = mixInputs[kk,:]
            # These two lines calculate the weight that is the nearest to the datapoint selected
            aux = np.absolute(np.array(np.kron(np.ones((N,1)),selectedInput)) - np.array(w[N - 2]))
            aux = np.sum(np.abs(aux)**2,axis=-1)
            ii = np.argmin(aux) # The node ii is the winner
            for jj in range(N):
                dist = min(np.absolute(ii-jj) , np.absolute(np.absolute(ii-jj)-N)) # Centering the neighborhood function in the winner
                w[N - 2][jj,:] = w[N - 2][jj,:] + eta * A(dist,sigma) * (selectedInput - w[N - 2][jj,:]) # Updating the weights
        print(N,iterCount)
        iterCount = iterCount + 1 

elapsedTime = time.time() - t

在 MATLAB 中，每次迭代(每次变量 iterCount 增加 1)几乎都是即时的。在 Python 中，每一个都需要很长时间。我不知道为什么它们的表现如此不同，但我想看看是否有可能加速 Python 版本。有什么建议吗？

编辑:按照评论中的要求，这里是代码的更快的 MATLAB 版本。

y = rand(2500,3) % Random dataset
A = @(d,s) exp(-d^2 / (2*s^2));
sigma_0 = 5;
eta_0 = 1;
iterMax = 250;
NMax = 10;
w = cell(NMax - 1,1);

%% KOHONEN 1D
tic();
for N = 2 : NMax 
    w{N - 1} = rand(N,size(y,2)) - 0.5;
    iterCount = 1; 
    while (iterCount < iterMax)
        mixInputs = y(randperm(size(y,1)),:);
        sigma = sigma_0 - (sigma_0/(iterMax + 1)) * iterCount;
        eta = eta_0 - (eta_0/(iterMax + 1)) * iterCount;
        for kk = 1 : size(mixInputs,1)
            input = mixInputs(kk,:);
            % [~,ii] = min(pdist2(input,w{N - 1}));
            aux = abs(repmat(input,N,1) - w{N - 1});
            [~,ii] = min((sum(aux.^2,2))); 
            for jj = 1 : N
                dist = min(abs(ii-jj) , abs(abs(ii-jj)-N));
                w{N - 1}(jj,:) = w{N - 1}(jj,:) + eta * A(dist,sigma) * (input - w{N - 1}(jj,:));
            end
        end
        N % Show N
        iterCount = iterCount + 1 % Show iterCount
    end
end
toc();

Answer 1

使用 profiling module 找出正在调用的函数，以及它们花费了多长时间。

在下面的输出中，各列的含义如下:

ncalls for the number of calls,

tottime for the total time spent in the given function (and excluding time made in calls to sub-functions)

percall is the quotient of tottime divided by ncalls

cumtime is the cumulative time spent in this and all subfunctions (from invocation till exit). This figure is accurate even for recursive functions.

percall is the quotient of cumtime divided by primitive calls

filename:lineno(function) provides the respective data of each function

---

看起来您正在调用 A() 很多很多次...通常使用相同的值。

python2.7 -m cProfile -s tottime ${YOUR_SCRIPT} 的输出

         5481855 function calls (5481734 primitive calls) in 4.986 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.572    1.572    4.986    4.986 x.py:1(<module>)
   214500    0.533    0.000    0.533    0.000 x.py:8(A)
   107251    0.462    0.000    1.986    0.000 shape_base.py:686(kron)
   107251    0.345    0.000    0.456    0.000 numeric.py:1015(outer)
   214502    0.266    0.000    0.563    0.000 {sorted}
...

尝试缓存值:

A_vals = {}

def A(d,s): # Neighborhood function
    t = (d,s)
    if t in A_vals:
        return A_vals[t]
    ret = np.exp(-d**2 / (2*s**2))
    A_vals[t] = ret
    return ret

现在我们看到:

         6206113 function calls (6205992 primitive calls) in 4.986 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.727    1.727    4.986    4.986 x.py:1(<module>)
   121451    0.491    0.000    2.180    0.000 shape_base.py:686(kron)
   121451    0.371    0.000    0.496    0.000 numeric.py:1015(outer)
   242902    0.293    0.000    0.621    0.000 {sorted}
   121451    0.265    0.000    0.265    0.000 {method 'reduce' of 'numpy.ufunc' objects}
...
   242900    0.091    0.000    0.091    0.000 x.py:7(A)
...

此时它变成了一个简单的优化练习!...

您列表中的下一个是 kron() - 享受吧!

---

您还会发现将脚本分解成更小的函数很有帮助(从样式的角度和 分析)。我纯粹出于分析原因做了以下操作 - 所以您最好使用合理的名称，并可能进行更好的分割。

import numpy as np
import time

y = np.random.rand(2500,3) # Create random dataset to test

A_vals = {}
def A(d,s): # Neighborhood function
    t = (d,s)
    if t in A_vals:
        return A_vals[t]
    ret = np.exp(-d**2 / (2*s**2))
    A_vals[t] = ret
    return ret

def a():
    sigma_0 = float(5) # Initial standard deviation for A
    eta_0 = float(1) # Initial learning rate
    iterMax = 250 # Maximum number of iterations
    NMax = 10 # Maximum number of neurons
    w = range(NMax - 1) # Initialize the size of the weight matrix (it will store NMax-2+1 sets of weights, each of varying size depending on the value of N)
    #%% KOHONEN 1D
    t = time.time() # Start time
    for N in np.arange(2,NMax + 1): # Size of the network
        b(w, N, sigma_0, eta_0, iterMax)

    elapsedTime = time.time() - t

def b(w, N, sigma_0, eta_0, iterMax):
    w[N - 2] = np.random.uniform(0,1,(N,np.size(y,axis=1))) - 0.5 # Initial weights
    for iterCount in range(1, iterMax):
        c(N, sigma_0, eta_0, iterMax, iterCount, w)

def c(N, sigma_0, eta_0, iterMax, iterCount, w):
    # Mix the datapoints to choose them in random order
    mixInputs = y[np.random.permutation(np.size(y,axis = 0)),:]
    # Decrease the value of the variance and the learning rate
    sigma = sigma_0 - (sigma_0/(iterMax + 1)) * iterCount
    eta = eta_0 - (eta_0/(iterMax + 1)) * iterCount
    for kk in range(np.size(mixInputs,axis = 0)): # Picking one datapoint at a time
        d(N, w, mixInputs, sigma, eta, kk)
    print(N,iterCount)

def d(N, w, mixInputs, sigma, eta, kk):
    selectedInput = mixInputs[kk,:]
    # These two lines calculate the weight that is the nearest to the datapoint selected
    aux = np.absolute(np.array(np.kron(np.ones((N,1)),selectedInput)) - np.array(w[N - 2]))
    aux = np.sum(np.abs(aux)**2,axis=-1)
    ii = np.argmin(aux) # The node ii is the winner
    for jj in range(N):
        e(N, w, sigma, eta, selectedInput, ii, jj)

def e(N, w, sigma, eta, selectedInput, ii, jj):
    dist = min(np.absolute(ii-jj) , np.absolute(np.absolute(ii-jj)-N)) # Centering the neighborhood function in the winner
    f(N, w, sigma, eta, selectedInput, jj, dist)

def f(N, w, sigma, eta, selectedInput, jj, dist):
    w[N - 2][jj,:] = w[N - 2][jj,:] + eta * A(dist,sigma) * (selectedInput - w[N - 2][jj,:]) # Updating the weights

a()

         6701974 function calls (6701853 primitive calls) in 4.985 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   238921    0.691    0.000    0.777    0.000 x.py:56(f)
   119461    0.613    0.000    4.923    0.000 x.py:43(d)
   119461    0.498    0.000    2.144    0.000 shape_base.py:686(kron)
   238921    0.462    0.000    1.280    0.000 x.py:52(e)
   119461    0.369    0.000    0.495    0.000 numeric.py:1015(outer)

这将 f() 标识为最大时间。