從零開始實現主成分分析(PCA)演算法

來自專欄 手把手教你實現機器學習演算法

聲明：版權所有，轉載請聯繫作者並註明出處: 風雪夜歸子 - CSDN博客

知乎專欄： https://www.zhihu.com/people/feng-xue-ye-gui-zi

由於知乎不支持markdown格式，所以有公式的地方都使用了截圖，若影響閱讀效果，可移步我的Blog：風雪夜歸子 - CSDN博客，文章會同步更新！

前面兩篇文章詳細講解了線性判別分析LDA，說到LDA，就不能不提到主成份分析，簡稱為PCA，是一種非監督學習演算法，經常被用來進行數據降維、有損數據壓縮、特徵抽取、數據可視化（Jolliffe, 2002）。它也被稱為Karhunen-Loève變換。

1. PCA原理

PCA的思想是將$n$維特徵映射到$k$維空間上$k<n$，這$k$維特徵是全新的正交特徵，是重新構造出來的$k$維特徵，而不是簡單地從$n$維特徵中去除其餘$n-k$維特徵。那麼如何衡量投影向量的優劣呢？在數學上有三種方法衡量投影的優劣！PCA可以被定義為數據在低維線性空間上的正交投影，這個線性空間被稱為主?空間（principal subspace），使得投影數據的?差被最?化（Hotelling, 1933），即最大方差理論。等價地，它也可以被定義為使得平均投影代價最?的線性投影，即最小誤差理論。平均投影代價是指數據點和它們的投影之間的平均平?距離（Pearson, 1901）。還有另一個理論也可以解釋PCA原理，即坐標軸相關度理論。這裡簡單探討前兩種，最後一種在討論PCA意義時簡單概述。

3. 代碼實現

from __future__ import print_functionfrom sklearn import datasetsimport matplotlib.pyplot as pltimport matplotlib.cm as cmximport matplotlib.colors as colorsimport numpy as np%matplotlib inlinedef shuffle_data(X, y, seed=None): if seed: np.random.seed(seed) idx = np.arange(X.shape[0]) np.random.shuffle(idx) return X[idx], y[idx]# 正規化數據集 Xdef normalize(X, axis=-1, p=2): lp_norm = np.atleast_1d(np.linalg.norm(X, p, axis)) lp_norm[lp_norm == 0] = 1 return X / np.expand_dims(lp_norm, axis)# 標準化數據集 Xdef standardize(X): X_std = np.zeros(X.shape) mean = X.mean(axis=0) std = X.std(axis=0) # 做除法運算時請永遠記住分母不能等於0的情形 # X_std = (X - X.mean(axis=0)) / X.std(axis=0) for col in range(np.shape(X)[1]): if std[col]: X_std[:, col] = (X_std[:, col] - mean[col]) / std[col] return X_std# 劃分數據集為訓練集和測試集def train_test_split(X, y, test_size=0.2, shuffle=True, seed=None): if shuffle: X, y = shuffle_data(X, y, seed) n_train_samples = int(X.shape[0] * (1-test_size)) x_train, x_test = X[:n_train_samples], X[n_train_samples:] y_train, y_test = y[:n_train_samples], y[n_train_samples:] return x_train, x_test, y_train, y_test# 計算矩陣X的協方差矩陣def calculate_covariance_matrix(X, Y=np.empty((0,0))): if not Y.any(): Y = X n_samples = np.shape(X)[0] covariance_matrix = (1 / (n_samples-1)) * (X - X.mean(axis=0)).T.dot(Y - Y.mean(axis=0)) return np.array(covariance_matrix, dtype=float)# 計算數據集X每列的方差def calculate_variance(X): n_samples = np.shape(X)[0] variance = (1 / n_samples) * np.diag((X - X.mean(axis=0)).T.dot(X - X.mean(axis=0))) return variance# 計算數據集X每列的標準差def calculate_std_dev(X): std_dev = np.sqrt(calculate_variance(X)) return std_dev# 計算相關係數矩陣def calculate_correlation_matrix(X, Y=np.empty([0])): # 先計算協方差矩陣 covariance_matrix = calculate_covariance_matrix(X, Y) # 計算X, Y的標準差 std_dev_X = np.expand_dims(calculate_std_dev(X), 1) std_dev_y = np.expand_dims(calculate_std_dev(Y), 1) correlation_matrix = np.divide(covariance_matrix, std_dev_X.dot(std_dev_y.T)) return np.array(correlation_matrix, dtype=float)class PCA(): """ 主成份分析演算法PCA，非監督學習演算法. """ def __init__(self): self.eigen_values = None self.eigen_vectors = None self.k = 2 def transform(self, X): """ 將原始數據集X通過PCA進行降維 """ covariance = calculate_covariance_matrix(X) # 求解特徵值和特徵向量 self.eigen_values, self.eigen_vectors = np.linalg.eig(covariance) # 將特徵值從大到小進行排序，注意特徵向量是按列排的，即self.eigen_vectors第k列是self.eigen_values中第k個特徵值對應的特徵向量 idx = self.eigen_values.argsort()[::-1] eigenvalues = self.eigen_values[idx][:self.k] eigenvectors = self.eigen_vectors[:, idx][:, :self.k] # 將原始數據集X映射到低維空間 X_transformed = X.dot(eigenvectors) return X_transformeddef main(): # Load the dataset data = datasets.load_iris() X = data.data y = data.target # 將數據集X映射到低維空間 X_trans = PCA().transform(X) x1 = X_trans[:, 0] x2 = X_trans[:, 1] cmap = plt.get_cmap(viridis) colors = [cmap(i) for i in np.linspace(0, 1, len(np.unique(y)))] class_distr = [] # Plot the different class distributions for i, l in enumerate(np.unique(y)): _x1 = x1[y == l] _x2 = x2[y == l] _y = y[y == l] class_distr.append(plt.scatter(_x1, _x2, color=colors[i])) # Add a legend plt.legend(class_distr, y, loc=1) # Axis labels plt.xlabel(Principal Component 1) plt.ylabel(Principal Component 2) plt.show()if __name__ == "__main__": main()

參考文獻：

《模式識別和機器學習》

主成分分析（Principal components analysis）-最小平方誤差解釋