Python機器學習的練習七：K-Means聚類和主成分分析
這部分練習涵蓋兩個吸引人的話題：K-Means聚類和主成分分析（PCA），K-Means和PCA都是無監督學習技術的例子，無監督學習問題沒有為我們提供任何標籤或者目標去學習做出預測，所以無監督演算法試圖從資料本身中學習一些有趣的結構，我們將首先實現k-means，並瞭解如何使用它來壓縮影象。我們還將用PCA進行實驗，以發現面部影象的低維度表示。
K-Means聚類
首先，我們在一個簡單的二維資料集上實現並應用k-means，以瞭解它如何工作。k-means是一種迭代的、無監督的聚類演算法，它將類似的例項組合成叢集。該演算法通過猜測每個叢集的初始centroid，反覆向最近的叢集分配例項，並重新計算該叢集的centroid。首先我們要實現一個函式，它為資料中的每個例項找到最接近的centroid。
import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sb from scipy.ioimport loadmat %matplotlib inline def find_closest_centroids(X, centroids): m= X.shape[0] k= centroids.shape[0] idx= np.zeros(m) for iin range(m): min_dist= 1000000 for jin range(k): dist= np.sum((X[i,:]- centroids[j,:])** 2) if dist < min_dist: min_dist= dist idx[i]= j return idx
測試函式確保它像預期的那樣工作，我們使用練習中的測試案例。
data= loadmat('data/ex7data2.mat')
X= data['X']
initial_centroids= initial_centroids= np.array([[3,3], [6,2], [8,5]])
idx= find_closest_centroids(X, initial_centroids)
idx[0:3]
array([0., 2., 1.])輸出與文字中的預期值相匹配(我們的陣列是zero-indexed而不是one-indexed，所以值比練習中的值要低1)。接下來，我們需要一個函式來計算叢集的centroid。centroid是當前分配給叢集的所有例子的平均值。
def compute_centroids(X, idx, k):
m, n= X.shape
centroids= np.zeros((k, n))
for iin range(k):
indices= np.where(idx== i)
centroids[i,:]= (np.sum(X[indices,:], axis=1)/ len(indices[0])).ravel()
return centroids
compute_centroids(X, idx,3)
array([[2.42830111, 3.15792418],
[5.81350331, 2.63365645],
[7.11938687, 3.6166844 ]])此輸出也與該練習的預期值相匹配。目前為止一切都很順利。下一部分涉及到實際執行演算法的迭代次數和視覺化結果。我們在練習中實現了這一步驟，它沒有那麼複雜，我將從頭開始構建它。為了執行這個演算法，我們只需要在分配到最近叢集的示例和重新計算叢集的centroids之間進行交替操作。
def run_k_means(X, initial_centroids, max_iters):
m, n= X.shape
k= initial_centroids.shape[0]
idx= np.zeros(m)
centroids= initial_centroids
for iin range(max_iters):
idx= find_closest_centroids(X, centroids)
centroids= compute_centroids(X, idx, k)
return idx, centroids
idx, centroids= run_k_means(X, initial_centroids,10)我們現在可以使用顏色編碼表示叢集成員。
cluster1= X[np.where(idx== 0)[0],:]
cluster2= X[np.where(idx== 1)[0],:]
cluster3= X[np.where(idx== 2)[0],:]
fig, ax= plt.subplots(figsize=(12,8))
ax.scatter(cluster1[:,0], cluster1[:,1], s=30, color='r', label='Cluster 1')
ax.scatter(cluster2[:,0], cluster2[:,1], s=30, color='g', label='Cluster 2')
ax.scatter(cluster3[:,0], cluster3[:,1], s=30, color='b', label='Cluster 3')
ax.legend()
我們跳過了初始化centroid的過程。這可能會影響演算法的收斂性。
接下來建立一個可以選擇隨機例子的函式，並將這些例子作為初始的centroid。
def init_centroids(X, k):
m, n= X.shape
centroids= np.zeros((k, n))
idx= np.random.randint(0, m, k)
for iin range(k):
centroids[i,:]= X[idx[i],:]
return centroids
init_centroids(X,3)
array([[1.15354031, 4.67866717],
[6.27376271, 2.24256036],
[2.20960296, 4.91469264]])我們的下一任務是應用K-means實現影象壓縮。我們可以使用叢集來查詢影象中最具有代表性的少量的顏色，並使用叢集分配將原來的24位顏色對映到一個低維度的顏色空間。這是我們要壓縮的影象。
原始畫素資料已經預載入了，把它輸入進來。
image_data= loadmat('data/bird_small.mat')
image_data
{'A': array([[[219,180,103],
[230,185,116],
[226,186,110],
...,
[14, 15, 13],
[13, 15, 12],
[12, 14, 12]],
...,
[[15, 19, 19],
[20, 20, 18],
[18, 19, 17],
...,
[65, 43, 39],
[58, 37, 38],
[52, 39, 34]]], dtype=uint8),
'__globals__': [],
'__header__':'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Tue Jun 5 04:06:24 2012',
'__version__':'1.0'}我們可以快速檢視資料的形狀，以驗證它是否像我們預期的影象。
A= image_data['A'] A.shape(128L,128L,3L)現在我們需要對資料進行預處理，並將它輸入到k-means演算法中。
# normalize value ranges
A= A/ 255.
# reshape the array
X= np.reshape(A, (A.shape[0]* A.shape[1], A.shape[2]))
# randomly initialize the centroids
initial_centroids= init_centroids(X,16)
# run the algorithm
idx, centroids= run_k_means(X, initial_centroids,10)
# get the closest centroids one last time
idx= find_closest_centroids(X, centroids)
# map each pixel to the centroid value
X_recovered= centroids[idx.astype(int),:]
# reshape to the original dimensions
X_recovered= np.reshape(X_recovered, (A.shape[0], A.shape[1], A.shape[2]))
plt.imshow(X_recovered)
我們在壓縮中建立了一些artifact，儘管將原始影象對映到僅16種顏色，但影象的主要特徵仍然存在。
這是關於k-means的部分，接下來我們來看關於主成分分析的部分。
主成分分析
PCA是一個可以在資料集中找到“主成分”或者最大方差方向的線性變換。它可以用於其他事物的維度減少。在這個練習中，我們需要實現PCA，並將其應用於一個簡單的二維資料集，觀察它是如何工作的。從載入和視覺化資料集開始。
data= loadmat('data/ex7data1.mat')
X= data['X']
fig, ax= plt.subplots(figsize=(12,8))
ax.scatter(X[:,0], X[:,1])
PCA的演算法相當簡單。在保證資料正規化後，輸出只是原始資料協方差矩陣的單值分解。由於numpy已經有內建函式來計算矩陣協方差和SVD，我們將利用這些函式而不是從頭開始。
def pca(X):
# normalize the features
X= (X- X.mean())/ X.std()
# compute the covariance matrix
X= np.matrix(X)
cov= (X.T* X)/ X.shape[0]
# perform SVD
U, S, V= np.linalg.svd(cov)
return U, S, V
U, S, V= pca(X)
U, S, V
(matrix([[-0.79241747,-0.60997914],
[-0.60997914, 0.79241747]]),
array([1.43584536, 0.56415464]),
matrix([[-0.79241747,-0.60997914],
[-0.60997914, 0.79241747]]))現在我們已經有了主成分（矩陣U），我們可以利用它把原始資料投入到一個更低維度的空間，對於這個任務，我們將實現一個函式，它計算投影並只選擇頂部K成分，有效地減少了維度的數量。
def project_data(X, U, k):
U_reduced= U[:,:k]
return np.dot(X, U_reduced)
Z= project_data(X, U,1)
Z
matrix([[-4.74689738],
[-7.15889408],
[-4.79563345],
[-4.45754509],
[-4.80263579],
...,
[-6.44590096],
[-2.69118076],
[-4.61386195],
[-5.88236227],
[-7.76732508]])我們也可以通過改變採取的步驟來恢復原始資料。
def recover_data(Z, U, k):
U_reduced= U[:,:k]
return np.dot(Z, U_reduced.T)
X_recovered= recover_data(Z, U,1)
X_recovered
matrix([[3.76152442, 2.89550838],
[5.67283275, 4.36677606],
[3.80014373, 2.92523637],
[3.53223661, 2.71900952],
[3.80569251, 2.92950765],
...,
[5.10784454, 3.93186513],
[2.13253865, 1.64156413],
[3.65610482, 2.81435955],
[4.66128664, 3.58811828],
[6.1549641 , 4.73790627]])如果我們嘗試去視覺化恢復的資料，會很容易的發現演算法的工作原理。
fig, ax= plt.subplots(figsize=(12,8))
ax.scatter(X_recovered[:,0], X_recovered[:,1])
注意這些點如何被壓縮成一條虛線。虛線本質上是第一個主成分。當我們將資料減少到一個維度時，我們切斷的第二個主成分可以被認為是與這條虛線的正交變化。由於我們失去了這些資訊，我們的重建只能將這些點與第一個主成分相關聯。
我們這次練習的最後一項任務是將PCA應用於臉部影象。通過使用相同降維技術，我們可以使用比原始影象少得多的資料來捕捉影象的“本質”。
faces= loadmat('data/ex7faces.mat')
X= faces['X']
X.shape
(5000L,1024L)該練習程式碼包含一個函式，它將在網格中的資料集中渲染前100個臉部影象。你可以在練習文字中找到它們，不需要重新生成。
face= np.reshape(X[3,:], (32,32))
plt.imshow(face)
只有32 x 32灰度影象。下一步我們要在臉部影象資料集上執行PCA，並取得前100個主成分。
U, S, V= pca(X)
Z= project_data(X, U,100)現在嘗試恢復原來的結構並重新渲染它。
X_recovered= recover_data(Z, U,100)
face= np.reshape(X_recovered[3,:], (32,32))
plt.imshow(face)
結果並沒有像預期的維度數量減少10倍，可能是因為我們丟失了一些細節部分。
本文為編譯文章，作者John Wittenauer，原網址
http://www.johnwittenauer.net/machine-learning-exercises-in-python-part-7/