谱聚类的实现

开始用accumulate做加和，结果发现laplacian矩阵求特征值后最小的特征值不是0，这是有问题的，聚类的准确率不是很高，原因也找不到。后来一步步查，发现diag矩阵有问题，最终查到accumulate是int相加，不准确，修改后准确率增加。

// spectral-cluster.cpp : 定义控制台应用程序的入口点。//#include "stdafx.h"#include <iostream> #include<vector>#include<set>#include <numeric>#include <cv.h>#include <highgui.h>/*********************author Marshall**************************//************************2015.10.20****************************//*********************version 1.0******************************/using namespace std;class spectral;class kmeans{friend class spectral;private:double*dataset;unsigned int k;unsigned int dim;unsigned int numofdata;typedef vector<double> Centroid;vector<Centroid> center;vector<set<int>>cluster_ID;vector<Centroid>new_center;vector<set<int>>new_cluster_ID;double threshold;int iter;private:void init();void assign();double distance(Centroid cen, int k2);void split(vector<set<int>>&clusters, int kk);void update_centers();bool isfinish();void show_result();public:kmeans(double*data, unsigned int datanums, const int dim, const int numofcluster):dataset(data), numofdata(datanums), dim(dim), k(numofcluster){threshold = 0.0000001;}void apply();};//template <typename T>void kmeans::init(){center.resize(k);set<int>bb;for (int i = 0; i < k; i++){int id = double(rand()) / double(RAND_MAX + 1.0)*numofdata;while (bb.find(id) != bb.end()){id = double(rand()) / double(RAND_MAX + 1.0)*numofdata;}bb.insert(id);center[i].resize(dim);for (int j = 0; j < dim; j++)center[i][j] = dataset[id*dim + j];}}bool kmeans::isfinish(){double error = 0;for (int i = 0; i < k; i++){for (int j = 0; j < dim; j++)error += pow(center[i][j] - new_center[i][j], 2);}return error < threshold ? true : false;}void kmeans::assign(){for (int j = 0; j < numofdata; j++){double mindis = 10000000;int belongto = -1;for (int i = 0; i < k; i++){double dis = distance(center[i], j);if (dis < mindis){mindis = dis;belongto = i;}}new_cluster_ID[belongto].insert(j);}for (int i = 0; i < k; i++){if (new_cluster_ID[i].empty()){split(new_cluster_ID, i);}}}double kmeans::distance(Centroid cen, int k2){double dis = 0;for (int i = 0; i < dim; i++)dis += pow(cen[i] - dataset[k2*dim + i], 2);return sqrt(dis);}void kmeans::split(vector<set<int>>&clusters, int kk){int maxsize = 0;int th = -1;for (int i = 0; i < k; i++){if (clusters[i].size() > maxsize){maxsize = clusters[i].size();th = i;}}#define DELTA 3vector<double>tpc1, tpc2;tpc1.resize(dim);tpc2.resize(dim);for (int i = 0; i < dim; i++){tpc2[i] = center[th][i] - DELTA;tpc1[i] = center[th][i] + DELTA;}for (set<int>::iterator it = clusters[th].begin(); it != clusters[th].end(); it++){double d1 = distance(tpc1, *it);double d2 = distance(tpc2, *it);if (d2 < d1){clusters[kk].insert(*it);}}_ASSERTE(!clusters[kk].empty());for (set<int>::iterator it = clusters[kk].begin(); it != clusters[kk].end(); it++)clusters[th].erase(*it);}void kmeans::update_centers(){for (int i = 0; i < k; i++){Centroid temp;temp.resize(dim);for (set<int>::iterator j = new_cluster_ID[i].begin(); j != new_cluster_ID[i].end(); j++){for (int m = 0; m < dim; m++)temp[m] += dataset[(*j)*dim + m];}for (int m = 0; m < dim; m++)temp[m] /= new_cluster_ID[i].size();new_center[i] = temp;}}void kmeans::apply(){init();new_center.resize(k);new_cluster_ID.resize(k);assign();update_centers();iter = 0;while (!isfinish()){center = new_center;cluster_ID = new_cluster_ID;new_center.clear();new_center.resize(k);new_cluster_ID.clear();new_cluster_ID.resize(k);assign();update_centers();iter++;}//new_cluster_ID.clear();}class spectral{private:unsigned int dim;unsigned int N;double**normalized_data;double*laplacian_matrix;double*diag_matrix;double*similarity_matrix;bool is_normalized;double sigma2;unsigned int k_nearest;unsigned int K;unsigned int numofcluster;private:void create_similarity_matrix();void create_diag_matrix();void create_laplacian_matrix();void eigen_of_laplacian(double*&newdata);//void kmeans();void normalize_data(double*&data, int datanums, int dim);double euclidean(const double*data1, const double*data2);bool is_knearest(const double*distance_matrix, const int kk, const int mm);public:spectral(double**data, const int N, const int dim, const int numofcluster);void spectral_cluster(int*pClusters);~spectral();};spectral::~spectral(){}spectral::spectral(double**data, const int N, const int dim, const int numofcluster):N(N), dim(dim), numofcluster(numofcluster){K = 4;k_nearest = 15;sigma2 = 600;normalized_data = data;is_normalized = false;}void spectral::normalize_data(double*&data, int datanums, int dim){double*mins = new double[dim];memset(mins, 0x3f3f3f3f, sizeof(double)*dim);double*maxs = new double[dim];memset(maxs, 0, sizeof(double)*dim);for (int i = 0; i < datanums; i++)for (int j = 0; j < dim; j++){if (data[i*dim + j] > maxs[j]){maxs[j] = data[i*dim + j];}if (data[i*dim + j] < mins[j])mins[j] = data[i*dim + j];}for (int i = 0; i < datanums; i++)for (int j = 0; j < dim; j++)data[i*dim + j] = 100 * (data[i*dim + j] - mins[j]) / (maxs[j] - mins[j]);delete[]maxs, mins;}double spectral::euclidean(const double*data1, const double*data2){double dis = 0;for (int i = 0; i < dim; i++)dis += pow(data1[i] - data2[i], 2);return dis;}bool spectral::is_knearest(const double*distance_matrix,const int kk, const int mm){double dis = distance_matrix[kk*N + mm];int count = 0;//还要考虑distance_matrix[kk*N + kk]=0;//因为k_nearest相比N很小（大多数情况），判断false好for (int i = 0; i < N; i++)if (i != kk&&distance_matrix[kk*N + i] < dis){count++;if (count > k_nearest)return false;}return true;}void spectral::create_similarity_matrix(){double*distance_matrix = new double[N*N];memset(distance_matrix, 0, N*N*sizeof(double));for (int i = 0; i < N - 1; i++)for (int j = i + 1; j < N; j++){distance_matrix[i*N + j] = euclidean(normalized_data[i], normalized_data[j]);distance_matrix[j*N + i] = distance_matrix[i*N + j];}if (similarity_matrix)delete[]similarity_matrix;similarity_matrix = new double[N*N];memset(similarity_matrix, 0, N*N*sizeof(double));//这里使用对称k最近邻for (int i = 0; i < N - 1; i++){for (int j = i + 1; j < N; j++){if (is_knearest(distance_matrix, i, j)|| is_knearest(distance_matrix, j, i)){similarity_matrix[i*N + j] = similarity_matrix[j*N + i]= 100 * exp(-distance_matrix[i*N + j] / sigma2);}//cout << similarity_matrix[i*N + j] << "     ";}//cout << endl << endl << endl << endl;}/*for (int i = 0; i < N; i++){for (int j = 0; j < N; j++){cout << similarity_matrix[i*N + j] << "     ";}cout << endl << endl << endl << endl;}*/delete[]distance_matrix;}void spectral::create_diag_matrix(){if (diag_matrix)delete[]diag_matrix;diag_matrix = new double[N*N];memset(diag_matrix, 0, N*N*sizeof(double));for (int i = 0; i < N; i++){//accumulate是int相加，不准确，修改后准确率增加//diag_matrix[i*N + i] = accumulate(similarity_matrix + i*N, similarity_matrix + i*N + N, 0);double sum = 0;for (int j = 0; j < N; j++)sum += similarity_matrix[i*N + j];diag_matrix[i*N + i] = sum;//cout << "diag" << i << "   " << diag_matrix[i*N + i] << endl;}}void spectral::create_laplacian_matrix(){if (laplacian_matrix)delete[]laplacian_matrix;laplacian_matrix = new double[N*N];for (int i = 0; i < N; i++)for (int j = i; j < N; j++){laplacian_matrix[i*N + j] = laplacian_matrix[j*N + i]= diag_matrix[i*N + j] - similarity_matrix[i*N + j];}if (is_normalized){for (int i = 0; i < N; i++)for (int j = 0; j < N; j++){int temp = i == j ? 1 : 0;laplacian_matrix[i*N + j] = temp - 1.0 / sqrt(diag_matrix[j*N + j]+ 1.0e-13)*laplacian_matrix[i*N + j] * 1.0 / sqrt(diag_matrix[i*N + i] + 1.0e-13);}}delete[]diag_matrix, similarity_matrix;}void spectral::eigen_of_laplacian(double*&newdata){//CvMat Ma = cvMat(N, N, CV_64FC1, laplacian_matrix);CvMat *pSrcMatrix = cvCreateMat(N, N, CV_64FC1);// E = 对应的特征向量 (每行)CvMat* E = cvCreateMat(N, N, CV_64FC1);CvMat* l = cvCreateMat(N, 1, CV_64FC1);for (int i = 0; i < N; i++){cvmSet(l, i, 0, 0);for (int j = 0; j < N; j++){cvmSet(pSrcMatrix, i, j, laplacian_matrix[i*N + j]);cvmSet(E, i, j, 0);}}cvEigenVV(pSrcMatrix, E, l);   // l = A的特征值 (降序排列)for (int i = 0; i < N; i++)cout << cvmGet(l, i, 0) << "      ";//半正定，最小的特征值应该是0，应该在此处验证for (int i = 0; i < N; i++){//cout << cvmGet(E, N - 1 , i);for (int j = 0; j < K; j++){newdata[i*K + j] = cvmGet(E, N - 2 - j, i);}}delete[]laplacian_matrix;// 释放矩阵所占内存空间cvReleaseMat(&pSrcMatrix);cvReleaseMat(&E);cvReleaseMat(&l);}void spectral::spectral_cluster(int*pClusters){create_similarity_matrix();create_diag_matrix();create_laplacian_matrix();double*newdata = new double[N*K];eigen_of_laplacian(newdata);//normalize_data(newdata, this->N, dim);kmeans km(newdata, N, K, this->numofcluster);km.apply();delete[]newdata;newdata = NULL;int count = 0;for (int i = 0; i < km.new_cluster_ID.size(); i++)for (set<int>::iterator it = km.new_cluster_ID[i].begin(); it != km.new_cluster_ID[i].end(); it++){pClusters[*it] = i;count++;}_ASSERTE(count == N);}int _tmain(int argc, _TCHAR* argv[]){#define MAX_CLUSTERS 5  CvScalar color_tab[MAX_CLUSTERS];IplImage* img = cvCreateImage(cvSize(500, 500), IPL_DEPTH_8U, 3);CvRNG rng = cvRNG(cvGetTickCount());CvPoint ipt;color_tab[0] = CV_RGB(255, 0, 0);color_tab[1] = CV_RGB(0, 255, 0);color_tab[2] = CV_RGB(100, 100, 255);color_tab[3] = CV_RGB(255, 0, 255);color_tab[4] = CV_RGB(255, 255, 0);int clusterNum = 4;int sampleNum = 500;int feaNum = 2;CvMat* sampleMat = cvCreateMat(sampleNum, 1, CV_32FC2);/* generate random sample from multigaussian distribution *///产生高斯随机数  for (int k = 0; k < clusterNum; k++){CvPoint center;int topIndex = k*sampleNum / clusterNum;int bottomIndex = (k == clusterNum - 1 ? sampleNum : (k + 1)*sampleNum / clusterNum);CvMat* localMat = cvCreateMatHeader(bottomIndex - topIndex, 1, CV_32FC2);center.x = cvRandInt(&rng) % img->width;center.y = cvRandInt(&rng) % img->height;cvGetRows(sampleMat, localMat, topIndex, bottomIndex, 1);  //此函数不会给localMat分配空间,不包含bottomIndex  cvRandArr(&rng, localMat, CV_RAND_NORMAL,cvScalar(center.x, center.y, 0, 0),cvScalar(img->width*0.1, img->height*0.1, 0, 0));}// shuffle samples 混乱数据  for (int i = 0; i < sampleNum / 2; i++){CvPoint2D32f* pt1 = (CvPoint2D32f*)sampleMat->data.fl + cvRandInt(&rng) % sampleNum;CvPoint2D32f* pt2 = (CvPoint2D32f*)sampleMat->data.fl + cvRandInt(&rng) % sampleNum;CvPoint2D32f temp;CV_SWAP(*pt1, *pt2, temp);}double** pSamples = new double*[sampleNum];for (int i = 0; i < sampleNum; i++)pSamples[i] = new double[feaNum];int* pClusters = new int[sampleNum];for (int i = 0; i < sampleNum; i++){//feaNum=2  pSamples[i][0] = (cvGet2D(sampleMat, i, 0)).val[0];//cout << pSamples[i][0] << "　　";pSamples[i][1] = (cvGet2D(sampleMat, i, 0)).val[1];//cout << pSamples[i][1] << endl;}cvReleaseMat(&sampleMat);cvZero(img);for (int i = 0; i < sampleNum; i++){//feaNum=2  ipt.x = pSamples[i][0];ipt.y = pSamples[i][1];cvCircle(img, ipt, 1, cvScalar(255, 255, 255), CV_FILLED, CV_AA, 0);}cvSaveImage("origin.jpg", img);//执行谱聚类spectral sp(pSamples, sampleNum, feaNum, clusterNum);sp.spectral_cluster(pClusters);cvZero(img);for (int i = 0; i < sampleNum; i++){//feaNum=2  int cluster_idx = pClusters[i];ipt.x = pSamples[i][0];ipt.y = pSamples[i][1];cvCircle(img, ipt, 1, color_tab[cluster_idx], CV_FILLED, CV_AA, 0);}delete[] pClusters;pClusters = NULL;delete[] pSamples;pSamples = NULL;cvNamedWindow("clusters");cvShowImage("clusters", img);cvWaitKey(0);cvSaveImage("clusters.jpg", img);cvReleaseImage(&img);system("pause");cvDestroyWindow("clusters");return 0;}

下面是聚类结果。

     

嫌c++冗长的看这个python的

#!/usr/bin/python# copyright (c) 2008 Feng Zhu, Yong Sunimport heapqfrom functools import partialfrom numpy import *from scipy.linalg import *from scipy.cluster.vq import *import pylabdef line_samples ():    vecs = random.rand (120, 2)    vecs [:,0] *= 3    vecs [0:40,1] = 1    vecs [40:80,1] = 2    vecs [80:120,1] = 3    return vecsdef gaussian_simfunc (v1, v2, sigma=1):    tee = (-norm(v1-v2)**2)/(2*(sigma**2))    return exp (tee)def construct_W (vecs, simfunc=gaussian_simfunc):    n = len (vecs)    W = zeros ((n, n))    for i in xrange(n):        for j in xrange(i,n):            W[i,j] = W[j,i] = simfunc (vecs[i], vecs[j])    return Wdef knn (W, k, mutual=False):    n = W.shape[0]    assert (k>0 and k<(n-1))    for i in xrange(n):        thr = heapq.nlargest(k+1, W[i])[-1]        for j in xrange(n):            if W[i,j] < thr:                W[i,j] = -W[i,j]    for i in xrange(n):        for j in xrange(i, n):            if W[i,j] + W[j,i] < 0:                W[i,j] = W[j,i] = 0            elif W[i,j] + W[j,i] == 0:                W[i,j] = W[j,i] = 0 if mutual else abs(W[i,j])vecs = line_samples()W = construct_W (vecs, simfunc=partial(gaussian_simfunc, sigma=2))knn (W, 10)D = diag([reduce(lambda x,y:x+y, Wi) for Wi in W])L = D - Wevals, evcts = eig(L,D)vals = dict (zip(evals, evcts.transpose()))keys = vals.keys()keys.sort()Y = array ([vals[k] for k in keys[:3]]).transpose()res,idx = kmeans2(Y, 3, minit='points')colors = [(1,2,3)[i] for i in idx]pylab.scatter(vecs[:,0],vecs[:,1],c=colors)pylab.show()

机器学习算法复习-谱聚类

以下转自 3月机器学习在线班第七课--聚类方法

相似度/距离计算方法总结

在我们正式讲解聚类算法之前，让我们来简单总结一下我们已有的距离计算方法。聚类，无法避免的就是需要计算距离，而不同计算距离方法的选择，则有可能导致差异较大的不同结果。

闵可夫斯基Minkowski距离： 

当时，Minkowski距离即为我们熟悉的欧式距离。

杰卡德相似系数(Jaccard): 

余弦相似度(cosine similarity): 

Pearson相似系数： 

相对熵(K-L距离): 

聚类的基本思想

在我们熟悉的机器学习算法和模型中，大部分情况下训练数据都是带有label（标签）的，这样的算法是有监督学习算法。当训练数据没有label（标签）时，便称之为无监督学习。聚类算法就是无监督学习中最常见的一种，给定一组数据，需要聚类算法去挖掘数据中的隐含信息。聚类算法的应用很广：顾客行为聚类，google新闻聚类等。

定义：给定一个有个对象的数据集，聚类将数据划分为个簇，而且这个划分满足两个条件：(1)每个簇至少包含一个对象；(2)每个对象属于且仅属于一个簇。

基本思想：对给定的，算法首先给出一个初始的划分方法，以后通过反复迭代的方法改变划分，使得每一次改进之后的划分方案都较前一次更好。

基础聚类方法--K-means算法

K-means算法，也被称为K-均值，是一种得到最广泛使用的聚类算法，它也可以成为其他聚类算法的基础。

基本思想： K-means算法首先随机地选择个对象，每个对象初始地代表一个簇的平均值或中心。对剩余的每个对象根据其与各个簇中心的距离（之前的各种距离的定义，在这里派上了用场，可以根据实际需求来选择不同的距离定义），将它赋给最近的簇。然后重新计算每个簇的平均值。这个过程不断重复，直到准则函数函数收敛（准则函数常常使用最小均方误差）。

算法流程:

设给定输入数据为,K-means算法如下： 
(1)选取初始的个聚类中心; 
(2)对每个样本数据而言，将其类别标签设为距离其最近的聚类中心的标签，即; 
(3)将每个聚类中心的值更新为与该类别中心相同类别的所有样本的平均值，即 

(4)重复(2)(3)步，直到聚类中心的变化()小于规定的阈值为止。

下图即为一个简单的示例（）: 

在K-means算法中，如何选择初始的聚类中心是一个普遍的问题(该算法是初值敏感的)。同时，K-means算法只能够保证收敛到一个局部极值，无法保证收敛到全局最优。一个较为简单的解决办法是随机初始化多次，以最优的聚类结果为最终结果。

在聚类结束后，如果一个中心没有得到任何样本(这是可能的)，那么需要去除这个中心点，或者重新初始化。

密度聚类方法

密度聚类方法的指导思想是：只要一个区域中的点的密度大于某个阈值，就把它加到与之相近的类簇中去。这类算法可以克服基于距离的算法只能够发现“类圆形”（凸）的聚类的缺点，可以发现任何形状的聚类，且对噪声数据不敏感。但是计算密度单元的计算复杂度大，需要建立空间索引来降低计算量。

两个主要的密度聚类方法：DBSCAN算法、密度最大值算法。

DBSCAN算法

DBSCAN（Density-Based Spatial Clustering of Applications with Noise，具有噪声的基于密度的聚类方法）是一种基于密度的空间聚类算法。该算法将具有足够密度的区域划分为簇，并在具有噪声的空间数据库中发现任意形状的簇，它将簇定义为密度相连的点的最大集合。 
可以想象一下人类聚居的过程，黄河（母亲河）边上首先会聚集更多的人，当人数足够多了之后，就会产生城市（类似于簇的概念）。可以想象一下城市的形状是不规则的，但城市是以空间上的人口密度作为划分基础的。

算法设计的基本定义： 
1.对象的-领域：在给定对象的半径内的区域。 
2.核心对象(core point)：对于给定的数目m(人为设定的标准)，如果一个对象的-领域至少包含了m个对象，那么称该对象为核心对象。 
3.直接密度可达(directly density-reachable)：给定一个对象集合,如果对象是在对象的-领域内，且是一个核心对象，那么我们说对象从对象出发是直接密度可达的。

Tips:以下图片帮助理解：如图=1cm，m=5，可以知道是一个核心对象，从对象出发到对象出发是直接密度可达的。 

4.密度可达(density-reachable)：如果存在一个对象链,令。假设对任意,有,且是从关于和m直接密度可达的，那么对象是从对象关于关于和m密度可达的。 
5.密度相连(density-connected)：如果对象集合中存在一个对象，使得对象p和q是从关于和m密度可达的，那么对象p和q是关于和m密度相连的。 
6.簇(cluster)：一个基于密度的簇是最大的密度相连对象的集合。 
7.噪声：不包含在任何簇中的对象我们称之为噪声。

Tips:以下图片帮助理解：如图m=3，红色部分的点都是核心对象，同样的，他们彼此之间是密度可达的。点B和点C不是核心对象，但是它们关于点A是密度可达的，于是它们属于同一个簇。点N是噪声点。 

通过以上定义，我们可以知道这些定义是一层层递进的，慢慢理解下来其实不难。我们也可以发现密度可达并不是双向关系的，而密度相连则是双向关系的。 
现在假设点是一个核心对象，那么通过不断搜索数据（这些数据可以是核心对象也可以不是核心对象）就可形成一个以为核心对象的簇。每一个簇至少包含一个核心对象，非核心对象可以是簇的一部分，它（非核心对象）构成的是簇的边缘(edge)。

DBSCAN算法流程 
DBSCAN需要两个参数：和m。该算法一开始选择任意一个未访问过的点p作为起始点，通过查看对象p的-领域中数据点的个数k，如果km,那么该点会被看成是异常点（值得注意的是该点有可能是另外一个核心对象的密度相连对象，到后来也可能是簇的一部分）；如果km，那么就会构造一个以点p为核心对象的簇（所有密度相连对象的集合）。在这个过程中还会有簇的合并的过程。 
不断重复以上过程，就可以将原始数据区分为多个簇以及噪声点。 
伪代码如下：

DBSCAN(D, eps, MinPts) {   C = 0   for each point P in dataset D {      if P is visited         continue next point      mark P as visited      NeighborPts = regionQuery(P, eps)      if sizeof(NeighborPts) < MinPts         mark P as NOISE      else {         C = next cluster         expandCluster(P, NeighborPts, C, eps, MinPts)      }   }}expandCluster(P, NeighborPts, C, eps, MinPts) {   add P to cluster C   for each point P' in NeighborPts {       if P' is not visited {         mark P' as visited         NeighborPts' = regionQuery(P', eps)         if sizeof(NeighborPts') >= MinPts            NeighborPts = NeighborPts joined with NeighborPts'      }      if P' is not yet member of any cluster         add P' to cluster C   }}regionQuery(P, eps)   return all points within P's eps-neighborhood (including P)

让我们来看看某一个实际的聚类效果图：

密度最大值聚类算法

密度最大值聚类是一种简洁优美的聚类算法，可以识别各种形状的类簇，并且参数很容易确定。这是2014年最新的论文，大家可以点击此处看看论文。 
首先，让我们给出以下定义：

局部密度: 

其中

• 是一个截断距离，即为所有与对象i的距离小于的对象的个数。由于该算法只对的相对值敏感，所以对的选择是稳健的，一种推荐做法就是选择使得平均每一个点的邻居数为所有点的1%-2%。

高局部密度点距离: 

在密度高于对象i的所有对象中，到对象i最近的距离，即为高局部密度点距离。对于密度最大的那个对象，设置。 
- 如果将局部密度看成是个体i的财富值，那么高局部密度点距离可以看做是比个体i有钱并且距离个体i最近的距离。 
- 只有那些密度是局部或者全局最大的点才会有远大于正常值的高局部密度点距离。

簇中心的识别： 
那些有着比较大的局部密度和很大的高局部密度点距离被认为是簇的中心；高局部密度点距离很大但是局部密度较小的点时异常点。 
下面的图可以给我们清晰演示簇中心与异常点的识别。左图为点的分布，右图横坐标为，纵坐标为。 

从右图可以识别出：点1、10都是簇中心，而点28、27、26都是异常点，再反过来看左图，确实如此！

簇中心找到后，剩余的每个点就会被划分到比该点具有更高密度的最近邻居点所属簇当中去，当然，如果该点的最近邻居点还未找到所属簇，那么就会以最近邻居点作为新的起始点，寻找具有更高密度的最近邻居点所属簇。这个过程就相当于：

现在有五个党派（簇中心），国内的每个人都要表明自己是哪个党派的，但是公民A并不清楚他应该属于哪个党派，于是他就去问他的大哥B（比A具有更高密度的最近邻居点），A觉得我大哥B属于哪个党派我就属于哪个党派吧；但是有可能公民B也不知道自己应该属于哪个党派呀，B就回去问他的大哥C，跟随C所属的党派。以此类推下去，所有人都会找到自己的党派（簇中心）。

可靠性分析：对边界和噪声的认识 
在聚类分析中，通常需要确定每个点划分给某个簇的可靠性。在该算法中，可以首先为每个簇定义一个边界区域（border region）,即划分给该簇但是距离其他簇的点的距离小于的点的集合。然后为每个簇找到其边界区域的局部密度最大的点，令其局部密度为. 
该簇中所有局部密度大于的点被认为是簇核心的一部分（即将该点划分给该类簇的可靠性很大），其他的点被认为是该类簇的光晕(halo),即可以认为这些点是噪声。

举例分析：下图中A图为生成数据的概率分布，B，C二图为分别从该分部中生成4000,1000个点。D,E分别是B,C两组数据的决策图(decision tree),可以看到两组数据都只有五个点有较大的较大的\delta_i$。这些点作为簇的中心，在确定了簇中心之后，每个点被划分到各个簇中（彩色点），或者是划分为簇的光晕（黑色点）。F图展示的是随着抽样点数量的增多，聚类的错误率在逐渐下降，说明该算法是鲁棒的。 

谱聚类

谱聚类算法涉及到较多的矩阵分析中的知识，在正式介绍算法之前，容我简单介绍两个知识点： 
1. 实对称矩阵的特征值是实数； 
2. 实对称矩阵的不同对称值对应的特征向量是正交的；

看到这两个都是关于实对称矩阵的知识点，我们头脑中的某些知识点就要翻出来了。那么我们来正式介绍谱和谱聚类吧~

谱：方阵作为线性算子，它的所有特征值统称为方阵的谱。 
谱半径：方阵的谱半径即为最大的特征值；矩阵A(非方阵)的谱半径：（AA）的最大特征值； 
谱聚类：它是一种基于图论的聚类方法，通过对样本数据的拉普拉斯矩阵（具体含义稍后解释）的特征向量进行聚类，从而达到对样本数据聚类的目的。

谱聚类包含概念

相似度图（similarity graph）：给定一组数据,记任意两个点之间的相似度（是距离的减函数）为,那么最初始的聚类目标就是将数据集中的数据依据相似度分组，同一组的数据相似度较高，不同组之间的数据相似度较低。假如我们没有获取到比相似度更多的信息，那么一个比较nice的展示数据的方式就是相似度图G=（V,E）【其中V表示点集，E表示边集，图模型的基本概念】。图模型中点集V=,顶点表示的是样本数据。当点和之间的相似度大于某一特定阈值的时候，那么我们说和是相连着的，和之间的边的权重值为。注：这里的相似度图是无向的，即相似度函数须满足。

此时，聚类问题可以转化为相似度图来表示：划分的原则是：组间的边的权重最小化（意味着组间的数据相似度要尽可能低），组内的边的权重最大化（意味着组内的数据相似度尽可能高）。下图即为一个图模型聚类的示例： 
 
在上图中，一共六个定点，定点之间的连线上的数值表示两个顶点的相似度。现在要将这六个点划分为2个类，怎么分割？根据谱聚类的思想，应该去掉用虚线表示的那条边，进而分为两类。

邻接权重矩阵：如果和之间的相似度大于一定的阈值（这个阈值需要先验给定），那么可以认为这两个点是连接的，满足条件的记为。如果即表示这两个点未发生连接。同样的，由于相似度图G是无向的，那么有。这样，由构成邻接权重矩阵W。

度矩阵：图G中某一个点V的度可以定义为： 

度矩阵D=是一个对角矩阵。

子图：图G的子图相当于全集的子集，即图G中的部分点构成的一个新的图。

• 子图的指示向量：设置子图A的指示向量，当A时候，，否则。
• 如果A和B是图G的不相交子图（无交集），那么定义子图的连接权为： 
  
• 子图A的size测量指标有两个：(1):子图中的顶点个数； (2)。

相似度图G的建立方法

有多种将给定数据集以及对称性相似度（或距离）代入到相似度图的方法。构建相似度图G时，主要是要明确数据点相邻点之间的关系。

1. 近邻图：在这种图中，我们将距离小于的点连接起来。两点间的权值可以取为距离的倒数。这种建立方法的难点在于的确定。

2. k近邻图：在这种图中，以点为例，如果是的k近邻之一，那么将与连起来。当然，这会造成相似度图G是一种有向图，因为是的k近邻之一，但是不一定是的k近邻之一。很像我们生活中：你是我的最好朋友，但是我不一定是你的最好朋友。有两种方法可以将此时的图化为无向图：（1）忽略边的方向性，即只要满足是的k近邻之一或者是的k近邻之一，那么我们就将与连起来。（2）只有满足是的k近邻之一而且是的k近邻之一，我们才将与连起来。依照（1）方法所得的图被称为k近邻图，依照（2）方法所得的图被称为互k近邻图。

3. 全连接图：在这种图中，我们将所有的点都相互连接起来，同时所有的边的权重设置为相似度。这种图建立的关键在于相似度函数的确立，相似度函数能够较好地反映实际中的相邻关系，那么效果较好。一个相似度函数的例子为高斯相似度函数：,其中参数的选取是个难点。该方法中的参数的作用与近邻图中的的作用是相似的。

让我们看看实际上这几种建立方法的效果如何呢？下图即为三种方法的对比。假设原图中的上半部分为"月牙部分"，下半部分为"高斯部分" 

• 近邻图:=0.3,"月牙部分"已经非常紧密连接了，但是"高斯部分"很多都没连接。当数据有不同的"密度"时，往往发生这种问题。
• k近邻图:可以解决数据存在不同密度时有些无法连接的问题，甚至低密度的"高斯部分"和高密度的"月牙部分"也能够连接。同时，虽然两个"月牙部分"的距离比较近，但k近邻还可以把它们区分开。

• 互k近邻图:它趋向于连接相同密度的部分，而不连接不同密度的部分。这种性质介于近邻图和k近邻图之间。如果需要聚类有不同的密度，这种性质非常有效。

• 经验总结：首先尝试使用k近邻图。

拉普拉斯矩阵及其性质

非正则化拉普拉斯矩阵及其性质

非正则化拉普拉斯矩阵定义如下： 

其中矩阵是度矩阵，为权重矩阵。拉普拉斯矩阵满足如下性质：

1. 对于任意向量,满足：
   .
2. 矩阵是对称半正定矩阵。
3. 矩阵的最小特征值为0，与之相对应的特征向量是常数向量.
4. 矩阵有个非负实特征值。

证明性质1： 

证明性质3：当向量为常数向量时，我们有,代入到性质1中有,由半正定矩阵的性质我们可以得到：。由此可以得到矩阵的特征值0对应的特征向量为常数向量。又因为半正定矩阵的特征值非负，可以知道0是矩阵的最小特征值。

性质2、4可以很容易得到。

定理：令相似度图G是权值非负的无向图，拉普拉斯矩阵的特征值为0的重根个数等于图G中连通分量数（子图数）。记G的连通分量为，则特征值0的特征向量即为指示向量。

定理的证明 
首先，让我们先来证明的情况，即特征值为0的个数只有一个。由定理的结论可以知道，此时的图是全连接的。假设此时特征值0所对应的特征向量表示为，那么有： 

由于此时的图是全连接的，那么,那么必然有。即此时的特征向量为常数向量。由于图是全连接的，此时的指示向量也是常数向量。得证。

接下来让我们考虑有k个连通分量的情况。不失一般性的前提下，我们假设顶点都是按照连通分量的归属来排序的。这种情况下的权值矩阵必然是对角块的形式（只有连通分量内的权值才不为0，其余全为0），由定义可以知道拉普拉斯矩阵也是对角块形式 

在连通分量内部，拉普拉斯矩阵是全连接的，特征值为0的个数只有一个（无重根，仅仅在连通分量内部而言），可知此时的特征向量（等于指示向量）为常数向量。 
拉普拉斯矩阵由众多小的拉普拉斯矩阵组成，的特征向量实际上就是方块的特征向量（常数向量）+其他方块的位置补零得到的（由于指示向量等于特征向量，考虑指示向量的含义即可理解）。

正则化拉普拉斯矩阵及其性质

有两个矩阵都被称为正则化拉普拉斯矩阵，两个矩阵的定义是相似的，具体定义如下： 

矩阵被称为对称拉普拉斯矩阵，被称为随机游走(random walk)拉普拉斯矩阵.

注：对角矩阵的(-1/2)次幂依然是对角矩阵，其对角线元素为原矩阵对角线元素的(-1/2)次幂。

正则化拉普拉斯矩阵满足如下性质：

1. 对于任意向量,满足：
   .
2. 是的特征值和特征向量，当且仅当是的特征值和特征向量;
3. 是的特征值和特征向量，当且仅当是的特征值和特征向量;
4. 矩阵和是半正定的，有个非负实特征值。

定理：令相似度图G是权值非负的无向图，拉普拉斯矩阵和的特征值为0的重根个数等于图G中连通分量数（子图数）。记G的连通分量为，则特征值0的特征向量又下列指示向量确定。 

谱聚类算法

在开始介绍算法流程之前，让我们来定义：给定数据集，相似度矩阵是对称矩阵（因为相似度图G是无向的）。

非正则化谱聚类算法流程： 
输入：相似度矩阵,需要聚类的数目.

• 构建相似度图G（具体方法参考前一部分"相似度图G的建立方法"），令为权重矩阵。
• 计算非正则化拉普拉斯矩阵。
• 计算矩阵的按照从小到大顺序排列的前个特征向量。
• 将这k个特征向量（列向量）构造成矩阵。
• 对于，将矩阵的第i行表示成向量。
• 在k维空间里，使用k均值聚类算法将n个数据点聚成k个簇。

输出：聚类结果,其中。

正则化谱聚类算法有两种，取决于采用哪一种正则化拉普拉斯矩阵。实际上，正则化谱聚类算法的流程和非正则化谱聚类算法的流程很相似。

正则化谱聚类算法流程(随机游走拉普拉斯矩阵)： 
输入：相似度矩阵,需要聚类的数目.

• 构建相似度图G（具体方法参考前一部分"相似度图G的建立方法"），令为权重矩阵。
• 计算非正则化拉普拉斯矩阵。
• 计算特征值问题的按照从小到大顺序排列的前个特征向量。
• 将这k个特征向量（列向量）构造成矩阵。
• 对于，将矩阵的第i行表示成向量。
• 在k维空间里，使用k均值聚类算法将n个数据点聚成k个簇。

输出：聚类结果,其中。

---

正则化谱聚类算法流程(对称拉普拉斯矩阵)： 
输入：相似度矩阵,需要聚类的数目.

• 构建相似度图G（具体方法参考前一部分"相似度图G的建立方法"），令为权重矩阵。
• 计算正则化拉普拉斯矩阵。
• 计算矩阵的按照从小到大顺序排列的前个特征向量。
• 将这k个特征向量（列向量）构造成矩阵。
• 将矩阵的每一行都归一化到1，得到矩阵,其中。
• 对于，将矩阵的第i行表示成向量。
• 在k维空间里，使用k均值聚类算法将n个数据点聚成k个簇。

输出：聚类结果,其中。

经验：首选游走拉普拉斯矩阵对应的谱聚类算法。

通过对比三种算法的具体过程，我们发现它们是非常相似的，只不过拉普拉斯矩阵不同从而导致了一丝不同而已。谱聚类算法的主要trick在于将原始数据转化为的过程（很多时候都是升维的过程）。这个过程之所以有效主要在于拉普拉斯矩阵的性质。经过变化后的数据，往往使用简单的k均值算法就可以实现很好的聚类效果。

那么为什么呢？算法的原理是什么呢？

从图的分割的角度来看待谱聚类

聚类的一开始，我们想的是根据数据点的相似性的大小来分成不同的组。进一步，将聚类问题转化为相似度图的划分来表示：划分的原则是：组间的边的权重最小化（意味着组间的数据相似度要尽可能低），组内的边的权重最大化（意味着组内的数据相似度尽可能高）。在这个维度来看，让我们来分析谱聚类问题如何近似为图的切割问题。

给定一个带权重矩阵的相似度图,最简单最直接对一个图划分的方法就是最小化切割问题，定义最小化切割问题：选择一组划分：，最小化下面的式子： 

其中为的补；表示子图的连接权（之前已有定义）。

在实际操作中，发现上述的目标函数存在问题：在很多情况下，mincut的解，将图分成了1个点和其余个点。这1个点很有可能是个异常点，这样的切割对于实际聚类来讲是没有什么意义的。为了避免这个问题，目标函数应该显式要求所得到的划分应该足够大。最常见的修正后的目标函数有两个：RatioCut 以及 Ncut,定义如下： 

其中的表示的顶点个数，表示的所有边的权值和，在之前章节都有定义。 
那么我们就来分析一下为什么RatioCut和Ncut就能够解决问题呢？上述目标函数以的点数或者权值和作为被除数，使得函数在所有相等的时候达到最小值，也就是所有顶点被平均分配到每个子图时最小；函数在所有相等的时候达到最小值。这样，RatioCut和Ncut在一定程度上能够使得最终得到的"簇"趋向于平衡。

上面只是讲了目标函数的修正，那么具体如何求解这个最值问题呢？我们重点分析RatioCut的解法，实际上RatioCut最终推导得到的最优解即为非正则化谱聚类算法的计算结果。而有兴趣的朋友可以去研究一下Ncut最终推导得到的最优解，实际上即为正则化谱聚类算法的计算结果。实际上RatioCut与Ncut的推导过程是相似的。

k=2时的RatioCut

当，这个事实目标函数就面成了求解如下问题： 

给定一个子图，其中表示全图，我们定义向量，其中 

向量表达式定下来以后，由下面的等式，我们可以将求解RatioCut的问题转化为求解拉普拉斯矩阵的问题： 

而经过研究发现，向量还满足一定的约束条件： 

上式表明向量和常数向量是正交的。另外 

将所有加在一起，我们可以将之前的RatioCut问题转化为： 

定义如上，且

该目标函数的自变量部分，是不同的子图划分A的选择，从而得到不同的向量。每进行一次子图划分，就会得到一个向量，进而可以得到不同的值。优化的目标，是使得该目标函数取值最小。但是由于只能够取2个值（离散化的定义域），而n个数据点的子图划分选择可以有指数级的选择，如果没有一个有效措施的话，这个问题的求解是一个NP问题。

为了解决这个困难，我们将目标条件做一个放松relaxation，将的严格定义用的性质来代替，向量各个分量的取值从离散化的若干个值拓展到整个实数集，那么目标函数变为 

如此放松之后，问题就变成了一个非常标准的线性代数问题，找到这样一个，它和常数向量是相交的，同时它的模等于，最终使得函数最小即可。

幸运的是，根据Rayleigh-Ritz定理，该目标函数的解即为矩阵的次小特征向量。

分为k个子图时的RatioCut

与划分为2个子图的做法相似，我们可以得到与是的RatiCut相似的推导过程与结论： 
需要将一个图划分为k个子图，我们定义k个指示向量如下： 

将这k个特征向量（列向量）构造成矩阵。发现矩阵中的每一列向量都是相互正交的，即有。

和之前的计算过程相似，我们可以得到： 

更进一步，有： 

于是可得： 

经过放松relaxation之后，原问题变成： 

根据Rayleigh-Ritz定理，该目标函数的解矩阵即为矩阵的从小到大排序的前个特征向量。实际上，细心的朋友会发现，此时的矩阵即为非正则化谱聚类算法流程中的矩阵。实际上，从RatioCut推导出问题最优解的过程正是在解释“为什么谱聚类的算法是可行的过程”。

后记

除了从图的切割的角度来去看待谱聚类算法以外，还可以用随机游走、扰动论等理论来解释。

前路漫漫！

---

参考文献： 
Clustering by fast search and find of densitypeak. Alex Rodriguez, Alessandro Laio 
A tutorial on spectral clustering，Ulrike von Luxburg, 2007

谱聚类算法详解