知識星球
作者：王天慶
如需轉載請聯絡大資料（ID：hzdashuju）
Numpy是Numerical Python extensions 的縮寫，字面意思是Python數值計算擴充套件。Numpy是Python中眾多機器學習庫的依賴，這些庫透過Numpy實現基本的矩陣計算，Python的OpenCV庫自然也不例外。
Numpy支援高階、大量計算的矩陣、向量計算，與此同時提供了較為豐富的函式。Numpy採用友好的BSD許可協議開放原始碼。它是一個跨平臺的科學計算庫，提供了與Matlab相似的功能和操作方法。
雖然科學計算領域一直是Matlab的天下，但是Numpy基於更加現代化的程式語言——Python。而且Python憑藉著開源、免費、靈活性、簡單易學、工程特性好等特點風靡技術圈，已經成為機器學習、資料分析等領域的主流程式語言。
雖然Matlab提供的包非常多，但是Python因其簡單靈活、擴充套件性強等特點，也誕生了一系列優秀的庫。例如Scipy具有大多數Matlab所具備的功能，Matplotlib能夠簡便地進行資料視覺化。雖然當前Matlab的地位仍然難以撼動，但是，隨著時間的推移，Python在科學計算上的生態系統也會越來越豐富。
安裝Numpy的方法也很簡單，使用Python的包管理工具pip或者anaconda便可以實現。例如在shell中輸入下列命令列便可以透過pip安裝Numpy：
pip install numpy
另外，Numpy是OpenCV的一個依賴庫，所以，我們使用pip工具安裝好OpenCV庫之後，Numpy一般也都已經安裝好了。
01 array型別
Numpy的array型別是該庫的一個基本資料型別，這個資料型別從字面上看是陣列的意思，也就意味著它最關鍵的屬性是元素與維度，我們可以這個資料型別來實現多維陣列。
因此，透過這個資料型別，我們可以使用一維陣列用來表示向量，二維陣列來表示矩陣，以此類推用以表示更高維度的張量。
我們透過下麵的例子來簡單體會一下在Numpy中array型別的使用。
1. Numpy中array型別的基本使用
import numpy as np
# 引入numpy庫，並將其別名為np
array = np.array([1,2,3,4])
# 透過np.array()方法建立一個名為array的array型別，引數是一個list
array
# 在shell中輸入，傳回的結果為：
# array([1, 2, 3, 4])
array.max()
# 獲取該array型別中最大的元素值，結果為:
# 4
array.mean()
# 求該array中元素的平均值，結果為:
# 2.5
array.min()
# 獲取該array中元素的最小值：
# 1
array * 2
# 直接將該array乘以2，Python將該運運算元多載，將每一個元素都乘以了2，其輸出結果為：
# array([2, 4, 6, 8])
array + 1
# 將每一個元素都加上1，輸出結果為：
# array([2, 3, 4, 5])
array / 2
# 將每一個元素都除以2，得到浮點數表示的結果為：
# array([ 0.5, 1. , 1.5, 2. ])
array % 2
# Numpy庫除了可以對array實現除法運算，還可以實現取模運算，結果為：
# array([1, 0, 1, 0], dtype=int32)
array.argmax()
# 獲取該組資料中元素值最大的一個資料的索引，下標從0開始，其結果為：
# 3
透過上面的程式碼片段，我們可以瞭解Numpy中array型別的基本使用方法。我們可以看到，array其實是一個類，透過傳入一個list引數來實體化為一個物件，也就實現了對資料的封裝。這個物件中包含對各個元素進行計算的基本方法，例如求平均值、求最大值等。除此之外，我們再看一下關於更高維度資料的處理。
2. Numpy對更高維度資料的處理
import numpy as np
array = np.array([[1,2],[3,4],[5,6]])
# 建立一個二維陣列，用以表示一個3行2列的矩陣，名為array
array
#在互動式程式設計介面中輸入array，傳回結果為：
# array([[1, 2],
# [3, 4],
# [5, 6]])
array.shape
# 檢視資料的維度屬性，下麵的輸出結果元組，代表的是3行2列
# (3, 2)
array.size
# 檢視array中的元素數量，輸出結果為：
# 6
array.argmax()
# 檢視元素值最大的元素索引，結果為：
# 5
array.flatten()
# 將shape為(3,2)的array轉換為一行表示，輸出結果為：
# array([1, 2, 3, 4, 5, 6])
# 我們可以看到，flatten()方法是將多維資料“壓平”為一維陣列的過程
array.reshape(2,3)
# 將array資料從shape為(3,2)的形式轉換為(2,3)的形式：
# array([[1, 2, 3],
# [4, 5, 6]])
除此之外，Numpy還包含了建立特殊類別的array型別的方法，例如。
3. Numpy建立特殊類別的array型別
import numpy as np
array_zeros = np.zeros((2,3,3))
#生成結果為：
# array([[[ 0., 0., 0.],
# [ 0., 0., 0.],
# [ 0., 0., 0.]],
#
# [[ 0., 0., 0.],
# [ 0., 0., 0.],
# [ 0., 0., 0.]]])
array_ones = np.ones((2,3,3))
# 生成所有元素都為1的array，其shape是(2,3,3)
array_ones.shape
# (2, 3, 3)
array_arange = np.arange(10)
# 生成一個array，從0遞增到10，步長為1，結果為：
# array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
array_linespace = np.linspace(0,10,5)
# 生成一個array從0到10遞增，步長為5，結果為：
# array([ 0. , 2.5, 5. , 7.5, 10. ])
Numpy作為Python的一款著名數值計算庫，其在基礎計算上的功能也是非常完備的，程式碼如下。
4. Numpy基礎計算演示
import numpy as np
np.abs([1,-2,-3,4])
# 取絕對值，結果為：array([1, 2, 3, 4])
np.sin(np.pi / 2)
# 求餘弦值，結果為：1.0
np.arctan(1)
# 求反正切值，結果為：0.78539816339744828
np.exp(2)
# 求自然常數e的2次方，結果為：7.3890560989306504
np.power(2,3)
# 求2的3次方，結果為：8
np.dot([1,2],[3,4])
# 將向量[1,2]與[3,4]求點積，結果為：11
np.sqrt(4)
# 將4開平方，結果為：2.0
np.sum([1,2,3,4])
# 求和，結果為：10
np.mean([1,2,3,4])
# 求平均值，結果為：2.5
np.std([1,2,3,4])
# 求標準差，結果為：1.1180339887498949
除此之外，Numpy所包含的基本計算功能還有很多，例如將array切分、拼接、倒序等。
02 線性代數相關
我們在前面介紹了array型別及其基本操作方法，瞭解到使用array型別可以表示向量、矩陣和多維張量。線性代數計算在科學計算領域非常重要，在機器學習和資料挖掘領域，線性代數相關函式的使用也是非常頻繁的。下麵，我們介紹一下Numpy為我們提供的線性代數操作。
5. Numpy提供的線性代數操作
import numpy as np
vector_a = np.array([1,2,3])
vector_b = np.array([2,3,4])
# 定義兩個向量vector_a與vector_b
np.dot(vector_a,vector_b)
# 將兩個向量相乘，在這裡也就是點乘，結果為20
vector_a.dot(vector_b)
# 將vector_a與vector_b相乘，結果為20
np.dot(vector_a,vector_b.T)
'''
將一個行向量與一個列向量叉乘的結果相當於將兩個行向量求點積，在這裡我們測試了dot()方法。其中array型別的T()方法表示轉置。
測試結果表明：
dot()方法對於兩個向量預設求其點積。對於符合叉乘格式的矩陣，自動進行叉乘。
我們看一下下麵這個例子：
'''
matrix_a = np.array([[1,2],
[3,4]])
# 定義一個2行2列的方陣
matrix_b = np.dot(matrix_a,matrix_a.T)
# 這裡將該方陣與其轉置叉乘，將結果賦予matrix_b變數
matrix_b
'''
結果為：
array([[ 5, 11],
[11, 25]])
'''
np.linalg.norm([1,2])
# 求一個向量的範數的值，結果為2.2360679774997898
# 如果norm()方法沒有指定第二個引數，則預設為求2範數。
np.linalg.norm([1,-2],1)
# 指定第二個引數值為1，即求1範數，我們在前面介紹過，1範數的結果即為向量中各元素絕對值之和，結果為3.0
np.linalg.norm([1,2,3,4],np.inf)
# 求向量的無窮範數，其中np.inf表示正無窮，也就是向量中元素值最大的那個，其結果為4.0
np.linalg.norm([1,2,3,4],-np.inf)
# 同理，求負無窮範數的結果為1，也就是向量中元素的最小值
np.linalg.norm(matrix_b)
# 除了向量可以求範數，矩陣也可以有類似的運算，即為F範數，結果為29.866369046136157
np.linalg.det(matrix_a)
# 求矩陣matrix_a的行列式，結果為-2.0000000000000004
np.trace(matrix_a)
# 求矩陣matrix_a的跡，結果為5
np.linalg.matrix_rank(matrix_a)
# 求矩陣的秩，結果為2
vector_a * vector_b
# 使用*符號將兩個向量相乘，是將兩個向量中的元素分別相乘，也就是前面我們所講到的哈達馬乘積，結果為array([ 2, 6, 12])
vector_a ** vector_b
# 使用二元運運算元**對兩個向量進行操作，結果為array([ 1, 8, 81], dtype=int32)
# 表示將向量vector_a中元素對應vector_b中的元素值求冪運算。例如最終結果[1,8,81]可以表示為：
# [1*1,2*2*2,3*3*3*3]
np.linalg.pinv(matrix_a)
'''
求矩陣的逆矩陣，方法pinv()求的是偽逆矩陣，結果為：
array([[-2. , 1. ],
[ 1.5, -0.5]])
不使用偽逆矩陣的演演算法，直接使用逆矩陣的方法是inv()，即
np.linalg.inv(matrix_a)
結果相同，也為：
array([[-2. , 1. ],
[ 1.5, -0.5]])
'''
03 矩陣的高階函式
我們在前面學習了Numpy的基本資料型別array，同時瞭解了一些基本的數學運算方法。其實除了前面我們所提到的對矩陣求逆、求秩、求轉置等基本運算之外，Numpy還為我們提供了矩陣的分解等更高階的函式。
-
矩陣分解
矩陣分解(decomposition, factorization)是將矩陣拆解為若干個矩陣的相乘的過程。在數值分析，常常被用來實現一些矩陣運算的快速演演算法，在機器學習領域有非常重要的作用。
例如我們在前面介紹過線性降維的PCA演演算法，其中就涉及矩陣分解的步驟。今日頭條、亞馬遜網上商城這類網際網路產品，總會根據我們的個人喜好給我們推送一些它認為我們會感興趣的資訊或商品，這類用於推送訊息的系統稱為推薦系統(Recommendation System)。
在推薦系統的實現過程中，就用到了矩陣分解演演算法。例如主流的開源大資料計算引擎Spark在ml機器學習庫中透過ALS演演算法實現了推薦系統，也有的推薦系統採用SVD演演算法來實現整套系統中的矩陣分解過程。
在Numpy中，為我們提供了基於SVD演演算法的矩陣分解，SVD演演算法即為奇異值分解法，相對於矩陣的特徵值分解法，它可以對非方陣形式的矩陣進行分解，將一個矩陣A分解為如下形式：
A = U∑VT
式中，A代表需要被分解的矩陣，設其維度是m×n。U矩陣是被分解為的三個矩陣之一，它是一個m×m的方陣，構成這個矩陣的向量是正交的，被稱為左奇異向量；∑是一個m×n的向量，它的特點是除了對角線中的元素外，其餘元素都為0。V是一個n×n的方陣，它的轉置也是一個方陣，與U矩陣類似，構成這個矩陣的向量也是正交的，被稱為右奇異向量。整個奇異值分解演演算法矩陣的形式如圖4-1所示，具體演演算法實現在此不再贅述。
▲圖4-1 SVD演演算法的矩陣形式
我們使用Numpy演示一下SVD演演算法的使用。
6. SVD演演算法演示
import numpy as np
matrix = np.array([
[1,2],
[3,4]])
another_matrix = np.dot(matrix,matrix.T)
# 生成一個矩陣 another_matrix
print(another_matrix)
'''
該矩陣為：
array([[ 5, 11],
[11, 25]])
'''
U,s,V = np.linalg.svd(another_matrix,2)
# 使用奇異值分解法將該矩陣進行分解，分解得到三個子矩陣U,s,V
# 在s矩陣的基礎上，生成S矩陣為：
S = np.array([[s[0],0],
[0,s[1]]])
# 我們在下麵看一下生成的幾個矩陣的樣子
print(U)
'''
[[-0.40455358 -0.9145143 ]
[-0.9145143 0.40455358]]
'''
print(s)
'''
[ 29.86606875 0.13393125]
'''
print(V)
'''
[[-0.40455358 -0.9145143 ]
[-0.9145143 0.40455358]]
'''
# 利用生成的U,S,V三個矩陣，我們可以重建回原來的矩陣another_matrix
np.dot(U,np.dot(S,V))
# 輸出結果為:
'''
array([[ 5., 11.],
[ 11., 25.]])
'''
在上面的程式碼片段中，s向量表示的是分解後的∑矩陣中對角線上的元素，所以我們在這裡面引入了一個S矩陣，將s向量中的元素放置在這個矩陣中，用以驗證分解後的矩陣重建回原先的矩陣A的過程。
仔細的讀者可能會註意到，為什麼這裡使用SVD演演算法生成的矩陣U與VT是相同的。大家可能會註意到在上面的程式碼片段中，為何多了一個生成矩陣another_matrix的過程。這是因為一個矩陣與其轉置相乘之後的矩陣是對稱矩陣（矩陣中的元素沿著對角線對稱），將對稱矩陣進行分解後的結果可以表示為：
A = V∑VT
透過觀察上式，我們不難發現U與V矩陣是相同的，因為這個例子中，U與V矩陣本身也是對稱矩陣，不論它的轉置與否形式都是一樣的。
我們在第2章介紹過用於線性降維的PCA演演算法，該演演算法中有一個步驟是將協方差矩陣分解然後重建，下麵我們演示一下使用Numpy的SVD演演算法來實現PCA演演算法的例子：
7. 基於SVD實現PCA演演算法
import numpy as np
# 零均值化，即中心化，是資料的預處理方法
def zero_centered(data):
matrix_mean = np.mean(data, axis=0)
return data - matrix_mean
def pca_eig(data, n):
new_data = zero_centered(data)
cov_mat = np.dot(new_data.T, new_data) # 也可以用 np.cov() 方法
eig_values, eig_vectors = np.linalg.eig(np.mat(cov_mat)) # 求特徵值和特徵向量，特徵向量是列向量
value_indices = np.argsort(eig_values) # 將特徵值從小到大排序
n_vectors = eig_vectors[:, value_indices[-1:-(n + 1):-1]] # 最大的n個特徵值對應的特徵向量
return new_data * n_vectors # 傳回低維特徵空間的資料
def pca_svd(data, n):
new_data = zero_centered(data)
cov_mat = np.dot(new_data.T, new_data)
U, s, V = np.linalg.svd(cov_mat) # 將協方差矩陣奇異值分解
pc = np.dot(new_data, U) # 傳回矩陣的第一個列向量即是降維後的結果
return pc[:, 0]
def unit_test():
data = np.array(
[[2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2], [3.1, 3.0], [2.3, 2.7], [2, 1.6], [1, 1.1], [1.5, 1.6],
[1.1, 0.9]])
result_eig = pca_eig(data, 1) # 使用常規的特徵值分解法，將2維資料降到1維
print(result_eig)
result_svd = pca_svd(data, 1) # 使用奇異值分解法將協方差矩陣分解，得到降維結果
print(result_svd)
if __name__ == '__main__':
unit_test()
經過降維的資料為：
[-0.82797019 1.77758033 -0.99219749 -0.27421042 -1.67580142 -0.9129491
0.09910944 1.14457216 0.43804614 1.22382056]
我們可以看到，資料已經從2維的變為1維的了，這兩個PCA演演算法的計算結果是相同的。其中pca_eig() 函式是使用常規的特徵值分解方法來求解的，讀者可以參照前面講述的PCA演演算法過程來理解這段程式碼。pca_svd() 函式使用奇異值分解法來求解的。這段程式碼雖然相對精簡，但是背後是經過複雜的數學推導的，下麵簡要闡述一下PCA演演算法中奇異值分解的步驟。
1) PCA演演算法中得到樣本的協方差矩陣是經過零均值化處理的，將其去掉常數部分，則也可表示為：
C = XTX
其中，X是經過中心化處理後的樣本矩陣X. 前面我們介紹過，一個矩陣與其轉置矩陣相乘的結果是一個對稱矩陣。觀察到協方差矩陣C便是一個對稱矩陣，那麼將其進行奇異值分解後則可以表示為：
C = V∑VT
2) 將經過中心化的樣本矩陣X進行奇異值分解，可以得到：
X = U∑VT
因此，我們可以得到：
XTX
= (U∑VT)T(U∑VT)
= V∑TUTU∑VT
= V∑2VT
奇異矩陣V中的列對應著PCA演演算法主成分中的主方向，因此可以得到主成分為：
XV = U∑VTV = U∑
關於更詳細的數學推倒過程，讀者可參考該網址：
https://stats.stackexchange.com/questions/134282/relationship-between-svd-and-pca-how-to-use-svd-to-perform-pca
-
隨機數矩陣
Numpy除了為我們提供常規的數學計算函式和矩陣相關操作之外，還提供了很多功能豐富的模組，隨機數模組就是其中一部分。利用隨機數模組可以生成隨機數矩陣，比Python自帶的隨機數模組功能要強大，我們看一下下麵這個例子。
8. Numpy的隨機數功能演示
import numpy as np
# 置隨機數種子：
np.random.seed()
# 從[1,3)中生成一個整數的隨機數，連續生成10個
np.random.randint(1,3,10)
# 傳回：array([2, 2, 1, 1, 1, 1, 2, 1, 2, 2])
# 若要連續產生[1,3)之間的浮點數，可以使用下述方法：
2*np.random.random(10) + 1
'''
傳回：
array([ 1.25705585, 2.38059578, 1.73232769, 2.12303283, 2.33946996,
2.28020734, 2.15724069, 1.32845829, 2.91361293, 1.78637408])
'''
np.random.uniform(1,3,10)
'''
傳回：
array([ 1.37993226, 1.38412227, 1.18063785, 1.75985962, 1.42775752,
1.62100074, 1.71768721, 1.50131522, 2.20297121, 1.08585819])
'''
# 生成一個滿足正太分佈(高斯分佈)的矩陣，其維度是4*4
np.random.normal(size=(4,4))
'''
傳回：
array([[-1.81525915, -2.02236963, 0.90969106, 0.25448426],
[-1.04177298, -0.35408201, 1.67850233, -0.70361323],
[-0.30710761, 0.57461312, -0.37867596, -0.74010685],
[-0.94046747, 2.37124816, -0.78503777, -0.33485225]])
'''
# 隨機產生10個，n=5,p=0.5的二項分佈資料:
np.random.binomial(n=5,p=0.5,size=10)
# 傳回：array([2, 0, 1, 3, 3, 1, 3, 3, 4, 2])
data = np.arange(10) # 產生一個0到9的序列
np.random.choice(data,5) # 從data資料中隨機採集5個樣本，採集過程是有放回的
# 傳回：array([0, 0, 1, 6, 2])
np.random.choice(data,5,replace=False)
#從data資料中隨機採集5個樣本，採集過程是沒有放回的
# 傳回：array([0, 4, 3, 9, 7])
np.random.permutation(data) # 對data進行亂序，傳回亂序結果
# 傳回：array([2, 8, 6, 4, 9, 1, 3, 5, 7, 0])
np.random.shuffle(data) # 對data進行亂序，並替換為新的data
print(data)
# 輸出：[1 2 8 4 3 6 9 0 5 7]
關於作者：王天慶，長期從事分散式系統、資料科學與工程、人工智慧等方面的研究與開發，在人臉識別方面有豐富的實踐經驗。現就職某世界100強企業的資料實驗室，從事資料科學相關技術領域的預研工作。
本文摘編自《Python人臉識別：從入門到工程實踐》，經出版方授權釋出。
延伸閱讀《Python人臉識別》
點選上圖瞭解及購買
轉載請聯絡微信：DoctorData
推薦語：華為資深AI工程師撰寫，全面講解人臉識別各項基礎技術、原理和演演算法，從零實現工程級人臉識別引擎。
朋友會在“發現-看一看”看到你“在看”的內容