Multi-Exemplar Affinity Propagation

AP算法存在的固有缺陷是它不能对包含很多子类的category进行建模，而在image categorization, face categorization, 多字体optical character recognition，手写数字分类，每个category可能包含很多子类。
比如说，在自然场景类别可能包含很多主题，比如说，街景中可能包含一些主题，比如road，car，pedestrian，building等等。

在OCR和手写数字分类问题中，代表letter或digit的类可能由多个子类组成，每个子类代表着不同的style或者字体。显然，我们用一个代表点来统一表示这些子类是不合理的。
这篇论文提出了MEAP算法。每一个cluster都是自动决定exemplars和superexemplars的数目，每个数据点都自动分配给最接近的exemplar,而每个exemplar都分配给最接近的superexemplar。
superexemplar定义为代表一类cluster的所有exemplars中的最具代表性的那一个。

所以目标函数是最大化数据点和代表点(exemplar)之间的相似度，以及exemplar和superexemplar之间的相似度。
直接求解这个问题是NP困难的。所以我们可以使用max-sum belief propagation，可以产生对初始化信息不敏感的算法。

AP算法复习

在之前的博客中已经有过具体的介绍，为了数学符号表示的一致性，这里用新的数学表达式再表示一次。

给定一个用户定义的N个数据点的相似矩阵[sij]N×N[s_{ij}]_{N\times N}，我们的目标是获得一个labels的valid configurationc=[c1,c2,...,cN]c=[c_1,c_2,...,c_N]，来优化目标函数：

S(c)=∑i=1Nsici+∑k=1Nδk(c)

S(c)=\sum_{i=1}^{N}s_{ic_i}+\sum_{k=1}^{N}\delta_k(c)
而 δk(c)\delta_k(c)是一个exemplar-consistency约束，就是说如果有数据点i选择k作为代表点，也就是满足 ci=kc_i=k,那么k必须同时也是自己的代表点，也就是 ck=kc_k=k

δk(c)={−∞,if ck≠k but ∃ i:ci=k0, otherwise

\delta_k(c)=\left\{ \begin{aligned} -\infty , if \ c_k\neq k \ but \ \exists \ i: c_i =k \\ 0, \ otherwise \end{aligned} \right.

更新迭代的过程可以浓缩到下面几行：

最后的类别向量c=[c1,...,cN]c=[c_1,...,c_N]的计算方法为：

ci=arg maxj[a(i,j)+r(i,j)]

c_i = arg \ max_j[a(i,j)+r(i,j)]

MEAP

假设[sij]N×N[s_{ij}]_{N\times N}是一个用户定义的相似度矩阵，sijs_{ij}代表着数据点i和代表点j之间的相似度，[lij]N×N[l_{ij}]_{N\times N}代表着代表点i和潜在的super-exemplar j之间的相似度。
我们的目标是最大化所有数据点和它们对应的代表点之间的相似度S1S_1，同时最大化代表点和superexemplar之间的相似度S2S_2。

模型

假设C是assignment 矩阵，非对角线元素代表数据点j是数据点i的代表点cij=1c_{ij}=1,也就是ψ1(i)=j\psi_1(i)=j, 对角线元素代表 cii∈{1,...,N}c_{ii} \in \{1,...,N\} 代表着ciic_{ii}是代表点i的superexemplar。

S1=∑i=1N∑j=1Nsij⋅[cij≠0],S2=∑i=1Nlicii⋅[cii≠0]

S_1 = \sum_{i=1}^{N}\sum_{j=1}^{N} s_{ij}\cdot [c_{ij}\neq 0], S_2 = \sum_{i=1}^{N} l_{ic_{ii}}\cdot [c_{ii}\neq 0]

其中[⋅][\cdot]是Iverson notation，当true时，取值为1。

我们可以得到S1+S2=∑Ni=1∑Nj=1Sij(cij)S_1+S_2 = \sum_{i=1}^{N}\sum_{j=1}^{N}S_{ij}(c_{ij})。同时必须满足下面3个约束:
1) 每个数据点i必须只能分配给一个代表点

2) 代表点一致性约束
如果数据点i选择了数据点j作为代表点，那么数据点j本身必须是代表点。

3) superexemplar一致性约束
如果数据点i选择了数据点k作为superexemplar, 那么k必须同时选择自己作为superexemplar。

所以MEAP的目标是最大化下面的目标函数：

上图描述的是一个多层的结构，S1评估的是within-subcluster 的compactness，高层的ψ2\psi_2描述的是exemplar 和superexemplar之间的关系。
根据single-exemplar的理论，最大化within-cluster 相似度会自动最大化between-cluster separation.
从最大化margin clustering的角度，MEAP比AP要好，见下图：

优化

MEAP的因子图如下图：

多代表点的模型是单代表点的模型的普及，对其进行优化是NP-hard的。因此，我们使用max-sum belief propagation，这是一个local-message-passing 算法，它会收敛到局部最大值。
（可以发现对角线的变量多连接了F的函数），也就是 superexemplar的一致性约束。

上图中圆形表示的变量节点，方形表示的是函数节点。
从变量到函数，将与这个变量连接的函数的信息求和（除了接收信息的函数以外）。

从函数到变量，包括除该变量外所有函数变量的maximization。

其中X=ne(f)X=ne(f)是函数ff的参数集。（或者我们可以理解成与该函数节点连接的所有变量节点）
这里我们可以发现对角线和非对角线的变量节点很不同，所以我们对它们分开进行讨论。

非对角线元素的Messages

左图有5种messages，与cijc_{ij}连接，i≠ji \neq j，如左下图：

非对角线元素

对角线元素

求解过程

非对角线变量

对于非对角线结点，m=cijm=c_{ij}

我们来讨论cjjc_{jj}，将这一项单独写出来,它的取值可以是１也可以是０
第一种情况:
如果是１，也就是m=1, 如果数据点i选择j作为代表点，那么j必须是一个代表点，这时Ej(c1j,...,cij)E_j(c_{1j},...,c_{ij}) 为0，而别的数据点i′i'可以选择ｊ作为代表点也可以不选，也就是ci′jc_{i'j}的值不定。
第二种情况:
m=0,也就是说ｉ没有选ｊ作为代表点，我们再讨论数据点ｊ的情况，如果数据点j本身是一个代表点，别的数据点i′i'要么选它作为数据点，要么不选，就如第一种情况的结论，但是如果ｊ本身不选自己，那么逆否命题成立我们可以推知，一定有cjj=0=>ci′j=0c_{jj}=0 => c_{i'j}=0,　所以可以得到∑i′:i′≠iρi′j(0)\sum_{i':i' \neq i} \rho_{i'j}(0)

对于ηij\eta_{ij},同样是从函数结点到变量结点。

ηij(m)=μIi→cij(m)=maxcij′:j′≠j[Ii(ci1,...,ciN)+∑j′′:j′≠jβij′′(cij′′)]

\eta_{ij}(m)=\mu_{I_i \rightarrow c_{ij}}(m)=max_{c_{ij' : j' \neq j}}[I_i(c_{i1},...,c_{iN})+\sum_{j'':j'\neq j}\beta_{ij''}(c_{ij''})]

第一种情况，cij=1c_{ij}=1,　由于唯一性约束，也就是它只能有一个代表点，那么β\beta参数cij′c_{ij'}必须是０。
第二种情况，如果是０，那么假设它的代表点是j′j',再讨论这个j′j'是不是ｉ本身，除去j′j'的代表点的所属权都必须是０。
(我一开始有点疑惑，因为大前提已经是i≠ji\neq j,后来我反应过来这只是对cijc_{ij}的讨论，而公式中覆盖的条件是全面的。)

对角线变量

对于对角线上的变量，如果它们本身是代表点，那么它们的取值为{1,...,N}\{1,...,N\},　如果它们不是代表点，那么取值为０。

(1)
如果ｉ是代表点，也就是cii=m,m∈{1,...,N}c_{ii}=m, m \in \{1,...,N\}那么别的数据点i′i'可以选择它或者不选它。
此时Ei(c1i,...,cii=m,...,cNi)=0E_i(c_{1i},...,c_{ii}=m,...,c_{Ni})=0
(2)
如果ｉ本身不是代表点，那么ci′i=0c_{i'i}=0

(1)
如果ｉ本身是一个代表点，也就是它选自己当代表点，所以它不能再选别人当代表点，由唯一性约束，那么别的cii′=0c_{ii'}=0
(2)
如果ｉ不是代表点，假如它选了i’,那么它不能选别的数据点作为代表点，所以可得

maxi′:i′≠i[βii′(1)+∑i′′:i′′∉{i,i′}βii′′(0)]

max_{i':i' \neq i}[\beta_{ii'}(1)+\sum_{i'':i''\notin \{i,i'\}}\beta_{ii''}(0)]

对于γ\gamma的讨论比较复杂，
(1) 如果cii=k=ic_{ii}=k=i，也就是说一个代表点ｋ，选了自己作为自己的superexemplar.
(2) 如果cii≠k=ic_{ii} \neq k = i
可以推出ckk≠kc_{kk} \neq k, 也就是说数据点ｋ并没有选自己作为superexemplar，那么别的数据点i′i'也不能选自己ｋ作为superexemplar。
(3) cii=k≠ic_{ii}=k \neq i
如果ｉ为代表点，选择ｋ作为超级代表点(superexemplar)，同时k≠ik\neq i，那么ｋ本身必须也是一个代表点同时选择自己作为超级代表点,也就是ckk=kc_{kk}=k,而别的除了i和ｋ的数据点并没有限制。
(4) cii≠k≠ic_{ii}\neq k \neq i
i没有选择ｋ作为超级代表点，我们可以对ｋ的情况进行讨论，如果它本身不是超级代表点，也就是ckk≠kc_{kk}\neq k,
这一项单独写出来可以得到

如果ｋ本身是超级代表点，ckk=kc_{kk}=k,单独写出来是：

最后的结果是两项的较大值。

Message summary与代码解读

S: similiarty matrix，n*n的矩阵
L: linkage matrix, n*n的矩阵
Rhoij : ρ~ij\tilde\rho_{ij} Rhoim: ρ~mi\tilde\rho_i^m
Alphaij: α~ij\tilde \alpha_{ij} (N*N) Alphai: α~ii\tilde\alpha_{ii}(N*1)
Betaij : β~ij\tilde \beta_{ij} Betaim : β~mi\tilde\beta_{i}^m
Etaij: η~ij\tilde \eta_{ij} (N*N) Etaii : η~ii\tilde\eta_{ii}(N*1)
Gammaik: γ~ik\tilde\gamma_{ik}
Phiik : ϕ~ik\tilde \phi_{ik}

Rhoij=zeros(N,N); Rhoim=zeros(N,N);
Alphaij=zeros(N,N); Alphai=zeros(N,1);
Betaij=zeros(N,N); Betaim=zeros(N,N);
Etaij=zeros(N,N); Etai=zeros(N,1);
Gammaik=zeros(N,N); Phiik=zeros(N,N);

更新ρ\rho

 %% rhoOldRhoij=Rhoij;Rhoij=S+Etaij;Rhoij=(1-lambda)*Rhoij+lambda*OldRhoij; %rho_ijOldRhoim=Rhoim;Rhoim=repmat(diag(S)+Etai,[1,N])+L+Gammaik;Rhoim=(1-lambda)*Rhoim+lambda*OldRhoim;

更新α\alpha

 OldAlphaij=Alphaij;OldAlphai=Alphai;Rp=max(Rhoij,0);for k=1:N, Rp(k,k)=max(Rhoim(k,:)); end;A=repmat(sum(Rp,1),[N,1])-Rp;dA=diag(A); Alphaij=min(A,0);for k=1:N, Alphai(k)=dA(k); end;Alphaij=(1-lambda)*Alphaij+lambda*OldAlphaij;Alphai=(1-lambda)*Alphai+lambda*OldAlphai;

更新β\beta

 Betaij=S+Alphaij;Betaim=repmat(diag(S)+Alphai,[1,N])+L+Gammaik;

更新η\eta

B=Betaij; [Y,I]=max(B,[],2);
for i=1:N, B(i,I(i))=-inf; end;%最大值赋值为负无穷
[Y2,I2]=max(B,[],2);% Y2代表次大值
R=repmat(Y,[1,N]);

当j的位置刚好对应最大值时，由于j′≠jj'\neq j所以取次大值
这一步得到

for i=1:N,R(i,I(i))=Y2(i);R(i,i)=-inf;
end;

接下来计算：

  T=zeros(N,N);for i=1:N, T(i,:)=max(Betaim(i,:)); end;

接下来计算：

  Etaij=-max(R,T);

接下来计算

  T=Betaij;for i=1:N, T(i,i)=-inf; end;%除去i'不等于ｉ的情况Etai=-max(T,[],2);

更新γ\gamma

这一步计算

    OldGammaik=Gammaik;RPhiik=max(Phiik,0); for k=1:N, RPhiik(k,k)=Phiik(k,k); end;Gammaik=repmat(sum(RPhiik,1),[N,1])-RPhiik;

对于对角线元素，通过上面的计算已经求出，也就是Gammaik的对角线,所不同的是对非对角线，还有一个min函数的截断。所以先存储对角线元素，随后对矩阵进行min函数的截断，再单独赋值。

    dGammaik=diag(Gammaik); Gammaik=min(Gammaik,0); for k=1:N, Gammaik(k,k)=dGammaik(k); end;Gammaik=(1-lambda)*Gammaik+lambda*OldGammaik;

更新ϕ\phi

计算

    OldPhiik=Phiik;LG=L+Gammaik; [Y,I]=max(LG,[],2);for i=1:N, LG(i,I(i))=-inf; end; [Y2,I2]=max(LG,[],2);R=repmat(Y,[1,N]);for i=1:N, R(i,I(i))=Y2(i); end;

计算sii+α~ii+η~iis_{ii}+\tilde\alpha_{ii}+\tilde\eta_{ii}

    SAE=repmat(diag(S)+Alphai+Etai,[1,N]);Phiik=L-max(R,-SAE);Phiik=(1-lambda)*Phiik+lambda*OldPhiik;

Assignment matrix

计算cijc_{ij}的所有输入信息的和，随后求解使得c^ij\hat c_{ij}最大的值。

    C=zeros(N,N);C(Alphaij+Rhoij>=0)=1;%非对角线[a,b]=max(Rhoim,[],2);for i=1:N, if a(i)+Alphai(i)>=0, C(i,i)=b(i); else C(i,i)=0; end; end;%对角线if isequal(OldC,C), stayiter=stayiter+1; else stayiter=1; OldC=C; end;if stayiter>convits% further check for validationexemplar_idx=zeros(N,1);validflag=1;for i=1:Na=find(C(i,:)~=0);if length(a)~=1, validflag=0; break; else exemplar_idx(i)=a(1); end; % points to exemplarsendif validflag==1a=find(diag(C)~=0);if length(a)<=1%只有一个代表点validflag=0;elsesuperexemplar_idx=zeros(length(a),2);superexemplar_idx(:,1)=a;superexemplar_idx(:,2)=diag(C(a,a));sa=unique(superexemplar_idx(:,2));if ~isequal(diag(C(sa,sa)),sa), validflag=0; end; % exemplars to super-exemplarsendendif validflag==0 % invalid results! continue iterating!stayiter=1; OldC=zeros(N,N);elsedisp(['iteration number ' num2str(iter)]);break;endend
end

计算net similarity

T=diag(L);
netSim=sum(S(C~=0))+sum(T(diag(C)~=0));