聚类算法，不是分类算法。
分类算法是给一个数据，然后判断这个数据属于已分好的类中的具体哪一类。
聚类算法是给一大堆原始数据，然后通过算法将其中具有相似特征的数据聚为一类。
这里的k-means聚类，是事先给出原始数据所含的类数，然后将含有相似特征的数据聚为一个类中。
聚类属于无监督学习，以往的回归、朴素贝叶斯、SVM等都是有类别标签y的，也就是说样例中已经给出了样例的分类。而聚类的样本中却没有给定y，只有特征x，比如假设宇宙中的星星可以表示成三维空间中的点集。聚类的目的是找到每个样本x潜在的类别y，并将同类别y的样本x放在一起。比如上面的星星，聚类后结果是一个个星团，星团里面的点相互距离比较近，星团间的星星距离就比较远了。
K-means算法是将样本聚类成k个簇（cluster），具体算法描述如下：
K是我们事先给定的聚类数，代表样例i与k个类中距离最近的那个类，的值是1到k中的一个。质心代表我们对属于同一个类的样本中心点的猜测，拿星团模型来解释就是要将所有的星星聚成k个星团，首先随机选取k个宇宙中的点（或者k个星星）作为k个星团的质心，然后第一步对于每一个星星计算其到k个质心中每一个的距离，然后选取距离最近的那个星团作为，这样经过第一步每一个星星都有了所属的星团；第二步对于每一个星团，重新计算它的质心（对里面所有的星星坐标求平均）。重复迭代第一步和第二步直到质心不变或者变化很小。1{Ci=j} ;这个被称作为示性函数，即使如果Ci=j；则最后1{Ci=j}的值为1，否则为0；上式中分母为统计第j个类的个数，分子为求每个样本对应维度“坐标”的和，相除后就是新的聚类中心了，其实就是每个样本的维度“坐标”求平均值（均值=总和/个数）
下图展示了对n个样本点进行K-means聚类的效果，这里k取2。
K-means面对的第一个问题是如何保证收敛，前面的算法中强调结束条件就是收敛，可以证明的是K-means完全可以保证收敛性。下面我们定性的描述一下收敛性，我们定义畸变函数（distortion function）如下：
J函数表示每个样本点到其质心的距离平方和。K-means是要将J调整到最小。假设当前J没有达到最小值，那么首先可以固定每个类的质心，调整每个样例的所属的类别来让J函数减少，同样，固定，调整每个类的质心也可以使J减小。这两个过程就是内循环中使J单调递减的过程。当J递减到最小时，和c也同时收敛。（在理论上，可以有多组不同的和c值能够使得J取得最小值，但这种现象实际上很少见）。
由于畸变函数J是非凸函数，意味着我们不能保证取得的最小值是全局最小值，也就是说k-means对质心初始位置的选取比较感冒，但一般情况下k-means达到的局部最优已经满足需求。但如果你怕陷入局部最优，那么可以选取不同的初始值跑多遍k-means，然后取其中最小的J对应的和c输出。
在kmeans初始化k个类别时，由于初始化具有随机性，如果选取的初始值点不同可能导致最后聚类的效果跟想象中的效果相差很远，这也就是kmeans的局部收敛问题。解决这个问题一般采用的方法是进行多次kmeans，然后计算每次kmeans的损失函数值，取损失函数最小对应的那个结果作为最终结果。
在kmeans中比较棘手的另一个问题是类别k的选择。因为有的数据集用不同的k来聚类都感觉比较合适，那么到底该用哪个k值呢？通常情况下的方法都是采用”elbow”的方法，即做一个图表，该图的横坐标为选取的类别个数k，纵坐标为kmeans的损失函数，通过观察该图找到曲线的转折点，一般这个图长得像人的手，而那个像人手肘对应的转折点就是我们最终要的类别数k，但这种方法也不一定合适，因为k的选择可以由人物确定，比如说我就是想把数据集分为10份（这种情况很常见，比如说对患者年龄进行分类），那么就让k等于10。
由于是初学k-means所以在看懂原理后，在百度上面搜索算法代码，结果搜索类好几种matlab 的k-means算法，都是各有千秋；有些代码真是可读性，变量命名，让人捉急看着头痛。最后还是看andrew ng发布的代码，感觉还是正规军牛啊，代码虽然理解起来有些难，因为里面用到了很多矩阵造作运算，但是代码的结构很清晰，效率很优化。而且还考虑了异常情况，即使当样本数较少，数据分类数较多，而造成出现空类型的时候把对应的聚类坐标变为0。
代码结构：
一：初始化聚类中心
方法1：常用的中心初始化方法就是随机地从样本中挑k个出来作为k个初始的聚类中心。但这不是个明智的选择。它有可能会导致图像趋于稠密聚集某些区域，因为如果随机选择的样本，也就是训练样本本身就在某个区域分布非常密，那么我们随机去选择聚类中心的时候，就会出现就在这个数据分布密集的地方被选出了很多的聚类中心（因为原数据分布密集，所以被选中的概率会大）。这是不好的，因为本身这个密集的区域，一个聚类中心才是好的，你硬搞出几个聚类中心，好好的队伍就会被强拆了，还搞得其他比较分散的数据点找不到归宿，因为聚类中心离它好遥远，大部分的聚类中心都让富二代给占有了，社会资源分配不当啊，哈哈。这样还会出现一种情况，就是如果样本分布密集的话，它们都几乎从属于一个聚类中心，那么其他聚类中心就几乎没有样本属于它。这样也不好。
方法2：一个比较好的初始化方式就是从一个正态分布中随机初始化聚类中心，然后归一化他们到单位长度。
Matlab中函数randn（m，n），用来产生均值为0，方差为1的符合正太分布放入m行n列的矩阵。此处需要把产生的数据通过变换到样本数据的维度。
Centrids=b+ randn(m，n)*a； a和b为系数，变换数据的范围
二：计算样本距离
方法1：单个循环处理
使用for循环，每次处理一个样本Xi，然后计算Xi到每个聚类中心的距离，找到最小的距离，最小距离对应的中心即为次样本的中心。
此代码清晰明确，每次处理一个样本，然后依次计算此样本到各个聚类中心点的距离。但是此代码处理速度慢，没有发挥matlab矩阵运算的优势。
方法2：矩阵批处理
采用批处理的方式，每次处理BATCH_SIZE（如果样本个数较少，可以让BATCH_SIZE=样本个数）个样本，通过矩阵运算，计算BATCH_SIZE个样本中每个样本的分类。此种方法计算快，效率高。
优化分析：
因为min(a-b)^2等价于求min(a^2+b^2-2*a*b)等价于求max(a*b-0.5*a^2-0.5*b^2),
a为x2 = sum(X.^2,2);每一个样本元素的平方和，x2指每个样本点与原点之间的欧式距离。
b为c2 = 0.5*sum(centroids.^2,2);c2表示类别中心点到原点之间的欧式距离
则每次从a中取出一个数据与b中所有中心点作比较时，此时a中的该数据可以忽略不计，只跟b有关。即原式等价于求：
[val,labels]=max(a*b-0.5*a^2)=max(bsxfun(@minus,centroids*X(i:lastIndex,:)',c2))
val为BATCH_SIZE大小的行向量，labels为每个样本经过一次迭代后所属的类别标号
Bxsfun函数为matlab自带的矩阵运算函数，可以快速实现举证的加减乘除运算。
bsxfun(@minus,centroids*X(i:lastIndex,:)',c2)=centroids*X(i:lastIndex,:) - c2
若c2和centroids*X(i:lastIndex,:)的矩阵维度不符，自动复制c2的行或者列，是两个矩阵的维数相等。
计算示性矩阵S。
矩阵S的维数为BATCH_SIZE*k（k为分类数）；每行为BATCH_SIZE中每个样本的分类情况，每行有k个元素，其中只有一个数为1，此样本属于此列对应的类别，其他都为0。由于矩阵S中大部分数据都为0，所以利用矩阵的稀疏表示可以加快矩阵运算，减少存储空间。
S = sparse(1:m,labels,1,m,k,m); %矩阵S的系数表示
summation = summation + S'*X(i:lastIndex,:);其结果为每个类别对应维度的坐标求和，
counts = counts+ sum(S,1)';统计每个类别的样本数。
三：更新聚类中心
centroids = bsxfun(@rdivide, summation, counts);计算新的聚类中心
run_kmeans.m:
输入：x为输入数据，k为类别数，iterations为迭代循环次数
输出：样本的聚类中心
function centroids = runkmeans(X, k, iterations)
x2= sum(X.^2,2);
centroids =randn(k,size(X,2))*0.1;%X(randsample(size(X,1), k), :);
BATCH_SIZE=1000;
for itr =1:iterations%iterations为kemans聚类迭代的次数
fprintf('K-means iteration %d / %d\n', itr, iterations);
c2 = 0.5*sum(centroids.^2,2);
summation = zeros(k, size(X,2));
counts = zeros(k, 1);
loss =0;
for i=1:BATCH_SIZE:size(X,1)
lastIndex=min(i+BATCH_SIZE-1, size(X,1));
m = lastIndex - i + 1;%m=1000
[val,labels] =max(bsxfun(@minus,centroids*X(i:lastIndex,:)',c2));
loss = loss + sum(0.5*x2(i:lastIndex) - val');
S = sparse(1:m,labels,1,m,k,m); % labels as indicator matrix
summation = summation + S'*X(i:lastIndex,:);
counts = counts + sum(S,1)';
end
centroids = bsxfun(@rdivide, summation, counts);
% just zap empty centroids so they don't introduce NaNseverywhere.
badIndex = find(counts == 0);
centroids(badIndex, :) = 0;%防止出现无穷大的情况
end
参考文献：
Deep Learning论文笔记之（三）单层非监督学习网络分析
Deep learning：二十五(Kmeans单层网络识别性能)