利用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
=============================================================
如果此文对您的学习有任何帮助，欢迎在下面回复告诉我！谢谢！
