8段代碼演示Numpy數據運算的神操作

作者|王天慶

來源|大數據(ID:hzdashuju)

導讀: 本文介紹一下在Python科學計算中非常重要的一個庫——Numpy。

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強企業的數據實驗室,從事數據科學相關技術領域的預研工作。

本文僅代表作者觀點,轉載請聯系原作者)*



?------------------------------------------------

原文地址:https://mp.weixin.qq.com/s/rXCKCm_P9pNB06CN-kInHQ

轉載請標明來之:阿貓學編程
更多教程:阿貓學編程-python基礎教程

所有評論

如果對文章有異議,請加qq:1752338621