機器學習練習（七）—— K-均值聚類與主成分分析
這篇文章是一系列 Andrew Ng 在 Coursera 上的機器學習課程的練習的一部分。這篇文章的原始程式碼，練習文字，資料檔案可從這裡獲得。
現在我們到了本系列最後兩篇文章了！在本部分，我們將會討論兩個有趣的主題： K-means 聚類和主成分分析（PCA）。K-means 和 PCA 都屬於無監督學習（ unsupervised learning）技術。 無監督學習問題沒有任何標籤或目標讓我們學習，從而不能從中作出預測，因此無監督演算法試圖從資料本身學習到一些感興趣的內容。 我們將首先實現 K-means，並看看如何使用它進行影象壓縮。 我們還將用 PCA 進行實驗，來找到面部影象的低維表示。 一如往常，使用練習文字將有助於我們更好的理解（上傳
K-Means Clustering
首先，我們將實現和應用 K-means 到一個簡單的 2 維資料集，以獲得一些關於它工作原理的直覺。 K-means 是一種迭代的，無監督的聚類演算法，將相似的例項組成簇。 該演算法通過猜測每個聚類的中心開始，然後重複地將例項分配給最近的聚類並重新計算該聚類的中心。 我們首先要實現是一個找出每個例項最近的中心的函式。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sb
from scipy.io import 讓我們來測試一下這個函式以確認它按照我們預期的方式執行。我們將使用本次練習中提供的測試樣例。
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.])該輸出和文字中的期望值相匹配（記住，我們的陣列索引是從 0 開始，而不是從 1 開始， 所以陣列中值比練習中的低一位）。 接下來，我們需要一個函式來計算聚類中心。中心僅僅意味著當前分配給聚類的所有樣本點的均值。
def compute_centroids(X, idx, k):
m, n = X.shape
centroids = np.zeros((k, n))
for i in 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 ]])此輸出也與練習中的預期值相匹配。 到現在為止，一切進展順利。 下一部分就是實際執行演算法來進行幾次迭代和視覺化結果。 這個步驟在練習中已經實現了，但由於它不是那麼複雜，我將從頭開始實現它。 為了執行演算法，我們只需交替進行分配樣本點到最近的聚類和重新計算聚類中心這兩步。
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 i in 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() 
在這裡，我們跳過了初始化聚類中心這一步。 這會影響演算法的收斂。 我們的任務是建立一個選擇隨機樣本點並將它們用作初始中心的函式。
def init_centroids(X, k):
m, n = X.shape
centroids = np.zeros((k, n))
idx = np.random.randint(0, m, k)
for i in 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 位顏色對映到較低維的顏色空間。 下面是我們要壓縮的影象：
原始畫素資料已經預先載入了，讓我們 pull 一下。
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'}
我們可以快速地看一下資料的 shape，以驗證它與我們所期望的影象一致。
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) 
Cool! 你可以看到我們在壓縮中創造了一些工件，但是影象的主要特徵還在，儘管原始影象被對映為了 16 位顏色。這就是 K-means。接下來，我們進行主成分分析（PCA）。
Principal Component Analysis
PCA 是一種線性變換，其目的是找到資料集中的“主分量”或最大方差方向。 除此之外，它還可以用於降維。 在本練習中，我們首先要實現並應用 PCA 於一個簡單的 2 維資料集，以瞭解它是如何工作的。 首先，讓我們載入和視覺化資料集。
data = loadmat('data/ex7data1.mat')
X = data['X']
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(X[:, 0], X[:, 1]) 
PCA 的演算法相當簡單。 在確保資料被歸一化之後，其輸出就是原始資料協方差矩陣的奇異值分解（singular value decomposition）。 由於 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×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 倍。
練習 7 結束！ 在最後一個練習中，我們將實現異常檢測的演算法，並使用協同過濾（collaborative filtering）構建一個推薦系統。