怎么在向量化NumPy數(shù)組上進(jìn)行移動(dòng)窗口

在一些編輯算法中，我們可能會(huì)需要對(duì)一個(gè)二維數(shù)組（更準(zhǔn)確的稱(chēng)呼應(yīng)該為：矩陣）的某個(gè)區(qū)域進(jìn)行取值計(jì)算，在特定區(qū)域取完值后我們需要移動(dòng)到下一個(gè)區(qū)域。我們把這個(gè)區(qū)域稱(chēng)之為窗口，移動(dòng)這個(gè)區(qū)域我們又叫它移動(dòng)窗口。numpy作為一個(gè)出名的科學(xué)計(jì)算的庫(kù)，它提供了矩陣運(yùn)算的支持，所以他也提供了移動(dòng)窗口這種操作接下來(lái)這篇文章我們就來(lái)了解一下python在numpy上移動(dòng)窗口怎么操作吧！

簡(jiǎn)介

今天很有可能你已經(jīng)做了一些使用滑動(dòng)窗口(也稱(chēng)為移動(dòng)窗口)的事情，而你甚至不知道它。例如：許多編輯算法都是基于移動(dòng)窗口的。

在GIS中做地形分析的大多數(shù)地形柵格度量(坡度、坡向、山坡陰影等)都基于滑動(dòng)窗口。很多情況下，對(duì)格式化為二維數(shù)組的數(shù)據(jù)進(jìn)行分析時(shí)，都很有可能涉及到滑動(dòng)窗口。

滑動(dòng)窗口操作非常普遍，非常有用。它們也很容易在Python中實(shí)現(xiàn)。學(xué)習(xí)如何實(shí)現(xiàn)移動(dòng)窗口將把你的數(shù)據(jù)分析和爭(zhēng)論技能提升到一個(gè)新的水平。

什么是滑動(dòng)窗？

下面的例子顯示了一個(gè)3×3(3×3)滑動(dòng)窗口。用紅色標(biāo)注的數(shù)組元素是目標(biāo)元素。這是滑動(dòng)窗口將計(jì)算的新度量的數(shù)組位置。例如，在下面的圖像中，我們可以計(jì)算灰色窗口中9個(gè)元素的平均值(平均值也是8)，并將其分配給目標(biāo)元素，用紅色標(biāo)出。你可以計(jì)算最小值(0)、最大值(16)或其他一些指標(biāo)，而不是平均值。對(duì)數(shù)組中的每個(gè)元素都這樣做。

就是這樣。這就是滑動(dòng)窗口的基本原理。當(dāng)然，事情可能變得更加復(fù)雜。有限差分方法可以用于時(shí)間和空間數(shù)據(jù)。邏輯可以實(shí)現(xiàn)?？梢允褂酶蟮拇翱诖笮』蚍钦叫未翱凇Ｄ愣?。但在其核心，移動(dòng)窗口分析可以簡(jiǎn)單地總結(jié)為鄰居元素的平均值。

需要注意的是，必須為邊緣元素設(shè)置特殊的調(diào)整，因?yàn)樗鼈儧](méi)有9個(gè)相鄰元素。因此，許多分析都排除了邊緣元素。為簡(jiǎn)單起見(jiàn)，我們將在本文中排除邊緣元素。

樣例數(shù)組

3x3的滑動(dòng)窗口

創(chuàng)建一個(gè)NumPy數(shù)組

為了實(shí)現(xiàn)一些簡(jiǎn)單的示例，讓我們創(chuàng)建上面所示的數(shù)組。首先，導(dǎo)入numpy。

import numpy as np

然后使用arange創(chuàng)建一個(gè)7×7的數(shù)組，值范圍從1到48。另外，創(chuàng)建另一個(gè)包含無(wú)數(shù)據(jù)值的數(shù)組，該數(shù)組的形狀和數(shù)據(jù)類(lèi)型與初始數(shù)組相同。在本例中，我使用-1作為無(wú)數(shù)據(jù)值。

a = np.arange(49).reshape((7, 7)) 
b = np.full(a.shape, -1.0)

我們將使用這些數(shù)組來(lái)開(kāi)發(fā)下面的滑動(dòng)窗口示例。

通過(guò)循環(huán)實(shí)現(xiàn)滑動(dòng)窗口

毫無(wú)疑問(wèn)，你已經(jīng)聽(tīng)說(shuō)過(guò)Python中的循環(huán)很慢，應(yīng)該盡可能避免。特別是在使用大型NumPy數(shù)組時(shí)。這是完全正確。盡管如此，我們將首先看一個(gè)使用循環(huán)的示例，因?yàn)檫@是一種簡(jiǎn)單的方法來(lái)概念化在移動(dòng)窗口操作中發(fā)生的事情。在你通過(guò)循環(huán)示例掌握了概念之后，我們將繼續(xù)使用更有效的向量化方法。

要實(shí)現(xiàn)移動(dòng)窗口，只需循環(huán)遍歷所有內(nèi)部數(shù)組元素，識(shí)別所有相鄰元素的值，并在特定的計(jì)算中使用這些值。

通過(guò)行和列偏移量可以很容易地識(shí)別相鄰值。3×3窗口的偏移量如下所示。

行偏移

列偏移

循環(huán)中NumPy移動(dòng)窗口的Python代碼

我們可以用三行代碼實(shí)現(xiàn)一個(gè)移動(dòng)窗口。這個(gè)例子在滑動(dòng)窗口內(nèi)計(jì)算平均值。首先，循環(huán)遍歷數(shù)組的內(nèi)部行。其次，循環(huán)遍歷數(shù)組的內(nèi)部列。第三，在滑動(dòng)窗口內(nèi)計(jì)算平均值，并將值賦給輸出數(shù)組中相應(yīng)的數(shù)組元素。

for i in range(1, a.shape[0]-1):
    for j in range(1, a.shape[1]-1): 
        b[i, j] = (a[i-1, j-1] + a[i-1, j] + a[i-1, j+1] + a[i, j-1] + a[i, j] + a[i, j+1] + a[i+1, j-1] + a[i+1, j] + a[i+1, j+1]) / 9.0

循環(huán)后結(jié)果

你將注意到結(jié)果與輸入數(shù)組具有相同的值，但是外部元素沒(méi)有被分配數(shù)據(jù)值，因?yàn)樗鼈儾话?個(gè)相鄰元素。

[[-1. -1. -1. -1. -1. -1. -1.]
 [-1. 8. 9. 10. 11. 12. -1.]
 [-1. 15. 16. 17. 18. 19. -1.]
 [-1. 22. 23. 24. 25. 26. -1.]
 [-1. 29. 30. 31. 32. 33. -1.] 
 [-1. 36. 37. 38. 39. 40. -1.]
 [-1. -1. -1. -1. -1. -1. -1.]]

向量化滑動(dòng)窗口

Python中的數(shù)組循環(huán)通常計(jì)算效率低下。通過(guò)對(duì)通常在循環(huán)中執(zhí)行的操作進(jìn)行向量化，可以提高效率。移動(dòng)窗口矢量化可以通過(guò)同時(shí)抵消數(shù)組內(nèi)部的所有元素來(lái)實(shí)現(xiàn)。

如下圖所示。每個(gè)圖像都有相應(yīng)的索引。你將注意到最后一張圖像索引了所有內(nèi)部元素，并且對(duì)應(yīng)的圖像索引了每個(gè)相鄰元素的偏移量。

從左到右的偏移索引：[1：-1，：-2]，[1：-1，2：]，[2 :, 2：]

從左到右的偏移索引：[2 :,：-2]，[2 :, 1：-1]，[：-2，1：-1]

從左到右的偏移索引：[：-2，2：]，[：-2，：-2]，[1：-1、1：-1]

Numpy數(shù)組上的向量化移動(dòng)窗口的Python代碼

有了上述偏移量，我們現(xiàn)在可以輕松地在一行代碼中實(shí)現(xiàn)滑動(dòng)窗口。 只需將輸出數(shù)組的所有內(nèi)部元素設(shè)置為根據(jù)相鄰元素計(jì)算所需輸出的函數(shù)。

b[1:-1, 1:-1] = (a[1:-1, 1:-1] + a[:-2, 1:-1] + a[2:, 1:-1] + a[1:-1, :-2] + a[1:-1, 2:] + a[2:, 2:] + a[:-2, :-2] + a[2:, :-2] + a[:-2, 2:]) / 9.0

矢量化滑動(dòng)窗口結(jié)果

如你所見(jiàn)，這將得到與循環(huán)相同的結(jié)果。

[[-1. -1. -1. -1. -1. -1. -1.]
 [-1. 8. 9. 10. 11. 12. -1.]
 [-1. 15. 16. 17. 18. 19. -1.]
 [-1. 22. 23. 24. 25. 26. -1.]
 [-1. 29. 30. 31. 32. 33. -1.]
 [-1. 36. 37. 38. 39. 40. -1.]
 [-1. -1. -1. -1. -1. -1. -1.]]

速度比較

上述兩種方法產(chǎn)生相同的結(jié)果，但哪一種更有效?我計(jì)算了從5行到100列的數(shù)組的每種方法的速度。每種方法對(duì)每個(gè)測(cè)試100次。下面是每種方法的平均時(shí)間。

很明顯，向量化的方法更加有效。隨著數(shù)組大小的增加，循環(huán)的效率呈指數(shù)級(jí)下降。另外，需要注意的是，一個(gè)包含10,000個(gè)元素(100行和100列)的數(shù)組非常小。

總結(jié)

移動(dòng)窗口計(jì)算在許多數(shù)據(jù)分析工作流程中非常常見(jiàn)。這些計(jì)算是非常有用的，非常容易實(shí)現(xiàn)。然而，使用循環(huán)來(lái)實(shí)現(xiàn)滑動(dòng)窗口操作是非常低效的。

向量化的移動(dòng)窗口實(shí)現(xiàn)不僅更高效，而且使用更少的代碼行。一旦掌握了實(shí)現(xiàn)滑動(dòng)窗口的向量化方法，就可以輕松有效地提高工作流程的速度。

補(bǔ)充：Python學(xué)習(xí)筆記——Numpy數(shù)組的移動(dòng)滑窗，使用as_strided實(shí)現(xiàn)

Numpy中移動(dòng)滑窗的實(shí)現(xiàn)

為何需要移動(dòng)滑窗

在量化投資分析過(guò)程中，對(duì)歷史數(shù)據(jù)進(jìn)行分析是一個(gè)必不可少的步驟?；霸跉v史數(shù)據(jù)分析中的重要性不言而喻。譬如移動(dòng)平均、指數(shù)平滑移動(dòng)平均、MACD、DMA等等價(jià)格指標(biāo)的計(jì)算都無(wú)一例外需要用到滑窗。

作為一種非常受歡迎的數(shù)據(jù)分析工具，pandas中提供了專(zhuān)門(mén)的滑窗類(lèi)：DataFrame.rolling()。通過(guò)這個(gè)滑窗類(lèi)，可以非常容易地實(shí)現(xiàn)移動(dòng)平均等等算法，但是，在某些情況下，Pandas的運(yùn)行速度還是不夠，需要借助Numpy的高效率進(jìn)一步提升速度，這時(shí)候就需要在Numpy中實(shí)現(xiàn)滑窗了。

Numpy中的移動(dòng)滑窗

可惜Numpy并沒(méi)有提供直接簡(jiǎn)單的滑窗方法，如果使用for-loop來(lái)實(shí)現(xiàn)滑窗，不僅效率打折扣，而且內(nèi)存占用也非常大。實(shí)際上，Numpy提供了一個(gè)非常底層的函數(shù)可以用來(lái)生成滑窗：Numpy.lib.stride_tricks.as_stried。

移動(dòng)滑窗的as_strided實(shí)現(xiàn)方法

舉一個(gè)例子，首先生成一個(gè)5000行200列的二維數(shù)組，我們需要在這個(gè)二維數(shù)組上生成一個(gè)寬度為200的滑窗，也就是說(shuō)，第一個(gè)窗口包含前0～199行數(shù)據(jù)，第二個(gè)窗口包含1～200行，第三個(gè)窗口包含2～201行，以此類(lèi)推，一共4801組：

In [106]: d = np.random.randint(100, size=(5000,200))

如果使用as_strided函數(shù)生成上述滑窗，需要用下面的代碼，它生成一個(gè)三維數(shù)組，包括4801組200X200的矩陣，每一組200X200的矩陣代表一組滑窗：

In [107]: %timeit sd = as_strided(d, (4801,200,200), (200*8, 200*8, 8))
5.97 μs ± 33.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

我們?cè)賴(lài)L試一下用for-loop的方法生成一個(gè)滑窗檢驗(yàn)一下前面生成的滑窗是否正確：

In [108]: %%timeit
     ...: sd2 = np.zeros((4801,200,200))
     ...: for i in range(4801):
     ...:     sd2[i] = d[i:i+200]
     ...: 
722 ms ± 98.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [109]: np.allclose(sd, sd2)
Out[109]: True

從上面的代碼可以看出，使用as_strided生成一組滑窗，速度竟然是for-loop的十萬(wàn)倍以上！那么as_strided是如何做到的呢？

關(guān)于as_strided函數(shù)的詳細(xì)解析

as_strided是怎么回事呢？看它的函數(shù)解釋?zhuān)?/p>

Signature: as_strided(x, shape=None, strides=None, subok=False, writeable=True)
Docstring:
Create a view into the array with the given shape and strides.

.. warning:: This function has to be used with extreme care, see notes.

Parameters
----------
x : ndarray
Array to create a new.
shape : sequence of int, optional
The shape of the new array. Defaults to "x.shape".
strides : sequence of int, optional
The strides of the new array. Defaults to "x.strides".
subok : bool, optional
If True, subclasses are preserved.
writeable : bool, optional
If set to False, the returned array will always be readonly. Otherwise it will be writable if the original array was. It is advisable to set this to False if possible (see Notes).

Returns
-------
view : ndarray

這個(gè)函數(shù)接受的第一個(gè)參數(shù)是一個(gè)數(shù)組，第二個(gè)參數(shù)是輸出的數(shù)據(jù)shape，第三個(gè)參數(shù)是stride。要控制數(shù)據(jù)的輸出，shape和stride都非常重要

shape的含義非常簡(jiǎn)單，就是指輸出的數(shù)據(jù)的行、列、層數(shù)，這個(gè)參數(shù)是一個(gè)元組，元組的元素?cái)?shù)量等于數(shù)組的維度。

而stride的含義就相對(duì)復(fù)雜一些，其實(shí)它的含義是指“步幅”，意思是每一個(gè)維度的數(shù)據(jù)在內(nèi)存上平移的字節(jié)數(shù)量。

因?yàn)閿?shù)組在內(nèi)存中的存放方式是一維線性方式存放的，因此要訪問(wèn)數(shù)組中的某個(gè)數(shù)字就需要知道平移到哪一個(gè)內(nèi)存單元，ndarray通過(guò)stride“步幅”來(lái)指定這個(gè)平移的幅度。

在as_strided函數(shù)中，stride也是一個(gè)元組，其元素的數(shù)量必須跟shape的元素?cái)?shù)量相同，每一個(gè)元素就代表該維度的每一個(gè)數(shù)據(jù)相對(duì)前一個(gè)數(shù)據(jù)的內(nèi)存間隔。

舉個(gè)例子：

In [188]: d = np.random.randint(10, size=(5,3))

In [189]: d
Out[189]: 
array([[4, 4, 6],
       [2, 9, 3],
       [5, 1, 1],
       [2, 0, 0],
       [9, 2, 3]])

我們之所以看到一個(gè)二維數(shù)組，是因?yàn)閚umpy數(shù)組的shape為(5, 3)，stride為(24, 8)，意思是說(shuō)，我們看到的數(shù)據(jù)有5行3列，對(duì)應(yīng)shape的(5, 3)，每一行與前一行間隔24個(gè)字節(jié)（其實(shí)就是三個(gè)數(shù)字，因?yàn)槊恳粋€(gè)int類(lèi)型占據(jù)8字節(jié)，而每一列數(shù)字比前一列相差8字節(jié)（1個(gè)數(shù)字）

理解上面的含義以后，也就能理解如何生成一個(gè)數(shù)據(jù)滑窗了，如果我們需要生成一個(gè)2X3的數(shù)據(jù)滑窗，在d上滑動(dòng)，實(shí)際上可以生成一個(gè)4組，2行3列的數(shù)據(jù)視圖，第一組覆蓋d的第0、1兩行，第二層覆蓋d的第1、2兩行，第三層覆蓋d的第2、3兩行……這樣就形成了數(shù)據(jù)滑窗的效果，我們只要在新的數(shù)據(jù)視圖上遍歷，就能遍歷整個(gè)滑窗。這樣做的好處是，在整個(gè)遍歷的過(guò)程中完全不需要對(duì)數(shù)據(jù)進(jìn)行任何移動(dòng)或復(fù)制的操作，因此速度飛快。

根據(jù)上面的思路，我們需要生成一個(gè)新的數(shù)據(jù)視圖，其shape為(4, 2, 3)代表4組（從頭到尾滑動(dòng)4次），2行3列（滑窗的尺寸）

接下來(lái)需要確定stride，如前所述stride同樣是一個(gè)包含三個(gè)元素的元組，第一個(gè)元素是兩層數(shù)據(jù)之間的內(nèi)存間隔，由于我們的滑窗每滑動(dòng)一次下移一行，因此層stride應(yīng)該是平移三個(gè)數(shù)字，也就是24個(gè)字節(jié)，行stride和列stride與原來(lái)的行列stride一致，因?yàn)槲覀冃枰瓨涌吹桨错樞虻臄?shù)字，因此，新的stride就是：(24, 24, 8)

我們來(lái)看看這個(gè)新的數(shù)據(jù)視圖是什么樣子：

In [190]: as_strided(d, shape=(4,2,3), strides=(24,24,8))
Out[190]: 
array([[[4, 4, 6],
        [2, 9, 3]],

       [[2, 9, 3],
        [5, 1, 1]],

       [[5, 1, 1],
        [2, 0, 0]],

       [[2, 0, 0],
        [9, 2, 3]]])

看！一個(gè)數(shù)據(jù)滑窗正確地出現(xiàn)了！

使用as_strided函數(shù)的危險(xiǎn)之處

使用s_strided函數(shù)的最大問(wèn)題是內(nèi)存讀取風(fēng)險(xiǎn)，在as_strided生成新的視圖時(shí)，由于直接操作內(nèi)存地址（這一點(diǎn)像極了C的指針操作），而且它并不會(huì)檢查內(nèi)存地址是否越界，因此如果稍有不慎，就會(huì)讀到別的內(nèi)存地址。關(guān)鍵是，如果不設(shè)置可讀參數(shù)，還能直接對(duì)內(nèi)存中的數(shù)據(jù)進(jìn)行操作，這樣就帶來(lái)了無(wú)比大的風(fēng)險(xiǎn)。了解這個(gè)風(fēng)險(xiǎn)對(duì)正確操作至關(guān)重要！

例如，使用下面的stride會(huì)直接溢出到其他的未知內(nèi)存地址上，并讀取它的值，甚至還可以直接修改它：

In [194]: as_strided(d, shape=(5,2,3), strides=(24,24,8))
Out[194]: 
array([[[               4,                4,                6],
        [               2,                9,                3]],

       [[               2,                9,                3],
        [               5,                1,                1]],

       [[               5,                1,                1],
        [               2,                0,                0]],

       [[               2,                0,                0],
        [               9,                2,                3]],

       [[               9,                2,                3],
        [2251799813685248,            18963,                0]]])

這時(shí)對(duì)象的第五組就映射到了三個(gè)未知的內(nèi)存地址上，如果不慎修改了這三個(gè)地址上的內(nèi)容，就可能造成難以預(yù)料的問(wèn)題，如程序崩潰等。

所以，官方才在文檔中鄭重地警告：如果有可能，盡量避免使用as_strided函數(shù)

以上就是python在numpy上移動(dòng)窗口的全部?jī)?nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持W3Cschool。