EX7：K-均值聚類和PCA
前言：本練習中，我們將利用K-均值算法壓縮一張圖片，第二部分中，將使用PCA為面部圖片找尋低維描述。
1、K-均值聚類
在第一個練習中，主要實現K-means算法並將其用於圖像壓縮。 首先從2D數據集樣本開始，目的了解K-means算法如何工作的直觀性。 之后，將使用K-means算法進行圖像壓縮，將圖像中發生的顏色數量減少為僅圖像中最常見的顏色。
K-means算法是一種自動將類似數據樣本聚類在一起的方法。 具體來說，給定一個訓練集{\(x^{(1)},x^{(2)},...,x^{(m)}\)}（其中\(x^{(i)} \in \mathbb{R}^n\)），希望將數據分組為幾個集群。K-means的直觀體現是一個迭代過程，它通過猜測初始聚類中心點（centroids）開始，然后重復地將樣本分配給最接近的中心，基於分配來重新計算中心點。
K-均值算法流程如下：
%初始化中心點
centroids = kMeansInitCentroids(x,K);
for iter = 1:iterations
%聚類分配步，將數據集分配到最鄰近的中心點.
%idx(i)對應於c^(i),為樣本點所屬中心點的索引
idx = findClosestCentroids(X,centroids);
%中心點移動步，計算所屬樣本下的均值點
centroids = computeMeans(X,idx,K);
end
算法的內循環重復執行兩個步驟：（1）將每個訓練樣本x（i）分配到其最近的中心點，（2）使用分配給它的點重新計算每個中心點的平均值。注意，融合解決方案可能並不總是理想的，並且取決於質心的初始設置。 因此，實際上，K-means算法通常用不同的隨機初始化運行幾次，從不同的隨機初始化中選擇這些不同解決方案的一種方法是選擇具有最低成本函數值（失真）的解決方案。
1.1找尋最近的中心點
在K-均值算法的聚類分配中，每個訓練樣本\(x^{(i)}\)被分配到最近的中心點，對於每一個樣本i，我們設置\(c^{(i)}\):= j (其中j是 \(||x{(i)}-\mu_j||^2\)取得最小值時)。即\(c^{(i)}\)為最接近點\(x^{(i)}\)處中心點的索引，程序中以idx(i)表述。程序名為：findClosestCentroids.m
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
% idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
% in idx for a dataset X where each row is a single example. idx = m x 1
% vector of centroid assignments (i.e. each entry in range [1..K])
%
% Set K
K = size(centroids, 1);
% You need to return the following variables correctly.
idx = zeros(size(X,1), 1)
% ====================== YOUR CODE HERE ======================
m = size(X,1);
temp = zeros(K, 1);
temp_p = [];
for i = 1:m
for j = 1:K
temp(j) = (X(i,:)-centroids(j,:))*(X(i,:)-centroids(j,:))';
end
temp_p = find(temp==min(temp));
idx(i) = temp_p(1);%避免出現同時兩個最小值
end
end
1.2計算中心點-centroid
對於中心點所屬的點，第二步就是通過對所屬的點計算均值重新定義中心點，對於每個中心點k，我們有：
，這里\(C_k\)為屬於中心點k的樣本集。
如：樣本\(x^{(3)}、x^{(5)}\)都屬於中心點k=2，這時我們應該更新樣本點：\(\mu_2=\frac{1}{2}(x^{(3)}+y^{(5)})\)。
程序名為：computeCentroids.m
function centroids = computeCentroids(X, idx, K)
% Useful variables
[m n] = size(X);
% You need to return the following variables correctly.
centroids = zeros(K, n);
% ====================== YOUR CODE HERE ======================
temp = [];
for i = 1:K
temp = X(find(idx==i),:);
centroids(i,:) = sum(temp)/(size(temp,1));
end
end
1.3樣本數據下的K-均值算法表現
完成上述兩個函數（findClosestCentroids和computeCentroids）后，下一步將在二維數據集上運行K-means算法，以了解K-means如何工作。 函數從runKmeans.m腳本中調用，具體如下：
function [centroids, idx] = runkMeans(X, initial_centroids, ...
max_iters, plot_progress)
% Set default value for plot progress
if ~exist('plot_progress', 'var') || isempty(plot_progress)
plot_progress = false;
end
% Plot the data if we are plotting progress
if plot_progress
figure;
hold on;
end
% Initialize values
[m n] = size(X);
K = size(initial_centroids, 1);
centroids = initial_centroids;
previous_centroids = centroids;
idx = zeros(m, 1);
% Run K-Means
for i=1:max_iters
% Output progress
fprintf('K-Means iteration %d/%d...\n', i, max_iters);
if exist('OCTAVE_VERSION')
fflush(stdout);
end
% 步一
idx = findClosestCentroids(X, centroids);
% Optionally, plot progress here
if plot_progress
plotProgresskMeans(X, centroids, previous_centroids, idx, K, i);
previous_centroids = centroids;
fprintf('Press enter to continue.\n');
pause;
end
% 步二
centroids = computeCentroids(X, idx, K);
end
% Hold off if we are plotting progress
if plot_progress
hold off;
end
end
此外在本腳本中，對k-均值計算的每一步可視化中都調用了函數plotProgressMeans.m，code如下：
function plotProgresskMeans(X, centroids, previous, idx, K, i)
%It is intended for use only with 2D data.
% Plot the examples
plotDataPoints(X, idx, K);
%plotDataPoints如下：
%function plotDataPoints(X, idx, K)
%% Create palette-創建調色板
%palette = hsv(K + 1);
%colors = palette(idx, :);
%
%% Plot the data
%scatter(X(:,1), X(:,2), 15, colors);%散點圖，又稱氣泡圖
% Plot the centroids as black x's
plot(centroids(:,1), centroids(:,2), 'x', ...
'MarkerEdgeColor','k', ...
'MarkerSize', 10, 'LineWidth', 3);
% Plot the history of the centroids with lines
for j=1:size(centroids,1)
drawLine(centroids(j, :), previous(j, :));%plot([p1(1) p2(1)], [p1(2) p2(2)], varargin{:});
end
% Title
title(sprintf('Iteration number %d', i))
end
在運行下一步時，K-means代碼將產生可視化，對每次迭代的進度和進度可視化表現。 按多次輸入以查看K-means算法的每個步驟如何更改質心和集群分配。 設置步長為10步，運行結果如下圖所示：
以上程序中的中心點的初始賦值是給定好的。實際上，初始化中心點的良好策略是從訓練集的樣本中隨機挑選。
函數kMensInitCentroids.m給出挑選過程，具體如下：
% Initialize the centroids to be random examples
% 重新隨機排序樣本的索引
randidx = randperm(size(X, 1));%randperm函數
% Take the first K examples as centroids
centroids = X(randidx(1:K), :);
上面的代碼首先隨機排列了樣本的索引（使用randperm）。 然后，它基於索引的隨機排列來選擇前K個樣本。保證了隨機選擇這些樣本的同時，不會選擇相同的樣本兩次的風險。
1.5基於k-均值的圖片壓縮
在下面的24位真彩色圖中（如，0x000000為黑色，0xFFFFFF為白色），每個像素可表示為指定紅色，綠色和藍色強度值的三個8位無符號整數（范圍從0到255）。 這種編碼通常被稱為RGB編碼。 我們的圖像包含數千種顏色，在這部分練習中，將減少為16種顏色的顏色數量。
通過這種減少，可以有效地表示（壓縮）照片。 具體來說，只需要存儲16個所選顏色的RGB值，並且對於圖像中的每個像素，現在只需要在該位置存儲顏色的索引（只需要4位來表示16種可能性）。
練習中利用K-means算法來選擇將用於表示壓縮圖像的16種顏色。 具體來說，將原始圖像中的每個像素視為數據示例，並使用K-means算法找到最佳組（聚類）三維RGB空間中的像素的16種顏色。 一旦計算了圖像上的聚類中心點，將使用16種顏色來替換原始圖像中的像素。
單個像素下的K-均值運用
首先，在Octave/MATLAB中，圖片可以通過以下的命令進行讀取：
% Load 128x128 color image (bird small.png)
A = imread('bird small.png');
% If you do not have the image package installed,
% you should instead change the following line to
% load('bird small.mat'); % Loads the image into the variable A
上述代碼所創建的一個三維矩陣A，其前兩個索引標識像素位置，其最后一個索引表示紅色，綠色或藍色。 例如，A（50,33,3）給出行50和列33處的像素的藍色強度。
ex7.m中的代碼首先加載圖像，然后重新構圖，創建一個m×3像素顏色矩陣（其中m = 16384 = 128×128），並在其上調用K-means函數。
在尋找到前K=16個顏色代表圖片后，現在可以使用函數findClosetCentroids將每個像素位置最近的中心點做替換。即使用每個像素的中心點來表示原始圖像。
請注意，此時已經顯着減少了描述圖像所需的位數。 原始圖像對於128×128像素位置中的每一個像素的描述需要24位，總大小為128×128×24 = 393216比特。 新的表示形式需要以字典形式的形式存儲16種顏色，每種顏色需要24位，但是圖像本身只需要每像素位置4位。 因此，使用的最終位數為16×24 + 128×128×4 = 65920比特，其對應於將原始圖像壓縮約6倍。運行程序，我們可以看到以下圖像，左側為原圖像，右圖為壓縮后的圖像。
代碼如下：
% Load an image of a bird
A = double(imread('bird_small.png'));
A = A / 255; % Divide by 255 so that all values are in the range 0 - 1
% Size of the image
img_size = size(A);
% Reshape the image into an Nx3 matrix where N = number of pixels.
% Each row will contain the Red, Green and Blue pixel values
% This gives us our dataset matrix X that we will use K-Means on.
X = reshape(A, img_size(1) * img_size(2), 3);%學習之處
% Run your K-Means algorithm on this data
% You should try different values of K and max_iters here
K = 16;
max_iters = 10;
% When using K-Means, it is important the initialize the centroids
% randomly.
% You should complete the code in kMeansInitCentroids.m before proceeding
initial_centroids = kMeansInitCentroids(X, K);
% Run K-Means
[centroids, idx] = runkMeans(X, initial_centroids, max_iters);
% Find closest cluster members
idx = findClosestCentroids(X, centroids);
% Essentially, now we have represented the image X as in terms of the
% indices in idx.
% We can now recover the image from the indices (idx) by mapping each pixel
% (specified by its index in idx) to the centroid value
X_recovered = centroids(idx,:);%直接寫：centroids(idx,:)
% Reshape the recovered image into proper dimensions
X_recovered = reshape(X_recovered, img_size(1), img_size(2), 3);%值得學習
% Display the original image
subplot(1, 2, 1);
imagesc(A);
title('Original');
%imagesc(A) 將矩陣A中的元素數值按大小轉化為不同顏色，並在坐標軸對應位置處以這種顏色染色
% Display compressed image side by side
subplot(1, 2, 2);
imagesc(X_recovered)
title(sprintf('Compressed, with %d colors.', K));
2.PCA-中心成分分析法
本節練習中，將使用PCA來完成維度縮減，將首先從二維的數據集開始，從而對PCA怎樣工作形成直觀的認識，接着對多達5000份的面部圖像進行處理。
2.1數據集樣本
為了從直觀上了解PCA的工作原理，首先從2D數據集開始，該數據集具有一個大的變化方向和一個較小的變化。 如下圖為腳本ex7 pca.m所繪制訓練數據。 在這部分練習中，將可以看出使用PCA將數據從二維縮小到一維時會發生什么。
2.2實施PCA
PCA由兩個計算步驟組成：首先，計算數據的協方差矩陣，然后使用Octave / MATLAB的SVD函數計算特征向量U1，U2，...，Un。
在使用PCA之前，首先通過從數據集中減去每個要素的平均值，並縮放每個維度使其處於相同的范圍來首先對數據進行歸一化。 在腳本ex7 pca.m中，使用featureNormalize函數為此進行了歸一化。歸一化數據之后，可以運行PCA來計算主要成分。
可以由下式計算數據的協方差矩陣，\(\Sigma = \frac{1}{m}X^TX\)：其中X是以行為單位的矩陣，m是樣本的數目，注意Σ是n×n矩陣，而不是求和運算符。
在計算過協方差矩陣之后，可以使用SVD函數：[U,S,V] = svd(Sigma)，這里U為主要成分，S為對角陣，之后運行可以得到下圖：
pca.m函數代碼如下：
% Useful values
[m, n] = size(X);
% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);
Sigma = X'*X./m;
[U,S,~] = svd(Sigma);
2.3使用PCA維度縮減
在計算中心成分之后，可以使用它們通過將每個樣本投影到較低維度空間，\(x^{(i)}->z^{(i)}\)（例如，將數據從2D投影到1D）。 在本節中將使用PCA返回的特征向量，並將示例數據集投影到一維空間中。
實際上，如果正在使用線性回歸或神經網絡等學習算法，則可以使用投影數據而不是原始數據。 通過使用投影數據，可以更快地訓練模型速度，因為輸入中的size較小。
將數據投影到主成分上
function Z = projectData(X, U, K)
% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);
Z = X*U(:,1:K);
重構數據
function X_rec = recoverData(Z, U, K)
% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));
X_rec = Z*U(:,1:K)';
在完成projectData和recoverData函數之后，ex7 pca.m將會把投影和近似值重建，以顯示投影如何影響數據。 原始數據點用藍色圓圈表示，而投影數據點用紅色圓圈表示。
2.4面部圖片數據集
ex7faces.mat中已經給出了32× 32的灰度面部圖片數據集X，其中X矩陣的每一行表示一張面部圖片，將數據集X的前100行即前100張面部圖片可視化，如下圖：
2.4.1 PCA on Faces
為了在面部圖片的數據集上運行PCA，首先應該先通過從數據矩陣X中減去每個特征的平均值來對數據集進行歸一化處理。
運行PCA后，將獲得數據集的中心成分。 請注意，U（每行）中的每個主成分是長度為n的矢量（面部圖片數據集的n = 1024）。這里我們可以通過將它們中的每一個重新形成-reshape與原始數據集中的像素對應的32×32矩陣來可視化這些計算出來的中心成分，可視化結果如下：(對比原圖，眼眶、鼻子和嘴巴突出的很明顯)
% Before running PCA, it is important to first normalize X by subtracting
% the mean value from each feature
[X_norm, mu, sigma] = featureNormalize(X);%
% Run PCA
[U, S] = pca(X_norm);
% 繪出前36個特征向量
displayData(U(:, 1:36)');
2.4.2維度縮減
現在可以使用上面所計算的面部中心成分來應用與面部數據集上，這樣就可以使用較小尺寸的輸入（例如，100維）來學習算法，而不是原始的1024維度，以加快學習算法的速度。
ex7 pca.m腳本中將面部數據集投影到前100個中心成分上。 具體地，現在每個臉部圖像為向量\(z^{(i))} \in \mathbb{R}^{100}\)。
K = 100;
Z = projectData(X_norm, U, K);
要了解維度降低中丟失的內容，可以使用投影數據集來恢復數據。 在腳本中，執行數據的近似恢復，原始和投影的臉部圖像顯示如下：
K = 100;
X_rec = recoverData(Z, U, K);
% Display normalized data
subplot(1, 2, 1);
displayData(X_norm(1:100,:));
title('Original faces');
axis square;
% Display reconstructed data from only k eigenfaces
subplot(1, 2, 2);
displayData(X_rec(1:100,:));
title('Recovered faces');
axis square;
2.5 PCA在繪圖的應用
在第一節中的K-means圖像壓縮練習中，已經通過對三維RGB空間中使用了K-means算法，這里在的腳本，我們使用scatter3函數（三維散點圖函數）可視化此3D空間中的最終像素分配。每個數據點根據已分配給它的集群進行着色，如上圖所示。
% Reload the image from the previous exercise and run K-Means on it
% For this to work, you need to complete the K-Means assignment first
A = double(imread('bird_small.png'));
% If imread does not work for you, you can try instead
% load ('bird_small.mat');
A = A / 255;
img_size = size(A);
X = reshape(A, img_size(1) * img_size(2), 3);
K = 16;
max_iters = 10;
initial_centroids = kMeansInitCentroids(X, K);
[centroids, idx] = runkMeans(X, initial_centroids, max_iters);
% Sample 1000 random indexes (since working with all the data is
% too expensive. If you have a fast computer, you may increase this.
sel = floor(rand(1000, 1) * size(X, 1)) + 1;
% Setup Color Palette
palette = hsv(K);
colors = palette(idx(sel), :);
% Visualize the data and centroid memberships in 3D
figure;
scatter3(X(sel, 1), X(sel, 2), X(sel, 3), 10, colors);
title('Pixel dataset plotted in 3D. Color shows centroid memberships');
事實證明，可視化數據集在3維及更高的維度可能很麻煩。因此，一般以丟失一些信息為代價，來僅顯示二維數據。在實踐中，PCA通常用於降低數據的維度以進行可視化目的。在ex7 pca.m的腳本中利用PCA實現三維數據的二維可視化，並在2D散點圖中顯示結果。 PCA投影可以被認為是選擇最大化數據擴展的視圖的旋轉，這通常對應於“最佳”視圖。