圖解 NumPy：常用函數(shù)的內(nèi)在機(jī)制！
?△點(diǎn)擊上方“Python貓”關(guān)注 ，回復(fù)“1”領(lǐng)取電子書
選自Medium
作者：Lev Maximov
機(jī)器之心
支持大量多維數(shù)組和矩陣運(yùn)算的 NumPy 軟件庫是許多機(jī)器學(xué)習(xí)開發(fā)者和研究者的必備工具，本文將通過直觀易懂的圖示解析常用的 NumPy 功能和函數(shù)，幫助你理解 NumPy 操作數(shù)組的內(nèi)在機(jī)制。
NumPy 是一個(gè)基礎(chǔ)軟件庫，很多常用的 Python 數(shù)據(jù)處理軟件庫都使用了它或受到了它的啟發(fā)，包括 pandas、PyTorch、TensorFlow、Keras 等。理解 NumPy 的工作機(jī)制能夠幫助你提升在這些軟件庫方面的技能。而且在 GPU 上使用 NumPy 時(shí)，無需修改或僅需少量修改代碼。
NumPy 的核心概念是 n 維數(shù)組。n 維數(shù)組的美麗之處是大多數(shù)運(yùn)算看起來都一樣，不管數(shù)組有多少維。但一維和二維有點(diǎn)特殊。本文分為三部分：
1. 向量：一維數(shù)組
2. 矩陣：二維數(shù)組
3. 三維及更高維
本文參考了 Jay Alammar 的文章《A Visual Intro to NumPy》并將其作為起點(diǎn)，然后進(jìn)行了擴(kuò)充，并做了一些細(xì)微修改。
NumPy 數(shù)組和 Python 列表
乍一看，NumPy 數(shù)組與 Python 列表類似。它們都可作為容器，能夠快速獲取和設(shè)置元素，但插入和移除元素會稍慢一些。
NumPy 數(shù)組完勝列表的最簡單例子是算術(shù)運(yùn)算：
除此之外，NumPy 數(shù)組的優(yōu)勢和特點(diǎn)還包括：
更緊湊，尤其是當(dāng)維度大于一維時(shí)；
當(dāng)運(yùn)算可以向量化時(shí)，速度比列表更快；
當(dāng)在后面附加元素時(shí)，速度比列表慢；
通常是同質(zhì)的：當(dāng)元素都是一種類型時(shí)速度很快。
這里 O(N) 的意思是完成該運(yùn)算所需的時(shí)間和數(shù)組的大小成正比，而 O*(1)（即所謂的「均攤 O(1)」）的意思是完成運(yùn)算的時(shí)間通常與數(shù)組的大小無關(guān)。
向量：一維數(shù)組
向量初始化
為了創(chuàng)建 NumPy 數(shù)組，一種方法是轉(zhuǎn)換 Python 列表。NumPy 數(shù)組類型可以直接從列表元素類型推導(dǎo)得到。
要確保向其輸入的列表是同一種類型，否則你最終會得到 dtype=’object’，這會影響速度，最終只留下 NumPy 中含有的語法糖。
NumPy 數(shù)組不能像 Python 列表一樣增長。數(shù)組的末端沒有留下任何便于快速附加元素的空間。因此，常見的做法是要么先使用 Python 列表，準(zhǔn)備好之后再將其轉(zhuǎn)換為 NumPy 數(shù)組，要么是使用 np.zeros 或 np.empty 預(yù)先留下必要的空間：
通常我們有必要?jiǎng)?chuàng)建在形狀和元素類型上與已有數(shù)組匹配的空數(shù)組。
事實(shí)上，所有用于創(chuàng)建填充了常量值的數(shù)組的函數(shù)都帶有 _like 的形式：
NumPy 中有兩個(gè)函數(shù)能用單調(diào)序列執(zhí)行數(shù)組初始化：
如果你需要類似 [0., 1., 2.] 這樣的浮點(diǎn)數(shù)數(shù)組，你可以修改 arange 輸出的類型：arange(3).astype(float)，但還有一種更好的方法。arange 函數(shù)對類型很敏感：如果你以整型數(shù)作為參數(shù)輸入，它會生成整型數(shù)；如果你輸入浮點(diǎn)數(shù)（比如 arange(3.)），它會生成浮點(diǎn)數(shù)。
但 arange 并不非常擅長處理浮點(diǎn)數(shù)：
在我們眼里，這個(gè) 0.1 看起來像是一個(gè)有限的十進(jìn)制數(shù)，但計(jì)算機(jī)不這么看。在二進(jìn)制表示下，0.1 是一個(gè)無限分?jǐn)?shù)，因此必須進(jìn)行約分，也由此必然會產(chǎn)生誤差。也因?yàn)檫@個(gè)原因，如果向 arange 函數(shù)輸入帶分?jǐn)?shù)部分的 step，通常得不到什么好結(jié)果：你可能會遇到差一錯(cuò)誤 (off-by-one error)。你可以使該區(qū)間的末端落在一個(gè)非整數(shù)的 step 數(shù)中（solution1），但這會降低代碼的可讀性和可維護(hù)性。這時(shí)候，linspace 就可以派上用場了。它不受舍入的影響，總能生成你要求的元素?cái)?shù)值。不過，使用 linspace 時(shí)會遇到一個(gè)常見的陷阱：它統(tǒng)計(jì)的是數(shù)據(jù)點(diǎn)的數(shù)量，而不是區(qū)間，因此其最后一個(gè)參數(shù) num 通常比你所想的數(shù)大 1。因此，上面最后一個(gè)例子中的數(shù)是 11，而不是 10。
在進(jìn)行測試時(shí)，我們通常需要生成隨機(jī)數(shù)組：
向量索引
一旦你的數(shù)組中有了數(shù)據(jù)，NumPy 就能以非常巧妙的方式輕松地提供它們：
除了「花式索引（fancy indexing）」外，上面給出的所有索引方法都被稱為「view」：它們并不存儲數(shù)據(jù)，也不會在數(shù)據(jù)被索引后發(fā)生改變時(shí)反映原數(shù)組的變化情況。
所有包含花式索引的方法都是可變的：它們允許通過分配來修改原始數(shù)組的內(nèi)容，如上所示。這一功能可通過將數(shù)組切分成不同部分來避免總是復(fù)制數(shù)組的習(xí)慣。
Python 列表與 NumPy 數(shù)組的對比
為了獲取 NumPy 數(shù)組中的數(shù)據(jù)，另一種超級有用的方法是布爾索引（boolean indexing），它支持使用各類邏輯運(yùn)算符：
any 和 all 的作用與在 Python 中類似，但不會短路。
不過要注意，這里不支持 Python 的「三元比較」，比如 3<=a<=5。
如上所示，布爾索引也是可寫的。其兩個(gè)常用功能都有各自的專用函數(shù)：過度重載的 np.where 函數(shù)和 np.clip 函數(shù)。它們的含義如下：
向量運(yùn)算
NumPy 在速度上很出彩的一大應(yīng)用領(lǐng)域是算術(shù)運(yùn)算。向量運(yùn)算符會被轉(zhuǎn)換到 C++ 層面上執(zhí)行，從而避免緩慢的 Python 循環(huán)的成本。NumPy 支持像操作普通的數(shù)那樣操作整個(gè)數(shù)組。
與 Python 句法一樣，a//b 表示 a 除 b（除法的商），x**n 表示 x?。
正如加減浮點(diǎn)數(shù)時(shí)整型數(shù)會被轉(zhuǎn)換成浮點(diǎn)數(shù)一樣，標(biāo)量也會被轉(zhuǎn)換成數(shù)組，這個(gè)過程在 NumPy 中被稱為廣播（broadcast）。
大多數(shù)數(shù)學(xué)函數(shù)都有用于處理向量的 NumPy 對應(yīng)函數(shù)：
標(biāo)量積有自己的運(yùn)算符：
執(zhí)行三角函數(shù)時(shí)也無需循環(huán)：
我們可以在整體上對數(shù)組進(jìn)行舍入：
floor 為舍、ceil 為入，around 則是舍入到最近的整數(shù)（其中 .5 會被舍掉）
NumPy 也能執(zhí)行基礎(chǔ)的統(tǒng)計(jì)運(yùn)算：
NumPy 的排序函數(shù)沒有 Python 的排序函數(shù)那么強(qiáng)大：
Python 列表與 NumPy 數(shù)組的排序函數(shù)對比
在一維情況下，如果缺少 reversed 關(guān)鍵字，那么只需簡單地對結(jié)果再執(zhí)行反向，最終效果還是一樣。二維的情況則會更困難一些（人們正在請求這一功能）。
搜索向量中的元素
與 Python 列表相反，NumPy 數(shù)組沒有索引方法。人們很久之前就在請求這個(gè)功能，但一直還沒實(shí)現(xiàn)。
Python 列表與 NumPy 數(shù)組的對比，index() 中的方括號表示可以省略 j 或同時(shí)省略 i 和 j。
一種查找元素的方法是 np.where(a==x)[0][0]，但這個(gè)方法既不優(yōu)雅，速度也不快，因?yàn)樗枰獧z查數(shù)組中的所有元素，即便所要找的目標(biāo)就在數(shù)組起始位置也是如此。
另一種更快的方式是使用 Numba 來加速 next((i[0] for i, v in np.ndenumerate(a) if v==x), -1)。
一旦數(shù)組的排序完成，搜索就容易多了：v = np.searchsorted(a, x); return v if a[v]==x else -1 的速度很快，時(shí)間復(fù)雜度為 O(log N)，但它需要 O(N log N) 時(shí)間先排好序。
事實(shí)上，用 C 來實(shí)現(xiàn)它進(jìn)而加速搜索并不是問題。問題是浮點(diǎn)比較。這對任何數(shù)據(jù)來說都不是一種簡單直接可用的任務(wù)。
比較浮點(diǎn)數(shù)
函數(shù) np.allclose(a, b) 能在一定公差下比較浮點(diǎn)數(shù)數(shù)組。
函數(shù) np.allclose(a, b) 的工作過程示例。并沒有萬能方法！
np.allclose 假設(shè)所有被比較的數(shù)都在典型的 1 的范圍內(nèi)。舉個(gè)例子，如果要在納秒級的速度內(nèi)完成計(jì)算，則需要用默認(rèn)的 atol 參數(shù)值除以 1e9：np.allclose(1e-9, 2e-9, atol=1e-17) == False.
math.isclose 則不會對要比較的數(shù)進(jìn)行任何假設(shè)，而是依賴用戶給出合理的 abs_tol 值（對于典型的 1 的范圍內(nèi)的值，取默認(rèn)的 np.allclose atol 值 1e-8 就足夠好了）：math.isclose(0.1+0.2–0.3, abs_tol=1e-8)==True.
除此之外，np.allclose 在絕對值和相對公差的公式方面還有一些小問題，舉個(gè)例子，對于給定的 a 和 b，存在 allclose(a, b) != allclose(b, a)。這些問題已在（標(biāo)量）函數(shù) math.isclose 中得到了解決，我們將在后面介紹它。對于這方面的更多內(nèi)容，請參閱 GitHub 上的浮點(diǎn)數(shù)指南和對應(yīng)的 NumPy 問題（https://floating-point-gui.de/errors/comparison/）。
矩陣：二維數(shù)組
NumPy 曾有一個(gè)專門的 matrix 類，但現(xiàn)在已經(jīng)棄用了，所以本文會交替使用「矩陣」和「二維數(shù)組」這兩個(gè)術(shù)語。
矩陣的初始化句法與向量類似：
這里必須使用雙括號，因?yàn)榈诙€(gè)位置參數(shù)是 dtype（可選，也接受整數(shù)）。
隨機(jī)矩陣生成的句法也與向量的類似：
二維索引的句法比嵌套列表更方便：
view 符號的意思是當(dāng)切分一個(gè)數(shù)組時(shí)實(shí)際上沒有執(zhí)行復(fù)制。當(dāng)該數(shù)組被修改時(shí)，這些改變也會反映到切分得到的結(jié)果上。
axis 參數(shù)
在很多運(yùn)算中（比如 sum），你需要告訴 NumPy 是在列上還是行上執(zhí)行運(yùn)算。為了獲取適用于任意維度的通用符號，NumPy 引入了 axis 的概念：事實(shí)上，axis 參數(shù)的值是相關(guān)問題中索引的數(shù)量：第一個(gè)索引為 axis=0，第二個(gè)索引為 axis=1，以此類推。因此在二維情況下，axis=0 是按列計(jì)算，axis=1 是按行計(jì)算。
矩陣算術(shù)運(yùn)算
除了逐元素執(zhí)行的常規(guī)運(yùn)算符（比如 +、-、、/、//、*），這里還有一個(gè)計(jì)算矩陣乘積的 @ 運(yùn)算符：
我們已在第一部分介紹過標(biāo)量到數(shù)組的廣播，在其基礎(chǔ)上進(jìn)行泛化后，NumPy 支持向量和矩陣的混合運(yùn)算，甚至兩個(gè)向量之間的運(yùn)算：
二維數(shù)組中的廣播
行向量和列向量
正如上面的例子所示，在二維情況下，行向量和列向量的處理方式有所不同。這與具備某類一維數(shù)組的 NumPy 實(shí)踐不同（比如二維數(shù)組 a— 的第 j 列 a[:,j] 是一個(gè)一維數(shù)組）。默認(rèn)情況下，一維數(shù)組會被視為二維運(yùn)算中的行向量，因此當(dāng)用一個(gè)矩陣乘以一個(gè)行向量時(shí)，你可以使用形狀 (n,) 或 (1, n)——結(jié)果是一樣的。如果你需要一個(gè)列向量，則有多種方法可以基于一維數(shù)組得到它，但出人意料的是「轉(zhuǎn)置」不是其中之一。
基于一維數(shù)組得到二維數(shù)組的運(yùn)算有兩種：使用 reshape 調(diào)整形狀和使用 newaxis 進(jìn)行索引：
其中 -1 這個(gè)參數(shù)是告訴 reshape 自動(dòng)計(jì)算其中一個(gè)維度大小，方括號中的 None 是用作 np.newaxis 的快捷方式，這會在指定位置添加一個(gè)空 axis。
因此，NumPy 共有三類向量：一維向量、二維行向量和二維列向量。下圖展示了這三種向量之間的轉(zhuǎn)換方式：
一維向量、二維行向量和二維列向量之間的轉(zhuǎn)換方式。根據(jù)廣播的原則，一維數(shù)組可被隱含地視為二維行向量，因此通常沒必要在這兩者之間執(zhí)行轉(zhuǎn)換——因此相應(yīng)的區(qū)域被陰影化處理。
矩陣操作
合并數(shù)組的函數(shù)主要有兩個(gè)：
這兩個(gè)函數(shù)適用于只堆疊矩陣或只堆疊向量，但當(dāng)需要堆疊一維數(shù)組和矩陣時(shí)，只有 vstack 可以奏效：hstack 會出現(xiàn)維度不匹配的錯(cuò)誤，原因如前所述，一維數(shù)組會被視為行向量，而不是列向量。針對這個(gè)問題，解決方法要么是將其轉(zhuǎn)換為行向量，要么是使用能自動(dòng)完成這一操作的 column_stack 函數(shù)：
堆疊的逆操作是拆分：
復(fù)制矩陣的方法有兩種：復(fù)制 - 粘貼式的 tile 和分頁打印式的 repeat：
delete 可以刪除特定的行和列：
刪除的逆操作為插入，即 insert：
append 函數(shù)就像 hstack 一樣，不能自動(dòng)對一維數(shù)組執(zhí)行轉(zhuǎn)置，因此同樣地，要么需要改變該向量的形狀，要么就需要增加一個(gè)維度，或者使用 column_stack：
事實(shí)上，如果你只需要向數(shù)組的邊緣添加常量值，那么（稍微復(fù)雜的）pad 函數(shù)應(yīng)該就足夠了：
網(wǎng)格
廣播規(guī)則使得我們能更簡單地操作網(wǎng)格。假設(shè)你有如下矩陣（但非常大）：
使用 C 和使用 Python 創(chuàng)建矩陣的對比
這兩種方法較慢，因?yàn)樗鼈儠褂?Python 循環(huán)。為了解決這樣的問題，MATLAB 的方式是創(chuàng)建一個(gè)網(wǎng)格：
使用 MATLAB 創(chuàng)建網(wǎng)格的示意圖
使用如上提供的參數(shù) I 和 J，meshgrid 函數(shù)接受任意的索引集合作為輸入，mgrid 只是切分，indices 只能生成完整的索引范圍，fromfunction 只會調(diào)用所提供的函數(shù)一次。
但實(shí)際上，NumPy 中還有一種更好的方法。我們沒必要將內(nèi)存耗在整個(gè) I 和 J 矩陣上。存儲形狀合適的向量就足夠了，廣播規(guī)則可以完成其余工作。
使用 NumPy 創(chuàng)建網(wǎng)格的示意圖
沒有 indexing=’ij’ 參數(shù)，meshgrid 會改變這些參數(shù)的順序：J, I= np.meshgrid(j, i)——這是一種 xy 模式，對可視化 3D 圖表很有用。
除了在二維或三維網(wǎng)格上初始化函數(shù)，網(wǎng)格也可用于索引數(shù)組：
使用 meshgrid 索引數(shù)組，也適用于稀疏網(wǎng)格。
獲取矩陣統(tǒng)計(jì)數(shù)據(jù)
和 sum 一樣，min、max、argmin、argmax、mean、std、var 等所有其它統(tǒng)計(jì)函數(shù)都支持 axis 參數(shù)并能據(jù)此完成統(tǒng)計(jì)計(jì)算：
三個(gè)統(tǒng)計(jì)函數(shù)示例，為了避免與 Python 的 min 沖突，NumPy 中對應(yīng)的函數(shù)名為 np.amin。
用于二維及更高維的 argmin 和 argmax 函數(shù)會返回最小和最大值的第一個(gè)實(shí)例，在返回展開的索引上有點(diǎn)麻煩。為了將其轉(zhuǎn)換成兩個(gè)坐標(biāo)，需要使用 unravel_index 函數(shù)：
使用 unravel_index 函數(shù)的示例
all 和 any 函數(shù)也支持 axis 參數(shù)：
使用 all 和 any 函數(shù)的示例
矩陣排序
axis 參數(shù)雖然對上面列出的函數(shù)很有用，但對排序毫無用處：
使用 Python 列表和 NumPy 數(shù)組執(zhí)行排序的比較
這通常不是你在排序矩陣或電子表格時(shí)希望看到的結(jié)果：axis 根本不能替代 key 參數(shù)。但幸運(yùn)的是，NumPy 提供了一些支持按列排序的輔助函數(shù)——或有需要的話可按多列排序：
1. a[a[:,0].argsort()] 可按第一列對數(shù)組排序：
這里 argsort 會返回原數(shù)組排序后的索引的數(shù)組。
這個(gè)技巧可以重復(fù)，但必須謹(jǐn)慎，別讓下一次排序擾亂上一次排序的結(jié)果：
a = a[a[:,2].argsort()]
a = a[a[:,1].argsort(kind='stable')]
a = a[a[:,0].argsort(kind='stable')]
2. lexsort 函數(shù)能使用上述方式根據(jù)所有列進(jìn)行排序，但它總是按行執(zhí)行，而且所要排序的行的順序是反向的（即自下而上），因此使用它時(shí)會有些不自然，比如
- a[np.lexsort(np.flipud(a[2,5].T))] 會首先根據(jù)第 2 列排序，然后（當(dāng)?shù)?2 列的值相等時(shí)）再根據(jù)第 5 列排序。
– a[np.lexsort(np.flipud(a.T))] 會從左向右根據(jù)所有列排序。
這里，flipud 會沿上下方向翻轉(zhuǎn)該矩陣（準(zhǔn)確地說是 axis=0 方向，與 a[::-1,...] 一樣，其中三個(gè)點(diǎn)表示「所有其它維度」，因此翻轉(zhuǎn)這個(gè)一維數(shù)組的是突然的 flipud，而不是 fliplr。
3. sort 還有一個(gè) order 參數(shù)，但如果一開始是普通的（非結(jié)構(gòu)化）數(shù)組，它執(zhí)行起來既不快，也不容易使用。
4. 在 pandas 中執(zhí)行它可能是更好的選擇，因?yàn)樵?pandas 中，該特定運(yùn)算的可讀性要高得多，也不那么容易出錯(cuò)：
– pd.DataFrame(a).sort_values(by=[2,5]).to_numpy() 會先根據(jù)第 2 列排序，然后根據(jù)第 5 列排序。
– pd.DataFrame(a).sort_values().to_numpy() 會從左向右根據(jù)所有列排序。
三維及更高維
當(dāng)你通過調(diào)整一維向量的形狀或轉(zhuǎn)換嵌套的 Python 列表來創(chuàng)建 3D 數(shù)組時(shí)，索引的含義是 (z,y,x)。第一個(gè)索引是平面的數(shù)量，然后是在該平面上的坐標(biāo)：
展示 (z,y,x) 順序的示意圖
這個(gè)索引順序很方便，舉個(gè)例子，它可用于保存一些灰度圖像：a[i] 是索引第 i 張圖像的快捷方式。
但這個(gè)索引順序不是通用的。當(dāng)操作 RGB 圖像時(shí)，通常會使用 (y,x,z) 順序：首先是兩個(gè)像素坐標(biāo)，最后一個(gè)是顏色坐標(biāo)（Matplotlib 中是 RGB，OpenCV 中是 BGR）：
展示 (y,x,z) 順序的示意圖
這樣，我們就能很方便地索引特定的像素：a[i,j] 能提供 (i,j) 位置的 RGB 元組。
因此，創(chuàng)建幾何形狀的實(shí)際命令取決于你所在領(lǐng)域的慣例：
創(chuàng)建一般的三維數(shù)組和 RGB 圖像
很顯然，hstack、vstack、dstack 這些函數(shù)不支持這些慣例。它們硬編碼了 (y,x,z) 的索引順序，即 RGB 圖像的順序：
NumPy 使用 (y,x,z) 順序的示意圖，堆疊 RGB 圖像（這里僅有兩種顏色）
如果你的數(shù)據(jù)布局不同，使用 concatenate 命令來堆疊圖像會更方便一些，向一個(gè) axis 參數(shù)輸入明確的索引數(shù)值：
堆疊一般三維數(shù)組
如果你不習(xí)慣思考 axis 數(shù)，你可以將該數(shù)組轉(zhuǎn)換成 hstack 等函數(shù)中硬編碼的形式：
將數(shù)組轉(zhuǎn)換為 hstack 中硬編碼的形式的示意圖
這種轉(zhuǎn)換的成本很低：不會執(zhí)行實(shí)際的復(fù)制，只是執(zhí)行過程中混合索引的順序。
另一種可以混合索引順序的運(yùn)算是數(shù)組轉(zhuǎn)置。了解它可能會讓你更加熟悉三維數(shù)組。根據(jù)你決定使用的 axis 順序的不同，轉(zhuǎn)置數(shù)組所有平面的實(shí)際命令會有所不同：對于一般數(shù)組，它會交換索引 1 和 2，對 RGB 圖像而言是 0 和 1：
轉(zhuǎn)置一個(gè)三維數(shù)據(jù)的所有平面的命令
不過有趣的是，transpose 的默認(rèn) axes 參數(shù)（以及僅有的 a.T 運(yùn)算模式）會調(diào)轉(zhuǎn)索引順序的方向，這與上述兩個(gè)索引順序慣例都不相符。
最后，還有一個(gè)函數(shù)能避免你在處理多維數(shù)組時(shí)使用太多訓(xùn)練，還能讓你的代碼更簡潔——einsum（愛因斯坦求和）：
它會沿重復(fù)的索引對數(shù)組求和。在這個(gè)特定的例子中，np.tensordot(a, b, axis=1) 足以應(yīng)對這兩種情況，但在更復(fù)雜的情況中，einsum 的速度可能更快，而且通常也更容易讀寫——只要你理解其背后的邏輯。
如果你希望測試你的 NumPy 技能，GitHub 有 100 道相當(dāng)困難的練習(xí)題：https://github.com/rougier/numpy-100。
近期熱門文章推薦：
