Pytorch和Tensorflow在相同資料規模規模下的降維PCA(Principal Component Analysis)演算法中的運算速度對比
技術標籤：tensorflowpytorch機器學習深度學習
Pytorch和Tensorflow在相同資料規模規模下的降維PCA(Principal Component Analysis)演算法中的運算速度對比
在網上找了好久都沒找到pytorch使用特徵值與特徵向量實現pca演算法的教程（即非svd演算法），本文中的程式碼是根據個人對pca的理解所寫的，致敬開源思想，將學習過程與程式放在這裡，希望對大家有幫助，共同學習，共同進步(不喜勿噴)
PCA演算法基本原理
PCA的基本原理是從大量的數中找到少數主要的成分點，在資料維度降低的情況下儘可能地儲存原始資料的資訊。舉個例子：假設有一個二維的資料，其分佈如下面第一個圖所示,藍色粗箭頭的方向就是第一個軸的方向（即離散程度最高大 方差最大），這就意味著資料點在藍色粗線上的投影代表了原始資料的絕大部分資訊。所以我們將藍色的兩條軸作為新的軸（兩條藍軸正交），然後將所有資料點沿著第二軸往第一軸投影，得到第三張圖的資料分佈。
那麼我們要如何找到我們想要的軸呢？（就是離散程度最大的軸）
協方差
我們知道我們知道數值的分散程度，可以用數學上的協方差來表述。協方差可以表示兩個變數之間的相關性。如果協方差為0那麼就說明這兩個變數之間是不相關的。否側他們相關。協方差公式為：
這時我們讓均值為 0，那麼就變為：
為了讓協方差為 0，我們選擇第二個基時只能在與第一個基正交的方向上進行選擇，因此最終選擇的兩個方向一定是正交的。
X
=
(
a
1
a
2
…
…
a
i
b
1
b
2
…
…
b
i
)
X=\begin{pmatrix} a_1 & a_2 & …… & a_i\\\\ b_1 & b_2 & …… & b_i\\\\ \end{pmatrix}
那麼：
1
m
−
1
X
X
T
=
(
1
m
−
1
∑
i
=
1
m
a
i
2
1
m
−
1
∑
i
=
1
m
a
i
b
i
1
m
−
1
∑
i
=
1
m
a
i
b
i
1
m
−
1
∑
i
=
1
m
b
i
2
)
\frac{1}{m-1}XX^T=\begin{pmatrix} \frac{1}{m-1}\sum_{i=1}^{m}a_{i}^2 & \frac{1}{m-1}\sum_{i=1}^{m}a_ib_{i} \\\\ \frac{1}{m-1}\sum_{i=1}^{m}a_ib_{i} & \frac{1}{m-1}\sum_{i=1}^{m}b_{i}^2 \\\\ \end{pmatrix}
=
(
C
o
v
(
a
,
a
)
C
o
v
(
a
,
b
)
C
o
v
(
a
,
b
)
C
o
v
(
b
,
b
)
)
=\begin{pmatrix} Cov(a,a) & Cov(a,b)\\\\ Cov(a,b) &Cov(b,b) \\\\ \end{pmatrix}
=⎝⎜⎜⎛Cov(a,a)Cov(a,b)Cov(a,b)Cov(b,b)⎠⎟⎟⎞
此時我們要主對角線以外的數都為0那麼我們假設
C
=
1
m
−
1
X
X
T
C= \frac{1}{m-1}XX^T
C=m−11XXT而設最終結果的矩陣為P那麼設我們要的變換後的結果為Y=PX，D為Y的協方差矩陣
D
=
1
m
−
1
Y
Y
T
=
1
m
−
1
P
X
(
P
X
)
T
D = \frac{1}{m-1} YY^T = \frac{1}{m-1}PX(PX)^T
D=m−11YYT=m−11PX(PX)T
=
1
m
−
1
P
X
X
T
P
T
=
P
(
1
m
−
1
X
X
T
)
P
T
=
P
C
P
T
= \frac{1}{m-1} PXX^TP^T = P(\frac{1}{m-1}XX^T)P^T =PCP^T
=m−11PXXTPT=P(m−11XXT)PT=PCPT
那麼我們只要求出C的特徵向量和特徵值就可以了，我們所需要的是最大的k（k為降維後的維度）個特徵值對應的特徵向量。
程式碼實現(使用的CPU為銳龍R5-3500U)
pytorch實現
iris資料集下
import torch
import numpy as np
import time
from sklearn import datasets
import matplotlib.pyplot as plt
data = datasets.load_iris(return_X_y=False)
def pca(x,dim=2):
m,n=x.shape[0],x.shape[1]
#print(m,n)
X=torch.from_numpy(x)
#去中心化 防止資料過分逼近大的值 而忽略小的值
X_mean = torch.mean(X,1).reshape(-1,1)
X=X-X_mean.expand_as(X)
# 無偏差的協方差矩陣
cov = torch.matmul(X.T,X)/(m - 1)
#print(cov)
# 計算特徵分解
e,v = torch.eig(cov,eigenvectors=True)
#print(e)
#print(v.shape[0])
e=torch.mean(e,1)*2
#print(e)
# 將特徵值從大到小排序，選出前dim個的index
sorted, e_index_sort = torch.sort(e,descending=True)
e_index_sort=torch.gather(e_index_sort,-1,torch.LongTensor([0,1]))
v_new = torch.index_select(v, 0, e_index_sort)
#print(v_new)
# 降維操作
pca = torch.matmul(X,v_new.T)
return(pca)
#執行程式碼
start = time.perf_counter()
pca_data = pca(data.data,dim=2)
elapsed = (time.perf_counter()-start)
print(f'time use :{elapsed}')
自定義10000*1000資料規模下
import torch
import numpy as np
import time
from sklearn import datasets
import matplotlib.pyplot as plt
n0=2*np.random.randn(5000,1000)+1
n1=-2*np.random.randn(5000,1000)-1
n=np.vstack((n0,n1))
def pca(x,dim=2):
m,n=x.shape[0],x.shape[1]
#print(m,n)
X=torch.from_numpy(x)
#去中心化 防止資料過分逼近大的值 而忽略小的值
X_mean = torch.mean(X,1).reshape(-1,1)
X=X-X_mean.expand_as(X)
# 無偏差的協方差矩陣
cov = torch.matmul(X.T,X)/(m - 1)
#print(cov)
# 計算特徵分解
e,v = torch.eig(cov,eigenvectors=True)
#print(e)
#print(v.shape[0])
e=torch.mean(e,1)*2
#print(e)
# 將特徵值從大到小排序，選出前dim個的index
sorted, e_index_sort = torch.sort(e,descending=True)
e_index_sort=torch.gather(e_index_sort,-1,torch.LongTensor([0,1]))
v_new = torch.index_select(v, 0, e_index_sort)
#print(v_new)
# 降維操作
pca = torch.matmul(X,v_new.T)
return(pca)
#執行程式碼
start = time.perf_counter()
pca_data = pca(data.data,dim=2)
elapsed = (time.perf_counter()-start)
print(f'time use :{elapsed}')
tensorflow實現
iris資料集下
import tensorflow.compat.v1 as tf
# 使用Eager Execution動態圖機制
tf.enable_eager_execution()
import numpy as np
from sklearn import datasets
import seaborn as sns
import matplotlib.pyplot as plt
import time
data = datasets.load_iris(return_X_y=False)
def pca(x,dim = 2):
'''
x:輸入矩陣
dim:降維之後的維度數
'''
with tf.name_scope("PCA"):#建立一個引數名空間
m,n= tf.to_float(x.get_shape()[0]),tf.to_int32(x.get_shape()[1])
assert not tf.assert_less(dim,n)
mean = tf.reduce_mean(x,axis=1)
# 去中心化 防止資料過分逼近大的值 而忽略小的值
x_new = x - tf.reshape(mean,(-1,1))
# 無偏差的協方差矩陣
cov = tf.matmul(x_new,x_new,transpose_a=True)/(m - 1)
#print(cov)
# 計算特徵分解
e,v = tf.linalg.eigh(cov,name="eigh")
#print(e)
#print(v)
# 將特徵值從大到小排序，選出前dim個的index
e_index_sort = tf.math.top_k(e,sorted=True,k=dim)[1]
#print(e_index_sort)
# 提取前排序後dim個特徵向量
v_new = tf.gather(v,indices=e_index_sort)
#print(v_new)
# 降維操作
pca = tf.matmul(x_new,v_new,transpose_b=True)
return pca
start = time.perf_counter()
pca_data = tf.constant(np.reshape(data.data,(data.data.shape[0],-1)),dtype=tf.float32)
pca_data = pca(pca_data,dim=2)
elapsed = (time.perf_counter()-start)
print(f'time use :{elapsed}')
自定義10000*1000資料規模下
import tensorflow.compat.v1 as tf
# 使用Eager Execution動態圖機制
tf.enable_eager_execution()
import numpy as np
from sklearn import datasets
import seaborn as sns
import matplotlib.pyplot as plt
import time
n0=2*np.random.randn(5000,1000)+1
n1=-2*np.random.randn(5000,1000)-1
n=np.vstack((n0,n1))
def pca(x,dim = 2):
'''
x:輸入矩陣
dim:降維之後的維度數
'''
with tf.name_scope("PCA"):#建立一個引數名空間
m,n= tf.to_float(x.get_shape()[0]),tf.to_int32(x.get_shape()[1])
assert not tf.assert_less(dim,n)
mean = tf.reduce_mean(x,axis=1)
# 去中心化 防止資料過分逼近大的值 而忽略小的值
x_new = x - tf.reshape(mean,(-1,1))
# 無偏差的協方差矩陣
cov = tf.matmul(x_new,x_new,transpose_a=True)/(m - 1)
#print(cov)
# 計算特徵分解
e,v = tf.linalg.eigh(cov,name="eigh")
#print(e)
#print(v)
# 將特徵值從大到小排序，選出前dim個的index
e_index_sort = tf.math.top_k(e,sorted=True,k=dim)[1]
#print(e_index_sort)
# 提取前排序後dim個特徵向量
v_new = tf.gather(v,indices=e_index_sort)
#print(v_new)
# 降維操作
pca = tf.matmul(x_new,v_new,transpose_b=True)
return pca
start = time.perf_counter()
pca_data = tf.constant(np.reshape(data.data,(data.data.shape[0],-1)),dtype=tf.float32)
pca_data = pca(pca_data,dim=2)
elapsed = (time.perf_counter()-start)
print(f'time use :{elapsed}')
因為pca演算法在計算的過程中運算量並不大，所以兩個架構在執行該演算法的時候時間相差並不大，但是也可以看出在不同的資料集下，兩個框架運算速度不同，說明資料對框架運算速度起到很大的影響（以上為個人觀點，不喜勿噴）