from - python graph histogram



numpy ufuncs速度vs循環速度 (1)

這是正常和預期的行為。 這太簡單了,無法應用“避免使用numpy循環”語句。 如果你正在處理內部循環(幾乎)總是如此。 但是如果是外部循環(比如你的情況),則會有更多的例外。 特別是如果替代方案是使用廣播,因為這通過使用更多的內存來加速您的操作。

只是為了避免使用numpy循環來添加一些背景:

NumPy數組存儲為具有c類型的連續數組。 Python的int不是一樣的C int ! 因此,無論何時迭代數組中的每個項目,都需要從數組中插入項目,將其轉換為Python int ,然後執行任何您想要的操作,最後可能需要再次將其轉換為ac整數(稱為裝箱和拆箱的價值)。 例如,你想要使用Python對數組中的項進行sum

import numpy as np
arr = np.arange(1000)
%%timeit
acc = 0
for item in arr:
    acc += item
# 1000 loops, best of 3: 478 µs per loop

你最好使用numpy:

%timeit np.sum(arr)
# 10000 loops, best of 3: 24.2 µs per loop

即使你將循環引入Python C代碼,你也遠離了Numpy的性能:

%timeit sum(arr)
# 1000 loops, best of 3: 387 µs per loop

這條規則可能有例外,但這些將會非常稀少。 至少只要有一些等效的numpy功能。 所以如果你要遍歷單個元素,那麼你應該使用numpy。

有時候一個簡單的Python循環就足夠了。 這不是廣告宣傳,但與Python功能相比,numpy函數有巨大的開銷。 例如考慮一個3元素數組:

arr = np.arange(3)
%timeit np.sum(arr)
%timeit sum(arr)

哪一個會更快?

解決方案:Python函數比numpy解決方案執行得更好:

# 10000 loops, best of 3: 21.9 µs per loop  <- numpy
# 100000 loops, best of 3: 6.27 µs per loop <- python

但是這與你的例子有什麼關係呢? 事實上並非如此,因為你總是在數組上使用numpy函數(而不是單個元素,甚至幾個元素),所以你的內部循環已經使用了優化的函數。 這就是為什麼兩者執行大致相同(+/-大約500個元素的元素為非常少的因素10的因素2)。 但是這不是真正的循環開銷,這是函數調用開銷!

你的循環解決方案

使用行剖析器resolution = 100

def fun_func(tim, prec, values):
    for i, ti in enumerate(tim):
        values[i] = np.sum(np.sin(prec * ti))
%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2       101          752      7.4      5.7      for i, ti in enumerate(tim):
     3       100        12449    124.5     94.3          values[i] = np.sum(np.sin(prec * ti))

95%被用在循環中,我甚至將循環體分成幾個部分來驗證這一點:

def fun_func(tim, prec, values):
    for i, ti in enumerate(tim):
        x = prec * ti
        x = np.sin(x)
        x = np.sum(x)
        values[i] = x
%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2       101          609      6.0      3.5      for i, ti in enumerate(tim):
     3       100         4521     45.2     26.3          x = prec * ti
     4       100         4646     46.5     27.0          x = np.sin(x)
     5       100         6731     67.3     39.1          x = np.sum(x)
     6       100          714      7.1      4.1          values[i] = x

消費者的時間是np.multiplynp.sinnp.sum ,因為您可以通過比較每次通話的時間和開銷來輕鬆檢查:

arr = np.ones(1, float)
%timeit np.sum(arr)
# 10000 loops, best of 3: 22.6 µs per loop

所以只要計算運行時間的計算函數調用開銷很小,就會有類似的運行時間。 即使有100個項目,你也非常接近開銷時間。 訣竅是知道他們在什麼時候盈虧平衡。 有1000個項目的呼叫開銷仍然很大:

%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2      1001         5864      5.9      2.4      for i, ti in enumerate(tim):
     3      1000        42817     42.8     17.2          x = prec * ti
     4      1000       119327    119.3     48.0          x = np.sin(x)
     5      1000        73313     73.3     29.5          x = np.sum(x)
     6      1000         7287      7.3      2.9          values[i] = x

但是在resolution = 5000的情況下,與運行時相比,開銷非常低:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2      5001        29412      5.9      0.9      for i, ti in enumerate(tim):
     3      5000       388827     77.8     11.6          x = prec * ti
     4      5000      2442460    488.5     73.2          x = np.sin(x)
     5      5000       441337     88.3     13.2          x = np.sum(x)
     6      5000        36187      7.2      1.1          values[i] = x

當你花費500us在每個np.sin調用你不關心20us的開銷了。

一個小心的話可能是為了: line_profiler可能包含一些額外的開銷每行,也許也是每個函數調用,所以函數調用開銷變得可忽略的點可能會更低!

您的廣播解決方案

我開始分析第一個解決方案,讓我們對第二個解決方案做同樣的事情:

def fun_func(tim, prec, values):
    x = tim[:, np.newaxis]
    x = x * prec
    x = np.sin(x)
    x = np.sum(x, axis=1)
    return x

再次使用resolution=100的line_profiler:

%lprun -f fun_func fun_func(tim, prec, values)
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2         1           27     27.0      0.5      x = tim[:, np.newaxis]
     3         1          638    638.0     12.9      x = x * prec
     4         1         3963   3963.0     79.9      x = np.sin(x)
     5         1          326    326.0      6.6      x = np.sum(x, axis=1)
     6         1            4      4.0      0.1      return x

這已經大大超過了開銷時間,因此我們比循環更快地結束10倍。

我也做了resolution=1000的分析:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2         1           28     28.0      0.0      x = tim[:, np.newaxis]
     3         1        17716  17716.0     14.6      x = x * prec
     4         1        91174  91174.0     75.3      x = np.sin(x)
     5         1        12140  12140.0     10.0      x = np.sum(x, axis=1)
     6         1           10     10.0      0.0      return x

precision=5000

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def fun_func(tim, prec, values):
     2         1           34     34.0      0.0      x = tim[:, np.newaxis]
     3         1       333685 333685.0     11.1      x = x * prec
     4         1      2391812 2391812.0    79.6      x = np.sin(x)
     5         1       280832 280832.0      9.3      x = np.sum(x, axis=1)
     6         1           14     14.0      0.0      return x

1000的尺寸仍然更快,但正如我們所見,在環路解決方案中,通話開銷仍然不可忽視。 但是對於resolution = 5000 ,每一步花費的時間幾乎是相同的(有些比較慢,其他比較快,但總體上相當類似)

另一個影響是,當你進行乘法運算時,實際的廣播變得很重要。 即使有非常聰明的numpy解決方案,這仍然包括一些額外的計算。 對於resolution=10000您可以看到廣播乘法開始佔用更多與循環解決方案有關的“%時間”:

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     1                                           def broadcast_solution(tim, prec, values):
     2         1           37     37.0      0.0      x = tim[:, np.newaxis]
     3         1      1783345 1783345.0    13.9      x = x * prec
     4         1      9879333 9879333.0    77.1      x = np.sin(x)
     5         1      1153789 1153789.0     9.0      x = np.sum(x, axis=1)
     6         1           11     11.0      0.0      return x


Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     8                                           def loop_solution(tim, prec, values):
     9     10001        62502      6.2      0.5      for i, ti in enumerate(tim):
    10     10000      1287698    128.8     10.5          x = prec * ti
    11     10000      9758633    975.9     79.7          x = np.sin(x)
    12     10000      1058995    105.9      8.6          x = np.sum(x)
    13     10000        75760      7.6      0.6          values[i] = x

但除了實際所花的時間之外,還有另一件事情是內存消耗。 你的循環解決方案需要O(n)內存,因為你總是處理n元素。 然而,廣播解決方案需要O(n*n)存儲器。 如果在循環中使用resolution=20000 ,則可能需要等待一段時間,但它仍然只需要8bytes/element * 20000 element ~= 160kB但在廣播中需要8bytes/element * 20000 element ~= 160kB 。 這忽略了常數因素(如臨時數組又名中間數組)。 假設你走得更遠,你會非常快地耗盡內存!

時間再次總結一下:

  • 如果你在一個numpy數組中的單個項目做一個python循環,那麼你做錯了。
  • 如果循環遍歷numpy數組的子陣列,請確保每個循環中的函數調用開銷與函數耗用的時間相比可以忽略不計。
  • 如果你廣播numpy數組,確保你沒有用完內存。

但關於優化的最重要的一點仍然是:

  • 如果速度太慢,只能優化代碼! 如果速度太慢,那麼只有在分析代碼後才進行優化。

  • 不要盲目信任簡化的陳述,而且不要盲目地優化沒有分析。

最後一個想法是:

如果在numpy或scipy中沒有已經存在的解決方案,則可以使用cython , numba或numexpr輕鬆實現這些需要循環或廣播的函數 。

例如,一個將循環解決方案的內存效率與resolutions的廣播解決方案的速度相結合的numba函數看起來像這樣:

from numba import njit

import math

@njit
def numba_solution(tim, prec, values):
    size = tim.size
    for i in range(size):
        ti = tim[i]
        x = 0
        for j in range(size):
            x += math.sin(prec[j] * ti)
        values[i] = x

正如在註釋中指出的那樣, numexpr也可以非常快速地評估廣播計算, 而不需要O(n*n)內存:

>>> import numexpr
>>> tim_2d = tim[:, np.newaxis]
>>> numexpr.evaluate('sum(sin(tim_2d * prec), axis=1)')

我讀了很多“避免與numpy循環”。 所以,我試了。 我正在使用這個代碼(簡化版)。 一些輔助數據:

 In[1]: import numpy as np
        resolution = 1000                             # this parameter varies
        tim = np.linspace(-np.pi, np.pi, resolution) 
        prec = np.arange(1, resolution + 1)
        prec = 2 * prec - 1
        values = np.zeros_like(tim)

我的第一個實現是for循環:

 In[2]: for i, ti in enumerate(tim):
           values[i] = np.sum(np.sin(prec * ti))

然後,我擺脫了明確for循環,並實現了這一點:

 In[3]: values = np.sum(np.sin(tim[:, np.newaxis] * prec), axis=1)

而這個解決方案對於小數組來說更快,但是當我放大時,我有這樣的時間依賴性:

我錯過了什麼或是正常的行為? 如果不是,在哪裡挖?

編輯 :根據意見,這裡是一些額外的信息。 用IPython的%timeit%%timeit來衡量時間,每次運行都在新鮮的內核上執行。 我的筆記本電腦是宏碁Aspire v7-482pg(i7,8GB) 。 我在用著:

  • python 3.5.2
  • numpy 1.11.2 + mkl
  • Windows 10




numpy-ufunc