作者：桂。
時間：2017-04-14 06:22:26
鏈接：http://www.cnblogs.com/xingshansi/p/6685811.html
聲明：歡迎被轉載，不過記得注明出處哦~
前言
之前梳理了一下非負矩陣分解（Nonnegative matrix factorization, NMF），主要有：
3）拉格朗日乘子法求解NMF（將含限定NMF的求解 一般化）
譜聚類可以參考之前的文章：
1）拉普拉斯矩陣（Laplace Matrix）與瑞利熵（Rayleigh quotient）
2）譜聚類（Spectral clustering）（1）：RatioCut
3）譜聚類（Spectral clustering）（2）:NCut
總感覺NMF跟聚類有聯系，這里試着從聚類角度分析一下非負矩陣分解，主要包括：
1）Kmeans與譜聚類
2）對稱非負矩陣分解（symmetric NMF，SyNMF）;
3）非對稱非負矩陣分解;
內容為自己的學習總結，如果有不對的地方，還請幫忙指出。文中多有借鑒他人的地方，最后一並給出鏈接。
一、Kmeans與譜聚類
A-Kmeans定義
，准則函數為：
可以重寫為：
定義h：
$n_k$為第k類樣本的個數，則准則函數變為：
從而Kmeans的優化問題，等價於：
B-Kmeans與譜聚類（Spectral clustering）的聯系
上文給出了h的定義：
回顧譜聚類中RatioCut定義h的思路：
回顧RatioCut求解問題：
可以看出求解的思路完全一致，不同的是RatioCut是拉普拉斯矩陣L，而Kmeans是矩陣$X^TX$。正因為這點不同，Kmeans在利用譜聚類的思路求解時，略有差別。利用Kmeans的思想求Kmeans，聽着這不是欠揍嗎？ 這里只是為了分析譜聚類一般方法（L）與譜聚類對應的Kmeans（$X^TX$）二者的不同。
再次寫出RatioCut的步驟：
步驟一：求解拉普拉斯矩陣L
步驟二：對L進行特征值分解，並取K個最小特征值對應的特征向量（K為類別數目）
步驟三：將求解的K個特征向量（並分別歸一化），構成新的矩陣，對該矩陣進行kmeans處理
kmeans得到的類別標簽，就是原數據的類別標簽，至此完成RatioCut聚類。
對應Kmeans呢？是不是把$X^TX$換成L就等價於對原數據Kmeans?
步驟一：求解$X^TX$
步驟二：對$X^TX$進行特征值分解，並取K個最大特征值對應的特征向量（K為類別數目）
步驟三：將求解的K個特征向量（並分別歸一化），構成新的矩陣，對該矩陣進行kmeans處理
kmeans得到的類別標簽，就是原數據的類別標簽，至此完成RatioCut聚類。
對應測試code:
L = X'*X;% matrix %%Step2:Eigenvalues decomposition K = 3; [Qini,V] = eig(L); %%Step3:New matrix Q [~,pos] = sort(diag(V),'ascend'); Q = Qini(:,pos(1:K)); Q = Q./repmat(sqrt(diag(Q'*Q)'),N,1); [idx,ctrs] = kmeans(Q,K);
我們可以利用數據測試一下：
再來看看直接對原數據Kmeans的結果：
利用譜聚類的思想求解，得出了錯誤的分類結果，而直接Kmeans效果是理想的。
原因何在？問題就出在矩陣L上。回顧之前提到的拉普拉斯矩陣L特性：
它表示的是不同數據點特征的差距，而$X^TX$呢？
更像是一種余弦距離的測度，普通的Kmeans自然不能很好解決聚類。其實對$X^TX$利用譜聚類，特征向量子空間是嚴謹的，以三類為例，將特征子空間投影到二維，如下圖中間所示，很容易看出子空間特征是分成三類的，但聚類之后呢？如右圖所示，就出現了判別錯誤。這也是為什么上面兩種結果不一致。
總而言之：Kmeans是與譜聚類的思想一致，但由於中間矩陣不同，二者思路略有差異。
C-Kernel Kmeans與譜聚類
這里就再啰嗦一下，數據分類效果不理想，映射到高維呢？也就是核函數（Kernel function）的思想。
重新寫出譜聚類框架下的Kmeans：
對X進行映射：
，得出核函數下的Kmeans(Kernel Kmeans)
重新給出code(以kernel 取gaussian為例):
%Kmeans
sigma2 = .02;
L = exp(-dist(X,X').^2*sigma2);
[Qini,V] = eig(L);
%%Step3:New matrix Q
[~,pos] = sort(diag(V),'descend');
Q = Qini(:,pos(1:K));
for i =1:K
Q(:,i) = Q(:,i)-min(Q(:,i));
end
Q = Q./repmat(sqrt(diag(Q'*Q)'),K*N,1);
[idx,ctrs] = kmeans(Q,K);
對應分類結果：
這個時候分類就理想了。
二、對稱非負矩陣分解（SyNMF）
A-原理介紹
Kmeans與RatioCut的理論框架是統一的，其實准則函數等價為：
現在將$H^TH = I$的約束去掉，泛化后的求解問題為：
這就是對稱非負矩陣分解（SyNMF）的思路。
B-算法求解
求解思路還是利用拉格朗日乘子+KKT，不再細說，給出結果：
泛化后：
得到H之后，如何實現數據的label判別？
可以看出得到的H分為三類，對應三種label, 即可實現數據分離。
給出SymNMF與譜聚類的對比：
對應的code可以點擊這里。
給出一個測試結果圖，測試數據為三類：
三、非對稱非負矩陣分解
SymNMF是Spectral clustering的泛化推廣。
上文分析的是對稱的譜聚類問題：
- Spectral clustering
- SymNMF
同樣的方式，分析非對稱的譜聚類問題：
- Spectral clustering
該問題可以轉化為：
同樣的，對於：
進行泛化：
這就是NMF的准則函數，即：
- NMF
現在來總結一下：
關系是不是一目了然了？如何添加更多約束項呢？比如希望H矩陣盡可能系數等等，就在上面這幾類問題的准則函數后添加約束，轉化成對偶問題求解即可。
題外話
A-非負矩陣NMF實現數據聚類
分析了這么多，已經解開了之前的困惑：譜聚類與NMF之間的聯系。
回顧上面分析Kmeans提到的三類數據聚類問題，對$XX^T$進行NMF處理：
更直觀地，將數據放在對應維度觀察：
可以看到，由於矩陣對稱，即figure對稱，對W/H聚類，都可以得到三類標簽，從而實現數據的聚類。對於非對稱呢？自然想到：如果W對應樣本數量維度，則對W進行聚類，如果H對應樣本數量維度，就對H進行聚類，同樣可以實現數據聚類。
給出一個利用（Kernel kmeans）聚類以及利用NMF聚類的對比示例。
利用$X^TX$是kmeans的思路，利用$L$是RatioCut思路，事實上，泛化之后，直接利用$X^TX$也可以利用是NMF的思路。
clc;clear all;close all;
N = 100;
K = 3;
X = [randn(N,2)+ones(N,2);...
randn(N,2)+3*ones(N,2);...
randn(N,2)-4*ones(N,2)];
posran = randperm(K*N);
X = X-min(min(X));
X = X(posran,:);
%%方法一 基於Data: 利用 XX' + 核心函數
sigma2 = 0.02;
L = exp(-dist(X,X')*sigma2);
%%Step2:Eigenvalues decomposition
[Qini,V] = eig(L);
%%Step3:New matrix Q
[~,pos] = sort(diag(V),'descend');
Q = Qini(:,pos(1:K));
for i =1:K
Q(:,i) = Q(:,i)-min(Q(:,i));
end
Q = Q./repmat(sqrt(diag(Q'*Q)'),K*N,1);
[idx,ctrs] = kmeans(Q,K);
figure
subplot 221
plot(X(:,1),X(:,2),'.');
title('原數據');grid on;
subplot 222
plot3(Q(:,1),Q(:,2),Q(:,3),'b.');
title('特征');grid on;
subplot 223
plot3(Q(idx==1,1),Q(idx==1,2),Q(idx==1,3),'b.');hold on;
plot3(Q(idx==2,1),Q(idx==2,2),Q(idx==2,3),'r.');hold on;
plot3(Q(idx==3,1),Q(idx==3,2),Q(idx==3,3),'g.');hold on;
title('特征聚類');grid on;
subplot 224
plot(X(idx==1,1),X(idx==1,2),'b.');hold on;
plot(X(idx==2,1),X(idx==2,2),'r.');hold on;
plot(X(idx==3,1),X(idx==3,2),'g.');hold on;
title('聚類結果');grid on;
%%方法2 基於NMF: 對X*X'聚類
%即NMF
Iter = 500;
K = 3;
[W,H] = nmf(X*X',K,Iter);
[idx,ctrs] = kmeans(W,K,'distance','cityblock');
figure
subplot 221
plot(X(:,1),X(:,2),'.');
title('原數據');grid on;
subplot 222
plot3(W(:,1),W(:,2),W(:,3),'b.');
title('W矩陣');grid on;
subplot 223
plot3(W(idx==1,1),W(idx==1,2),W(idx==1,3),'b.');hold on;
plot3(W(idx==2,1),W(idx==2,2),W(idx==2,3),'r.');hold on;
plot3(W(idx==3,1),W(idx==3,2),W(idx==3,3),'g.');hold on;
title('W矩陣聚類');grid on;
subplot 224
plot(X(idx==1,1),X(idx==1,2),'b.');hold on;
plot(X(idx==2,1),X(idx==2,2),'r.');hold on;
plot(X(idx==3,1),X(idx==3,2),'g.');hold on;
title('聚類結果');grid on;
這里的代碼只是為了說明理論問題，代碼本身不足以支持應用。給出對應結果圖：
基於$X^TX$的Kmeans:
基於$X^TX$的NMF:
事實上，直接對$X$進行NMF聚類，只要將$X^TX$替換為$X$即可，對應結果圖：
這里可以理解為：H是對應的字典信息，W是線性變換。
B-圖的聚類（不再是數據的聚類）
上文分析的是對數據點進行分類，如果直接對鄰接矩陣、拉普拉斯矩陣，也就是圖片信息進行分類呢？（可能有點繞，但數據點聚類，並不代表圖就是聚類，圖明顯可以聚類，對應的數據點也未必可以聚類）。
圖對應的矩陣，如拉普拉斯矩陣、鄰接矩陣等等，說到底都是figure的表達，如上圖所示，因此上文分析的矩陣W/B，可以用圖片的信息替代。
對於圖的聚類，一種思路是將圖中像素點讀取，轉化成數據格式，再進行聚類。
Data角度就是前面分析的種種類型，Figure角度呢？
給出NMF對圖片的處理結果（分四類）：
對應code:
img = imread('2.png');
V = 1-im2bw(rgb2gray(img));
figure
K = 4;
Iter = 500;
[W,H] = nmf(V, K, Iter);
subplot 331
mesh(V);
title('原數據')
for i =1:K
subplot(3,3,i+1);
mesh(W(:,i)*H(i,:));
title(['分解',num2str(i)])
end
NMF對於Figure是一種線性表達的思路，可以分離的基礎是數據不共線（橫/縱），再給出之前音樂分離的語譜圖：
之所以可以分離正是因為音樂頻譜的周期性，所以對於非周期的說話人，直接應用NMF應該是無效的，對應的時域波形：
參考：
- On the Equivalence of Nonnegative Matrix Factorization and Spectral Clustering
- Symmetric Nonnegative Matrix Factorization for Graph Clustering