用EM算法求解高斯混合模型

来源：互联网发布：家具三维设计软件编辑：程序博客网时间：2024/05/19 18:14

本文从高斯混合模型出发，引出EM算法，最后回归到利用EM算法求解高斯混合模型。理论部分力求详尽不留证明疑点，所以略显冗长。实验部分给出了生成高斯混合分布样本和利用EM算法求解高斯混合模型的matlab代码。

理论部分

高斯混合模型(GMM)

顾名思义，高斯混合模型就是由多个高斯分布混合构成的模型。K高斯混合分布的概率密度为：

p (x) = \sum k = 1 K ϕ k N (x | μ k, Σ k) .

这里，

∑Kk=1ϕk=1为混合系数，

N (x | μ, Σ) = 1 ( 2 π ) D / 2 1 | Σ | 1 / 2 exp {- 1 2 (x - μ) T Σ - 1 (x - μ)}

为

D维高斯分布，其中

μ为其均值向量，

Σ为其协方差矩阵。

直观来看，高斯混合分布可以看做下面分步过程的整合：
第一步，以ϕk概率选择第k个高斯模型；
第二步，利用第k个高斯模型生成一个样本x。

为了讨论方便，记θ={ϕ1,…,ϕK,μ1,…,μK,Σ1,…,ΣK}为高斯混合分布的参数集合。现在要解决的问题是，对于给定服从高斯混合分布的独立同分布样本集X={x1,…,xn}，最大化其对数似然函数：

max θ ln p (X | θ) = \sum i = 1 n ln \sum k = 1 K ϕ k N (x i | μ k, Σ k) .

这里，

p(X|θ)=∏i=1np(xi|θ)。由于

ln函数里的求和项，我们无法直接求得闭式解，而利用EM算法可以得到一个局部最优数值解。

EM算法

从分步角度来看，求解高斯混合模型的难点在于，我们不知道一个样本xi具体是由K个高斯模型中的哪一个生成的。所以，对于第i个样本xi来说，我们构造一个隐变量zi用来表示xi来自于哪个高斯模型。也即，zi=k当且仅当xi来自于第k个高斯模型。注意，虽然这个隐变量的取值是客观确定的，但对我们来说是不可见，因此仍将其看作随机变量。记Z={z1,…,zn}，对于固定的模型参数θ¯，下面的不等式给出了EM算法的框架。

ln p (X | θ ¯) = ln \sum Z p (X, Z | θ ¯) = ln \sum Z p ( X , Z | θ ¯ ) p ( Z | X , θ ¯ ) p (Z | X, θ ¯) = \sum Z p (Z | X, θ ¯) ln p ( X , Z | θ ¯ ) p ( Z | X , θ ¯ ) \leq max θ \sum Z p (Z | X, θ ¯) ln p ( X , Z | θ ) p ( Z | X , θ ¯ ) \leq ln \sum Z p ( X , Z | θ m a x ) p ( Z | X , θ ¯ ) p (Z | X, θ ¯) = ln p (X | θ m a x) (p ( X , Z | θ ¯ ) p ( Z | X , θ ¯ ) = p (X | θ ¯) 为 一 常 数 ， 故 等 号 成 立) (1) (J e n s e n 不 等 式)

如果将

θ¯和

θmax分别看成EM算法一次迭代前后的模型参数，上面的推导实际上证明了其会收敛到局部最优解。对于上面的式(1)，我们有

θ m a x = arg max θ \sum Z p (Z | X, θ ¯) ln p ( X , Z | θ ) p ( Z | X , θ ¯ ) = arg max θ \sum Z p (Z | X, θ ¯) ln p (X, Z | θ) . (2)

这里式(2)中的优化问题就对应了EM算法中的M步。另一方面，E步指的就是利用旧的模型参数

θ¯将目标函数

∑Zp(Z|X,θ¯)lnp(X,Z|θ)显示的表达出来。
对于EM算法中的M步的目标函数，在《PRML》和Ng的讲义中分别用两种形式来表达，下面的命题表明两种形式是等价的。注意到

p(Z|X,θ¯)=∏j=1np(zj|xj,θ¯)以及

p(X,Z|θ)=∏i=1np(xi,zi|θ)，则有

Proposition 1.

\sum Z p (Z | X, θ ¯) ln p (X, Z | θ) = \sum i = 1 n \sum z i p (z i | x i, θ ¯) ln p (x i, z i | θ) . (3)

Pf.

\sum Z p (Z | X, θ ¯) ln p (X, Z | θ) = \sum Z \prod j = 1 n p (z j | x j, θ ¯) ln \prod i = 1 n p (x i, z i | θ) = \sum z 1 \sum z 2 \dots \sum z n \prod j = 1 n p (z j | x j, θ ¯) \sum i = 1 n ln p (x i, z i | θ) = \sum i = 1 n \sum z i p (z i | x i, θ ¯) ln p (x i, z i | θ) \sum z 1 \sum z 2 \dots \sum z t, t \neq i \dots \sum z n \prod j = 1, j \neq i n p (z j | x j, θ ¯) = \sum i = 1 n \sum z i p (z i | x i, θ ¯) ln p (x i, z i | θ) \prod j = 1, j \neq i n \sum z j p (z j | x j) = \sum i = 1 n \sum z i p (z i | x i, θ ¯) ln p (x i, z i | θ) . (多 项 式 乘 积 展 开 式 的 逆) (概 率 归 一 性) □

窃以为式(3)中左式(PRML中的写法)简洁抽象，适合做理论推导；右式(Ng讲义中的写法)表达具体，适合指导算法实践。

EM算法求解高斯混合模型

具体应用到GMM模型中，EM算法中的M步我们要求解

max θ \sum i = 1 n \sum z i p (z i | x i, θ ¯) ln p (x i, z i | θ) = \sum i = 1 n \sum z i = 1 K p (z i | x i, θ ¯) ln ϕ z i N (x i | μ z i, Σ z i) = \sum i = 1 n \sum k = 1 K p (k | x i, θ ¯) ln ϕ k N (x i | μ k, Σ k)

注意到

p(k|xi,θ¯)不包含优化变量

θ，因此在此问题中可以看做关于

i和

k的常数，故我们暂记其为

γik，留到最后再讨论其具体取值。忽略掉优化无关常数，M步的优化问题正式写作：

max ϕ k, μ k, Σ k f (ϕ k, μ k, Σ k) = \sum i = 1 n \sum k = 1 K γ i k (ln ϕ k - 1 2 ln | Σ k | - 1 2 (x i - μ k) T Σ - 1 k (x i - μ k)) (4)

s.t. ϕ k > 0 \sum k = 1 K ϕ k = 1 Σ k ≻ 0 (k = 1, \dots, K) (k = 1, \dots, K)

这里，

Σk≻0表示

Σk正定。观察到对

μk没有任何约束，所以我们先让目标函数对

μk求偏导并令其为0

\partial \partial μ k = - \sum i = 1 n γ i k Σ - 1 k (x i - μ k) = 0,

令

nk=∑i=1nγik，利用

Σk满秩约束解得

μ * k = 1 n k \sum i = 1 n γ i k x i .

接下来，我们在没有

Σk≻0限制下求使目标函数极大的

Σ∗k，最后再验证

Σ∗k的确是一个正定矩阵。因为直接对

Σk求偏导比较困难，因此我们转而对

Σ−1k求偏导并令其为0

\partial \partial Σ - 1 k = \sum i = 1 n γ i k (adj ( Σ - 1 k ) T | Σ - 1 k | - (x i - μ k) (x i - μ k) T) = 0,

解得

Σ * k = 1 n k \sum i = 1 n γ i k (x i - μ k) (x i - μ k) T .

这里

adj(M)=|M|M−1表示矩阵

M的伴随矩阵。容易验证

Σ∗k为一个半正定矩阵。

Σ∗k不可逆的情况只会在样本集的防射维数小于数据空间维数情况下发生，对于此情况有两种解决办法。一是将原始数据利用PCA降维；二是对数据加上微小随机扰动。两种方式本质上都是保证数据集的防射维数等于数据空间维数。

然后，我们对ϕk在条件∑Kk=1ϕk=1下利用拉格朗日乘子法求解，然后再来验证ϕk>0。构造拉格朗日函数

L (ϕ k, λ) = f (ϕ k, μ k, Σ k) + λ (\sum k = 1 K ϕ k - 1),

利用KKT条件：

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ \partial \partial ϕ k L (ϕ k, λ) = 1 ϕ k \sum i = 1 n γ i k + λ = 0 \sum k = 1 K ϕ k = 1. k = 1, \dots, K

解得

ϕ * k = n k n,

可以看到

ϕk确实是大于0的。这样就验证了

μ∗k,Σ∗k,ϕ∗k确实是满足条件的一个局部最优解。最后再来看

γik=p(k|xi,θ¯)的具体表达式，利用贝叶斯公式可得

p (k | x i, θ ¯) = ϕ ¯ k N ( x i | μ ¯ k , Σ ¯ k ) \sum k = 1 K ϕ ¯ k N ( x i | μ ¯ k , Σ ¯ k ) .

综上所述，可得整个算法如下：
这里写图片描述

实验部分

生成数据

此步用matlab生成服从高斯混合分布的数据。
这里写图片描述

%参数区，可自行调整K = 3;      %高斯模型个数dim = 3;    %数据维度N = 1000;    %生成的样本个数step = 4;   %不同高斯模型均值之间的距离%变量区Phi = cell(K,1);         %混合系数构成的cellArrayMu=cell(K,1);           %均值向量构成的cellArraySigma = cell(K,1);      %协方差矩阵构成的cellArrayDat = zeros(N,dim+1);   %数据矩阵，其中最后一个维度为标签for i =1:K    Phi(i)={1/K};    %默认每个模型的选择系数相同，可自行调整endMu(1) = {ones(dim,1)};for i =2:K    t = rand(dim,1);    t = t/norm(t,2)*step;    Mu(i)={Mu{i-1}+t};endfor i =1:K    t = rand(dim);    t = orth(t);    tt = rand(dim,1);    Sigma(i) = {t'*diag(tt)*t};end%这里直接按系数比例生成样本，与直接按照高斯混合分布生成的样本大致是类似的cnt = 0;for i = 1:K    if i == K        n = N-cnt;    else        n = round(N*Phi{i});    end    t = mvnrnd(Mu{i}, Sigma{i}, n);    Dat(cnt+1:cnt+n,:) = [t ones(n, 1)*i];    cnt = cnt+n;end%存储数据部分save('Dat.mat', 'Dat');%画图部分close all;view(3);hold on;grid on;mycolor = colorcube(20);cnt = 0;for i = 1:K    if i == K        n = N-cnt;    else       n = round(N*Phi{i});    end    t =  Dat(cnt+1:cnt+n,:);    plot3(t(:,1), t(:,2), t(:,3), '.', 'Color', mycolor(i,:));    cnt = cnt+n;endhold off;

EM算法求解GMM

利用EM算法求解GMM的过程中涉及到参数初始化问题，比较高效的做法是利用Kmeans先做聚类，再把每一类样本的比例，均值和协方差矩阵作为初始参数，再利用EM算法进行微调。如果想观察EM算法求解GMM的收敛行为，可以使用随机初始化参数。在下面的代码中，我们利用迭代前后对数似然函数的相对偏差作为判断其是否收敛。
这里写图片描述

function []=main()    close all;    Eps = 1e-4;  %当迭代前后两次对数似然函数值的相对偏差小于Eps时算法结束    flag = false; %是否使用kmeans初始化参数    load('Dat.mat'); %导入数据，Dat为N*(dim+1)维的矩阵，每行最后一个元素为其隐变量取值(实验中没有用到)    N = size(Dat, 1);   %样本个数    dim = size(Dat, 2)-1;   %数据维度    K = max(Dat(:,dim+1));  %高斯模型个数    %初始化参数部分    if flag        %用kmeans初始化参数        [PhiBar, MuBar, SigmaBar] = initByKmeans(Dat(:,1:dim), K);        else        %随机初始化参数        PhiBar = cell(K,1);                MuBar=cell(K,1);                  SigmaBar = cell(K,1);          rPhi = rand(3,1);        rPhi = rPhi/(ones(1, 3)*rPhi);        for k=1:K            PhiBar(k) = {rPhi(k)};            MuBar(k) = {mean(Dat(:,1:dim))'.*(1+rand(dim, 1)*1)};            SigmaBar(k) = {eye(dim)};        end    end    logLikeBar = calLogLike(Dat(:,1:dim), K, PhiBar, MuBar, SigmaBar); %计算初始对数似然函数值    Phi= cell(K,1);            Mu=cell(K,1);              Sigma = cell(K,1);     %EM算法    cnt = 0;    while true        %E step:        Gamma = zeros(N, K);        for i = 1:N            Xi = Dat(i,1:dim)';            Pxi = 0;            for k = 1:K                Pxi = Pxi+PhiBar{k}*mvnpdf(Xi, MuBar{k}, SigmaBar{k});            end            for k = 1:K                Gamma(i,k) = PhiBar{k}*mvnpdf(Xi, MuBar{k}, SigmaBar{k})/Pxi;            end        end        %M step:        N_k = ones(1,N)*Gamma;        for k=1:K            Phi(k) = {N_k(k)/N};            mu = zeros(dim, 1);            for i = 1:N                mu = mu+Gamma(i,k)*Dat(i,1:dim)';            end            Mu(k) = {mu/N_k(k)};            sigma = zeros(dim);            for i = 1:N                sigma = sigma+Gamma(i,k)*(Dat(i,1:dim)'-Mu{k})*(Dat(i,1:dim)-Mu{k}');            end            Sigma(k) = {sigma/N_k(k)};        end        cnt = cnt+1;        logLike = calLogLike(Dat(:,1:dim), K, Phi, Mu, Sigma);        relGap = abs(logLike-logLikeBar)/abs(logLikeBar);        disp(['Iteration: ', num2str(cnt), ' RelativeGap: ', sprintf('%.6f LogLikelihood %.6f', relGap, logLike)]);        logLikeBar = logLike;        %绘图部分%         figure(cnt);%         view(3);%         hold on;%         grid on;%         mycolor = colorcube(20);%         for i = 1:N%             postP = 0;%             choose = 0;%             for k = 1:K %利用后验概率最大化准则%                 temp = Phi{k}*mvnpdf(Dat(i,1:dim)', Mu{k}, Sigma{k});%                 if temp > postP%                     postP = temp;%                     choose = k;%                 end%             end%             plot3(Dat(i,1), Dat(i,2), Dat(i,3),'.', 'Color', mycolor(choose,:), 'markersize', 4);%         end%         for k =1:K%             plot3(Mu{k}(1), Mu{k}(2), Mu{k}(3), 'k.', 'markersize', 15);%         end%         saveas(gcf,sprintf('%d.jpg', cnt));        if relGap < Eps            break;        end        PhiBar = Phi;        MuBar = Mu;        SigmaBar = Sigma;    endendfunction [Phi, Mu, Sigma] = initByKmeans(X, K) %Kmeans随机初始化参数    N = size(X,1);    dim = size(X,2);    Phi= cell(K,1);            Mu=cell(K,1);              Sigma = cell(K,1);      for k=1:K        Phi(k) = {0};        Mu(k) = {zeros(dim, 1)};        Sigma(k) = {zeros(dim)};    end    [idx, C] = kmeans(X, K);    for k = 1:K        Mu{k} = C(k,:)';    end    for i = 1:N        k = idx(i);        Phi{k} = Phi{k}+1;        Sigma{k} = Sigma{k}+(X(i,:)'-Mu{k})*(X(i,:)-Mu{k}');    end    for k=1:K        Sigma{k}= Sigma{k}/Phi{k}+1e-6*ones(dim);        Phi{k} = Phi{k}/N;    endendfunction [logLike] = calLogLike(X, K, Phi, Mu, Sigma) %计算对数似然函数值    logLike = 0;    N = size(X,1);    for i =1:N        p = 0;        for k = 1:K            p = p+Phi{k}*mvnpdf(X(i,:)', Mu{k}, Sigma{k});        end        logLike = logLike+log(p);    endend

参考文献

Christopher M.. Bishop. Pattern recognition and machine learning. pp. 430-459

Andrew Ng. Machine Learning Lecture Notes.

2 0