Deep learning：二十五(Kmeans單層網絡識別性能)

 

　　前言：

　　本文是用kmeans方法來分析單層網絡的性能，主要是用在CIFAR-10圖像識別數據庫上。關於單層網絡的性能可以參考前面的博文：Deep learning：二十(無監督特征學習中關於單層網絡的分析)。當然了，本文依舊是參考論文An Analysis of Single-Layer Networks in Unsupervised Feature Learning, Adam Coates, Honglak Lee, and Andrew Y. Ng. In AISTATS 14, 2011.只是重點在分析4個算法中的kemans算法（因為作者只提供關於kmeans的demo，呵呵，當然了另一個原因是sparse autoencoder在前面的博文中也介紹很多了）本文的代碼可以在Honglak Lee主頁中下載：http://web.eecs.umich.edu/~honglak/

 

　　實驗基礎：

　　Kmeans相關：

　　Kmeans可以分為2個步驟，第一步是cluster assignment step，就是完成各個樣本的聚類。第二步是move centroid，即重新選定類別中心點。Kmeans聚類不僅可以針對有比較明顯類別的數據，還可以針對不具有明顯類別的數據（即人眼看起來根本就沒有區別），即使是沒明顯區分的數據用kmeans聚類時得到的結果也是可以進行解釋的，因為有時候在某種原因下類別數是人定的。

　　既然kmeans是一種機器學習算法，那么它肯定也有一個目標函數需要優化，其目標函數如下所示：

 　　

　　在kmeans初始化k個類別時，由於初始化具有隨機性，如果選取的初始值點不同可能導致最后聚類的效果跟想象中的效果相差很遠，這也就是kmeans的局部收斂問題。解決這個問題一般采用的方法是進行多次kmeans，然后計算每次kmeans的損失函數值，取損失函數最小對應的那個結果作為最終結果。

　　在kmeans中比較棘手的另一個問題是類別k的選擇。因為有的數據集用不同的k來聚類都感覺比較合適，那么到底該用哪個k值呢？通常情況下的方法都是采用”elbow”的方法，即做一個圖表，該圖的橫坐標為選取的類別個數k，縱坐標為kmeans的損失函數，通過觀察該圖找到曲線的轉折點，一般這個圖長得像人的手，而那個像人手肘對應的轉折點就是我們最終要的類別數k，但這種方法也不一定合適，因為k的選擇可以由人物確定，比如說我就是想把數據集分為10份（這種情況很常見，比如說對患者年齡進行分類），那么就讓k等於10。

　　在本次試驗中的kmeans算法是分為先求出每個樣本的聚類類別，然后重新計算中心點這2個步驟。但是在求出每個樣本的聚類類別是不是簡單的計算那2個向量的歐式距離。而是通過內積實現的。我們要A矩陣中a樣本和B矩陣中所有樣本（此處用b表示）距離最小的一個求，即求min(a-b)^2，等價於求min(a^2+b^2-2*a*b)，等價於求max(a*b-0.5*a^2-0.5*b^2),假設a為輸入數據中固定的一個， b為初始化中心點樣本中的某一個，則固定的a和不同的b作比較時，此時a中的該數據可以忽略不計，只跟b有關。即原式等價於求max(a*b-0.5*b^2)。也就是runkmeans函數的核心思想。（這個程序一開始沒看懂，后面慢慢推算總算弄明白了，應該是它這樣通過矩陣操作進行kmeans距離的速度比較快吧！）

　　當通過聚類的方法得到了樣本的k個中心以后就要開始提取樣本的特征了，當然了這些樣本特征的提取是根據每個樣本到這k個類中心點的距離構成的，最簡單的方法就是取最近鄰，即取於這k個類別中心距離最近的那個類為類標簽1，其它都為0，其計算公式如下：

 　　

　　因為那樣計算就有很高的稀疏性（只有1個為1，其它 都為0），而如果需要放松條件則可以這樣考慮：先計算出對應樣本與k個類中心點的平均距離d，然后如果那些樣本與類別中心點的距離大於d的話都設置為0，小於d的則用d與該距離之間的差來表示。這樣基本能夠保證一半以上的特征都變成0了，也是具有稀疏性的，且考慮了更多那些距類別中心距離比較近的值。此時的計算公式如下：

 　　

　　首先是關於CIFAR-10的數據庫，到網站上http://www.cs.toronto.edu/~kriz/下載的CIFAR-10數據庫解壓后如下：

 　　

　　其中的每個data_batch都是10000x3072大小的，即有1w個樣本圖片，每個圖片都是32*32且rgb三通道的，這里的每一行表示一個樣本，與前面博文程序中的剛好相反。因為總共有5個data_batch，所以共有5w張訓練圖片。而測試數據test_batch則有1w張，是分別從10類中每類隨機選取1000張。

　　關於均值化的一點總結：

　　給定多張圖片構成的一個矩陣（其中每張圖片看成是一個向量，多張圖片就可以看做是一個矩陣了）。要對這個矩陣進行whitening操作，而在這之前是需要均值化的。在以前的實驗中，有時候是對每一張圖片內部做均值，也就是說均值是針對每張圖片的所有維度，而有的時候是針對矩陣中圖片的每一維做均值操作，那么是不是有矛盾呢？其實並不矛盾，主要是這兩種均值化的目的不同。如果是算該均值的協方差矩陣，或者將一些訓練樣本輸入到分類器訓練前，則應該對每一維采取均值化操作（因為協方差均值是描述每個維度之間的關系）。如果是為了增強每張圖片亮度的對比度，比如說在進行whitening操作前，則需要對圖片的內部進行均值化（此時一般還會執行除以該圖像內部的標准差操作）。

　　另外，一般輸入svm分類器中的樣本都是需要標准化過。

　　Matlab相關：

　　Matlab中function函數內部並不需要針對function有個end語句。

　　svd(),eig()：

　　其實按照道理這2者之間應該是完全不同的。相同之處是這2個函數的輸入參數必須都是方陣。

　　cov:

　　cov(x)是求矩陣x的協方差矩陣。但對x是有要求，即x中每一行為一個樣本，也就是說每一列為數據的一個維度值，不要求x均值化過。

　　var:

　　該函數是用來求方差的，求方差時如果是無偏估計則分母應該除以N-1,否則除以N即可。默認情況下分母是除以N-1的，即默認采用的是無偏估計。

　　b1 = var(a); % 按默認來求
　　b2 = var(a, 0); % 默認的公式（除以N-1）
　　c1 = var(a, 1); % 另外的公式（除以N）
　　d1 = var(a, 0, 1); % 對每列操作（除以N-1）
　　d2 = var(a, 0, 2); % 對每行操作（除以N-1）。

　　Im2col:

　　該函數是將一個大矩陣按照小矩陣取出來，並把取出的小矩陣展成列向量。比如說B = im2col(A,[m n],block_type):就是把A按照m*n的小矩陣塊取出，取出后按照列的方式重新排列成向量,然后多個列向量組成一個矩陣。而參數block_type表示的是取出小矩形框的方式，有兩種值可以取，分別為’distinct’和’sliding’。Distinct方式是指在取出的各小矩形在原矩陣中是沒有重疊的，元素不足的補0。而sliding是每次移動一個元素，即各小矩形之間有元素重疊，但此時沒有補0元素的說法。如果該參數不給出，則默認的為’sliding’模式。

　　random:

　　該函數和常見的rand,randi,randn不同，random可以產生各種不同的分布，其不同分布由參賽name決定，比如二項分布，泊松分布，指數分布等，其一般的調用形式為：   Y = random(name,A,B,C,[m,n,...])

　　rdivide：

　　在bsxfun(@rdivide,A,B)中，其中A是一個矩陣，B是一個行向量，則該函數的意思是將A中每個元素分別除以在B中對應列的值。

　　sum:

　　這里主要是想說進行多維矩陣的求sum操作，比如矩陣X為m*n*p維的，則sum(X,1)計算出的結果是1*n*p維的，而sum(x,2)后得到的尺寸是m*1*p維，sum(x,3) 后得到的尺寸是m*n*1維，也就是說，對哪一維求sum，則計算得到的結果后的那一維置1即可，其它可保持不變。

 

 

　　實驗結果：

　　kemans學習到的類中心點圖片顯示如下：

　　

　　用kmeans方法對CIFAR-10訓練圖片的識別效果如下

　　Train accuracy 86.112000%

　　對測試圖片識別的效果如下：

　　Test accuracy 77.350000%

 

　　實驗主要部分代碼：

kmeans_demo.m:

CIFAR_DIR='cifar-10-batches-mat/';

assert(strcmp(CIFAR_DIR, 'cifar-10-batches-mat/'), ...%strcmp相等時為1
       ['You need to modify kmeans_demo.m so that CIFAR_DIR points to ' ...
        'your cifar-10-batches-mat directory.  You can download this ' ...
        'data from:  http://www.cs.toronto.edu/~kriz/cifar-10-matlab.tar.gz']);

%% Configuration
addpath minFunc;
rfSize = 6;
numCentroids=1600;%類別總數
whitening=true;
numPatches = 400000;%40w張圖片，不少啊！
CIFAR_DIM=[32 32 3];

%% Load CIFAR training data
fprintf('Loading training data...\n');
f1=load([CIFAR_DIR '/data_batch_1.mat']);
f2=load([CIFAR_DIR '/data_batch_2.mat']);
f3=load([CIFAR_DIR '/data_batch_3.mat']);
f4=load([CIFAR_DIR '/data_batch_4.mat']);
f5=load([CIFAR_DIR '/data_batch_5.mat']);

trainX = double([f1.data; f2.data; f3.data; f4.data; f5.data]);%50000*3072
trainY = double([f1.labels; f2.labels; f3.labels; f4.labels; f5.labels]) + 1; % add 1 to labels!,變成1到10
clear f1 f2 f3 f4 f5;%及時清除變量

% extract random patches
patches = zeros(numPatches, rfSize*rfSize*3);%400000*108
for i=1:numPatches
    i=1;
  if (mod(i,10000) == 0) fprintf('Extracting patch: %d / %d\n', i, numPatches); end
  
  r = random('unid', CIFAR_DIM(1) - rfSize + 1);%符合均一分布
  c = random('unid', CIFAR_DIM(2) - rfSize + 1);
  %使用mod(i-1,size(trainX,1))是因為對每個圖片樣本，提取出numPatches/size(trainX,1)個patch
  patch = reshape(trainX(mod(i-1,size(trainX,1))+1, :), CIFAR_DIM);%32*32*3
  patch = patch(r:r+rfSize-1,c:c+rfSize-1,:);%6*6*3
  patches(i,:) = patch(:)';%patches的每一行代表一個小樣本
end

% normalize for contrast，亮度對比度均一化，減去每一行的均值然后除以該行的標准差（其實是標准差加10）
%bsxfun(@rdivide,A,B)表示A中每個元素除以B中對應行（或列）的值。
patches = bsxfun(@rdivide, bsxfun(@minus, patches, mean(patches,2)), sqrt(var(patches,[],2)+10));

% whiten
if (whitening)
  C = cov(patches);%計算patches的協方差矩陣
  M = mean(patches);
  [V,D] = eig(C);
  P = V * diag(sqrt(1./(diag(D) + 0.1))) * V';%P是ZCA Whitening矩陣
  %對數據矩陣白化前，應保證每一維的均值為0
  patches = bsxfun(@minus, patches, M) * P;%注意patches的行列表示的意義不同時，白化矩陣的位置也是不同的。
end

% run K-means
centroids = run_kmeans(patches, numCentroids, 50);%對樣本數據patches進行聚類，聚類結果保存在centroids中
show_centroids(centroids, rfSize); drawnow;

% extract training features
if (whitening)
  trainXC = extract_features(trainX, centroids, rfSize, CIFAR_DIM, M,P);%M為均值向量，P為白化矩陣，CIFAR_DIM為patch的維數，rfSize為小patch的大小
else
  trainXC = extract_features(trainX, centroids, rfSize, CIFAR_DIM);
end

% standardize data，保證輸入svm分類器中的數據都是標准化過了的
trainXC_mean = mean(trainXC);
trainXC_sd = sqrt(var(trainXC)+0.01);
trainXCs = bsxfun(@rdivide, bsxfun(@minus, trainXC, trainXC_mean), trainXC_sd);
trainXCs = [trainXCs, ones(size(trainXCs,1),1)];%每一個特征后面都添加了一個常量1

% train classifier using SVM
C = 100;
theta = train_svm(trainXCs, trainY, C);

[val,labels] = max(trainXCs*theta, [], 2);
fprintf('Train accuracy %f%%\n', 100 * (1 - sum(labels ~= trainY) / length(trainY)));

%%%%% TESTING %%%%%

%% Load CIFAR test data
fprintf('Loading test data...\n');
f1=load([CIFAR_DIR '/test_batch.mat']);
testX = double(f1.data);
testY = double(f1.labels) + 1;
clear f1;

% compute testing features and standardize
if (whitening)
  testXC = extract_features(testX, centroids, rfSize, CIFAR_DIM, M,P);
else
  testXC = extract_features(testX, centroids, rfSize, CIFAR_DIM);
end
testXCs = bsxfun(@rdivide, bsxfun(@minus, testXC, trainXC_mean), trainXC_sd);
testXCs = [testXCs, ones(size(testXCs,1),1)];

% test and print result
[val,labels] = max(testXCs*theta, [], 2);
fprintf('Test accuracy %f%%\n', 100 * (1 - sum(labels ~= testY) / length(testY)));

　　

run_kmeans.m:

function centroids = runkmeans(X, k, iterations)

  x2 = sum(X.^2,2);%每一個樣本元素的平方和，x2這里指每個樣本點與原點之間的歐式距離。
  centroids = randn(k,size(X,2))*0.1;%X(randsample(size(X,1), k), :); 程序中傳進來的k為1600，即有1600個聚類類別
  BATCH_SIZE=1000;
  
  
  for itr = 1:iterations%iterations為kemans聚類迭代的次數
    fprintf('K-means iteration %d / %d\n', itr, iterations);
    
    c2 = 0.5*sum(centroids.^2,2);%c2表示類別中心點到原點之間的歐式距離

    summation = zeros(k, size(X,2));
    counts = zeros(k, 1);
    
    loss =0;
    
    for i=1:BATCH_SIZE:size(X,1) %X輸入的參數為50000，所以該循環能夠進行50次
      lastIndex=min(i+BATCH_SIZE-1, size(X,1));%lastIndex=1000,2000,3000,...
      m = lastIndex - i + 1;%m=1000,2000,3000,...
      %這種算法也是求每個樣本的標簽，因為求min(a-b)^2等價於求min(a^2+b^2-2*a*b)等價於求max(a*b-0.5*a^2-0.5*b^2),假設a為輸入數據矩陣，而b為初始化中心點樣本
      %則每次從a中取出一個數據與b中所有中心點作比較時，此時a中的該數據可以忽略不計，只跟b有關。即原式等價於求max(a*b-0.5*a^2)
      [val,labels] = max(bsxfun(@minus,centroids*X(i:lastIndex,:)',c2));%val為BATCH_SIZE大小的行向量（1000*1），labels為每個樣本經過一次迭代后所屬的類別標號
      loss = loss + sum(0.5*x2(i:lastIndex) - val');%求出loss沒什么用
      
      S = sparse(1:m,labels,1,m,k,m); % labels as indicator matrix，最后一個參數為最大非0個數
      summation = summation + S'*X(i:lastIndex,:);%1600*108
      counts = counts + sum(S,1)';%1600*1的列向量，每個元素代表屬於該類樣本的個數
    end


    centroids = bsxfun(@rdivide, summation, counts);%步驟2，move centroids
    
    % just zap empty centroids so they don't introduce NaNs everywhere.
    badIndex = find(counts == 0);
    centroids(badIndex, :) = 0;%防止出現無窮大的情況
  end

 

extract_features.m:

function XC = extract_features(X, centroids, rfSize, CIFAR_DIM, M,P)
  assert(nargin == 4 || nargin == 6);
  whitening = (nargin == 6);
  numCentroids = size(centroids,1);%numCentroids中心點的個數
  
  % compute features for all training images
  XC = zeros(size(X,1), numCentroids*4);%為什么是4呢？因為后面是分為4個區域來pooling的
  for i=1:size(X,1)
    if (mod(i,1000) == 0) fprintf('Extracting features: %d / %d\n', i, size(X,1)); end
    
    % extract overlapping sub-patches into rows of 'patches'
    patches = [ im2col(reshape(X(i,1:1024),CIFAR_DIM(1:2)), [rfSize rfSize]) ;%類似於convolution一樣取出小的patches,patches中每一行都對應原圖中一個小圖像塊的rgb
                im2col(reshape(X(i,1025:2048),CIFAR_DIM(1:2)), [rfSize rfSize]) ;%因此patches中每一行也代表一個rgb樣本，每一行108維，每一張大圖片在patches中占27*27行
                im2col(reshape(X(i,2049:end),CIFAR_DIM(1:2)), [rfSize rfSize]) ]';

    % do preprocessing for each patch
    
    % normalize for contrast，whitening前對每一個樣本內部做均值
    patches = bsxfun(@rdivide, bsxfun(@minus, patches, mean(patches,2)), sqrt(var(patches,[],2)+10));
    % whiten
    if (whitening)
      patches = bsxfun(@minus, patches, M) * P;
    end
    
    % compute 'triangle' activation function
    xx = sum(patches.^2, 2);
    cc = sum(centroids.^2, 2)';
    xc = patches * centroids';
    
    z = sqrt( bsxfun(@plus, cc, bsxfun(@minus, xx, 2*xc)) ); % distances = xx^2+cc^2-2*xx*cc;
    [v,inds] = min(z,[],2);%中間的那個中括號不能少，否則會認為是將z中元素同2比較，現在的2表示z中的第2維
    mu = mean(z, 2); % average distance to centroids for each patch
    patches = max(bsxfun(@minus, mu, z), 0);%patches中每一行保存的是：小樣本與這1600個類別中心距離的平均值減掉與每個類別中心的距離，限定最小距離為0
    % patches is now the data matrix of activations for each patch
    
    % reshape to numCentroids-channel image
    prows = CIFAR_DIM(1)-rfSize+1;
    pcols = CIFAR_DIM(2)-rfSize+1;
    patches = reshape(patches, prows, pcols, numCentroids);
    
    % pool over quadrants
    halfr = round(prows/2);
    halfc = round(pcols/2);
    q1 = sum(sum(patches(1:halfr, 1:halfc, :), 1),2);%求區域內像素之和，是個列向量，1600*1
    q2 = sum(sum(patches(halfr+1:end, 1:halfc, :), 1),2);
    q3 = sum(sum(patches(1:halfr, halfc+1:end, :), 1),2);
    q4 = sum(sum(patches(halfr+1:end, halfc+1:end, :), 1),2);
    
    % concatenate into feature vector
    XC(i,:) = [q1(:);q2(:);q3(:);q4(:)]';%類似於pooling操作
  end

　　

train_svm.m:

function theta = train_svm(trainXC, trainY, C)
  
  numClasses = max(trainY);
  %w0 = zeros(size(trainXC,2)*(numClasses-1), 1);
  w0 = zeros(size(trainXC,2)*numClasses, 1);
  w = minFunc(@my_l2svmloss, w0, struct('MaxIter', 1000, 'MaxFunEvals', 1000), ...
              trainXC, trainY, numClasses, C);

  theta = reshape(w, size(trainXC,2), numClasses);
  
% 1-vs-all L2-svm loss function;  similar to LibLinear.
function [loss, g] = my_l2svmloss(w, X, y, K, C)
  [M,N] = size(X);
  theta = reshape(w, N,K);
  Y = bsxfun(@(y,ypos) 2*(y==ypos)-1, y, 1:K);

  margin = max(0, 1 - Y .* (X*theta));
  loss = (0.5 * sum(theta.^2)) + C*mean(margin.^2);
  loss = sum(loss);  
  g = theta - 2*C/M * (X' * (margin .* Y));
  g = g(:);

  %[v,i] = max(X*theta,[],2);
  %sum(i ~= y) / length(y)

 

 

　　參考資料：

     Deep learning：二十(無監督特征學習中關於單層網絡的分析)

　　An Analysis of Single-Layer Networks in Unsupervised Feature Learning, Adam Coates, Honglak Lee, and Andrew Y. Ng. In AISTATS 14, 2011.

      http://www.cs.toronto.edu/~kriz/

　　http://ai.stanford.edu/~ang/papers.php