谱聚类的实现
Posted zbxzc
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了谱聚类的实现相关的知识，希望对你有一定的参考价值。
开始用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 3
vector<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 Sun
import heapq
from functools import partial
from numpy import *
from scipy.linalg import *
from scipy.cluster.vq import *
import pylab
def 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 vecs
def 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 W
def 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 - W
evals, 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最近的距离，即为高局部密度点距离。对于密度最大的那个对象，设置 Spectral clustering（谱聚类）算法的实现