利用GMM进行无监督face recognition的MATLAB代码及分析

   呵呵,第一次写技术类博文,请多指教。因为做GMM时遇到些小麻烦,所以干脆把问题解决的过程都写下来,帮助学习GMM的童鞋省些工夫。闲话少说,开始。

   使用的数据库是http://www.zjucadcg.cn/dengcai/Data/FaceData.html所提供的。这个数据是YaleB中的10个subject加上Extended YaleB的28个subject,从原来的9个pose里面取出frontal pose,64种光照变化保留,所以一共有38*64 = 2432个。有些图象在获取时损坏被去掉,所以还剩2414个。YaleB_32x32.mat里面包含了所有的脸部图象(即fea,2414x1024的那个,因为每个图象被chop为32x32,所以reshape为1x1024)。YaleB_32x32_corrected.mat里面是经过光照修正的脸部图象。此网站还提供了一些split,就是随机序号,用来生成训练集和测试集,20Train就是每个subject的64个图象中的20个用于训练,依此类推。

   read_data是读入数据并归一化(除以255)。注意归一化很重要,如果不做结果会很差。

   做GMM之前先使用PCA把1024维的特征降到150维,否则难以收敛。coefs是主成分,scores每个sample是投射到低维空间后的坐标。

   显示GMM center时要把mu乘以coefs还原为1024维。

  

   我的GMM是采用了general的假设,即协方差未知、均值未知,但是假设所有的协方差矩阵均为对角阵(即除了对角线元素都为0)。也可以弱化为假设所有类协方差均相等,但是我没有试。关于GMM的原理和计算,请参考CMU教授Andrew Moore的讲义:http://www.autonlab.org/tutorials/gmm.html。我觉得既有趣又好懂。GMM基本上是和EM(Expectation Maximazation)绑在一起的,所以如果不知道什么是EM也请看此文。

   几个关键之处:

   1. 使用kmeans初始化,即GMM的各类中心设为kmeans给出的c个中心;

   2. 各协方差矩阵初始化为单位阵,不能设为kmeans给出的协方差,否则会不收敛(在这个上面忙了很久);

   3. sigma = sigma + 1e-4是防止sigma中出现0,只要加上一个小正数即可,不一定为1e-4,这个数是经验;

   4. thresh,即判定收敛与否的阈值,也是试出来的……但是MATLAB的gmdistribution.fit没有让用户设定阈值,它会自己判断,我想可能是判断error占mu的norm加sigma的norm的百分比吧,不知道它具体怎么实现的。

 

  总而言之,GMM的收敛问题很难办,一个细节没搞好就不收敛了……最后效果还不错,有几张脸和subject很像,但是有些脸基本上是全黑的,无法辨认;这个问题大概是因为GMM把所有光照不足黑乎乎的脸都认为是一类了。无监督学习容易出现这种问题。

 

  下面附GMM代码:

  主程序:

  main_gmm.m

===============================================

[fea_Train fea_Test gnd_Train gnd_Test] = read_data(50, 1);
fea_All = cat(1, fea_Train, fea_Test);

[coefs scores] = princomp(fea_All);  % princomp是MATLAB的PCA函数
coefs_gmm = coefs(:,1:150);     % 降到150维
scores_gmm = scores(:,1:150);
[mu sigma] = my_gmm(scores_gmm, 38);

gmm_face = mu*coefs_gmm';

for i = 1:38
    subplot(4,10,i);
    Y = reshape(gmm_face(i,:),[32 32]);
    imagesc(Y);axis off;colormap(gray)  % colormap会自动做直方图均衡化
    title(['Subject ', num2str(i)]);
    hold on
end

=============================================

  函数:

  my_gmm.m

 

function [mu sigma] = my_gmm(X,c)
% Fit Gaussian Mixture Model for input data X, c classes. For simplicity
% and stability covariance are set to be diagonal.
% Input:
% X:[nxd]. n observations by d variables.
% c: no. of classes.
% Output:
% mu: [cxd]. Each row is the mean of a class.
% sigma:[dxc].

% Initialization
% vl_kmeans can be replaced by MATLAB's kmeans, but this one is
% significantly faster. The toolbox can be downloaded at vlfeat.org.
[mu idx] = vl_kmeans(X',c);
clear idx
mu = mu';
mu_prev = mu;

% A priori probability of each class.
[n d] = size(X);
prior = ones(c,1)/c;
prior_prev = prior;

% A posterior probability of each sample belonging to each class.
post = ones(n,c)/c;
% set sigma as an array representing diagonal matrix
sigma = ones(1,d,c);sigma_prev = sigma;
likeli = zeros(n,c);

maxIter = 100;
thresh = 1;
error = 0;

for t=1:maxIter
    % E-step
    for i = 1:c
        post(:,i)= mvnpdf(X,mu_prev(i,:),sigma_prev(:,:,i))*prior_prev(i);
        likeli(:,i) = mvnpdf(X,mu_prev(i,:),sigma_prev(:,:,i));
    end
    for j = 1:n
        post(j,:) = post(j,:)/sum(post(j,:));
    end   
    % M-step   
    mu = post'*X./repmat(sum(post)',1,d);
    for i = 1:c
       sigma(:,:,i) = diag((X'-repmat(mu(i,:)',1,n))*diag(post(:,i))*(X'-repmat(mu(i,:)',1,n))'/sum(post(:,i)));  % 这个地方当时卡了很久T_T乘后验概率矩阵的时候一定要将其对角化!否则sigma不正定!diag(post(:,i))就是取每个数据属于第i类的后验概率向量(2414x1)将其变为2414x2414维的对角阵,这样X‘PX才与公式中的相等
    end
   
    % Regularize
    sigma = sigma + 1e-4;
    for i = 1:c
        prior(i) = sum(post(:,i))/n;
    end
   
    for i = 1:c
        error = error + norm(mu(i,:)-mu_prev(i,:)) + norm(sigma(:,:,i)-sigma_prev(:,:,i));
    end
   
    if error < thresh
        break
    end
    mu_prev = mu;
    sigma_prev = sigma;
    prior_prev = prior;
    error = 0;
end

fprintf('Iteration times: %d\n',t);

============================================================

  利用gmdistribution.fit函数进行参数估计的结果:

[fea_Train fea_Test gnd_Train gnd_Test] = read_data(50, 1);
fea_All = cat(1, fea_Train, fea_Test);
options = statset('Display', 'final');
obj = gmdistribution.fit(fea_All, 38, 'Options', options, 'CovType', 'diagonal', 'Regularize', 1e-4);
for i = 1:38
    subplot(4,10,i);
     Y = reshape(obj.mu(i,:),[32 32]);
     imagesc(Y);axis off;colormap(gray);
     title(['Subject ', num2str(i)]);
     hold on

end

=============================================================

 

  如果此文对您的学习有任何帮助,欢迎在下面回复告诉我!谢谢! 

posted @ 2010-12-01 21:50  倩倩  阅读(6597)  评论(14编辑  收藏  举报