如何在向量化NumPy數組上進行移動窗口

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

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

滑動窗口操作非常普遍，非常有用。它們也很容易在Python中實現。學習如何實現移動窗口將把你的數據分析和爭論技能提升到一個新的水平。

什麼是滑動窗？

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

就是這樣。這就是滑動窗口的基本原理。當然，事情可能變得更加復雜。有限差分方法可以用於時間和空間數據。邏輯可以實現。可以使用更大的窗口大小或非正方形窗口。你懂的。但在其核心，移動窗口分析可以簡單地總結為鄰居元素的平均值。

需要註意的是，必須為邊緣元素設置特殊的調整，因為它們沒有9個相鄰元素。因此，許多分析都排除瞭邊緣元素。為簡單起見，我們將在本文中排除邊緣元素。

樣例數組

3×3的滑動窗口

創建一個NumPy數組

為瞭實現一些簡單的示例，讓我們創建上面所示的數組。首先，導入numpy。

import numpy as np

然後使用arange創建一個7×7的數組，值范圍從1到48。另外，創建另一個包含無數據值的數組，該數組的形狀和數據類型與初始數組相同。在本例中，我使用-1作為無數據值。

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

我們將使用這些數組來開發下面的滑動窗口示例。

通過循環實現滑動窗口

毫無疑問，你已經聽說過Python中的循環很慢，應該盡可能避免。特別是在使用大型NumPy數組時。這是完全正確。盡管如此，我們將首先看一個使用循環的示例，因為這是一種簡單的方法來概念化在移動窗口操作中發生的事情。在你通過循環示例掌握瞭概念之後，我們將繼續使用更有效的向量化方法。

要實現移動窗口，隻需循環遍歷所有內部數組元素，識別所有相鄰元素的值，並在特定的計算中使用這些值。

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

行偏移

列偏移

循環中NumPy移動窗口的Python代碼

我們可以用三行代碼實現一個移動窗口。這個例子在滑動窗口內計算平均值。首先，循環遍歷數組的內部行。其次，循環遍歷數組的內部列。第三，在滑動窗口內計算平均值，並將值賦給輸出數組中相應的數組元素。

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

循環後結果

你將註意到結果與輸入數組具有相同的值，但是外部元素沒有被分配數據值，因為它們不包含9個相鄰元素。

[[-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.]]

向量化滑動窗口

Python中的數組循環通常計算效率低下。通過對通常在循環中執行的操作進行向量化，可以提高效率。移動窗口矢量化可以通過同時抵消數組內部的所有元素來實現。

如下圖所示。每個圖像都有相應的索引。你將註意到最後一張圖像索引瞭所有內部元素，並且對應的圖像索引瞭每個相鄰元素的偏移量。

從左到右的偏移索引：[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數組上的向量化移動窗口的Python代碼

有瞭上述偏移量，我們現在可以輕松地在一行代碼中實現滑動窗口。 隻需將輸出數組的所有內部元素設置為根據相鄰元素計算所需輸出的函數。

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

矢量化滑動窗口結果

如你所見，這將得到與循環相同的結果。

[[-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.]]

速度比較

上述兩種方法產生相同的結果，但哪一種更有效?我計算瞭從5行到100列的數組的每種方法的速度。每種方法對每個測試100次。下面是每種方法的平均時間。

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

總結

移動窗口計算在許多數據分析工作流程中非常常見。這些計算是非常有用的，非常容易實現。然而，使用循環來實現滑動窗口操作是非常低效的。

向量化的移動窗口實現不僅更高效，而且使用更少的代碼行。一旦掌握瞭實現滑動窗口的向量化方法，就可以輕松有效地提高工作流程的速度。

補充：Python學習筆記——Numpy數組的移動滑窗，使用as_strided實現

Numpy中移動滑窗的實現

為何需要移動滑窗

在量化投資分析過程中，對歷史數據進行分析是一個必不可少的步驟。滑窗在歷史數據分析中的重要性不言而喻。譬如移動平均、指數平滑移動平均、MACD、DMA等等價格指標的計算都無一例外需要用到滑窗。

作為一種非常受歡迎的數據分析工具，pandas中提供瞭專門的滑窗類：DataFrame.rolling()。通過這個滑窗類，可以非常容易地實現移動平均等等算法，但是，在某些情況下，Pandas的運行速度還是不夠，需要借助Numpy的高效率進一步提升速度，這時候就需要在Numpy中實現滑窗瞭。

Numpy中的移動滑窗

可惜Numpy並沒有提供直接簡單的滑窗方法，如果使用for-loop來實現滑窗，不僅效率打折扣，而且內存占用也非常大。實際上，Numpy提供瞭一個非常底層的函數可以用來生成滑窗：Numpy.lib.stride_tricks.as_stried。

移動滑窗的as_strided實現方法

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

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

如果使用as_strided函數生成上述滑窗，需要用下面的代碼，它生成一個三維數組，包括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)

我們再嘗試一下用for-loop的方法生成一個滑窗檢驗一下前面生成的滑窗是否正確：

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的十萬倍以上！那麼as_strided是如何做到的呢？

關於as_strided函數的詳細解析

as_strided是怎麼回事呢？看它的函數解釋：

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

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

shape的含義非常簡單，就是指輸出的數據的行、列、層數，這個參數是一個元組，元組的元素數量等於數組的維度。

而stride的含義就相對復雜一些，其實它的含義是指“步幅”，意思是每一個維度的數據在內存上平移的字節數量。

因為數組在內存中的存放方式是一維線性方式存放的，因此要訪問數組中的某個數字就需要知道平移到哪一個內存單元，ndarray通過stride“步幅”來指定這個平移的幅度。

在as_strided函數中，stride也是一個元組，其元素的數量必須跟shape的元素數量相同，每一個元素就代表該維度的每一個數據相對前一個數據的內存間隔。

舉個例子：

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]])

我們之所以看到一個二維數組，是因為numpy數組的shape為(5, 3)，stride為(24, 8)，意思是說，我們看到的數據有5行3列，對應shape的(5, 3)，每一行與前一行間隔24個字節（其實就是三個數字，因為每一個int類型占據8字節，而每一列數字比前一列相差8字節（1個數字）

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

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

接下來需要確定stride，如前所述stride同樣是一個包含三個元素的元組，第一個元素是兩層數據之間的內存間隔，由於我們的滑窗每滑動一次下移一行，因此層stride應該是平移三個數字，也就是24個字節，行stride和列stride與原來的行列stride一致，因為我們需要原樣看到按順序的數字，因此，新的stride就是：(24, 24, 8)

我們來看看這個新的數據視圖是什麼樣子：

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]]])

看！一個數據滑窗正確地出現瞭！

使用as_strided函數的危險之處

使用s_strided函數的最大問題是內存讀取風險，在as_strided生成新的視圖時，由於直接操作內存地址（這一點像極瞭C的指針操作），而且它並不會檢查內存地址是否越界，因此如果稍有不慎，就會讀到別的內存地址。關鍵是，如果不設置可讀參數，還能直接對內存中的數據進行操作，這樣就帶來瞭無比大的風險。瞭解這個風險對正確操作至關重要！

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

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]]])

這時對象的第五組就映射到瞭三個未知的內存地址上，如果不慎修改瞭這三個地址上的內容，就可能造成難以預料的問題，如程序崩潰等。

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

以上為個人經驗，希望能給大傢一個參考，也希望大傢多多支持WalkonNet。

推薦閱讀：

• 聊聊Numpy.array中[:]和[::]的區別在哪
• Numpy ndarray 多維數組對象的使用
• Python Numpy中ndarray的常見操作
• python數學建模之Numpy 應用介紹與Pandas學習
• Python機器學習三大件之一numpy