Python实现PCA降维算法及其应用

问题导读：
1. PCA 降维算法原理是什么？
2. 如何使用Python实现PCA降维算法？
3. PCA降维后的数据常用于什么场景？
4. Kmeans算法原理是什么？
5. 如何使用Spark Mllib自带的Kmeans算法?



前言

首先我们先确定一个事实，就是我们在做ML(机器学习)的时候，绝不是算法第一的。

我们在很多时候选择一个或者说决定一个模型开始训练，我们首先应该想的是：

• 数据来源(ETL的T)
• 数据的格式化(ETL的T)
• 数据采用的训练模型
• 模型的展示
  

PCA降维算法

PCA简单的说，它是一种通用的降维工具。在我们处理高维数据的时候，
可以（原文“了能”）降低后续计算的复杂度，在“预处理”阶段通常要先对原始数据进行降维，
而PCA就是干这个事的 本质上讲，PCA就是将高维的数据通过线性变换投影到低维空间上去
具体的数学原理我推荐这个（http://blog.csdn.net/xiaojidan2011/article/details/11595869）

降维当然意味着信息的丢失，不过鉴于实际数据本身常常存在的相关性，我们可以想办法在降维的同时将信息的损失尽量降低。

总结一下PCA的算法步骤：
设有m条n维数据。
1）将原始数据按列组成n行m列矩阵X
2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）求出协方差矩阵C=
4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P
6）Y=PX即为降维到k维后的数据

我们用python试着实现下：
[mw_shl_code=python,true]# -*- coding: utf-8 -*-

from math import *
import random as rd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def zeroMean(dataMat):
    meanVal = np.mean(dataMat,axis = 0)#计算该轴上的统计值（0为列，1为行）
    newData = dataMat - meanVal
    return newData,meanVal

def pca(dataMat,percent=0.99):
    '''求协方差矩阵
    若rowvar=0，说明传入的数据一行代表一个样本，若非0
    说明传入的数据一列代表一个样本。因为newData每一行代表一个样本，所以将rowvar设置为0 '''
    newData,meanVal=zeroMean(dataMat)
    covMat=np.cov(newData,rowvar=0)
    eigVals,eigVects = np.linalg.eig(np.mat(covMat))
    n=percentage2n(eigVals,percent)          #要达到percent的方差百分比，需要前n个特征向量
    print str(n) + u"vectors"
    eigValIndice=np.argsort(eigVals)            #对特征值从小到大排序  
    n_eigValIndice=eigValIndice[-1:-(n+1):-1]   #最大的n个特征值的下标  
    n_eigVect=eigVects[:,n_eigValIndice]        #最大的n个特征值对应的特征向量  
    lowDDataMat=newData * n_eigVect               #低维特征空间的数据  
    reconMat=(lowDDataMat * n_eigVect.T) + meanVal  #重构数据  
    return reconMat,lowDDataMat,n

def percentage2n(eigVals,percentage):  
    sortArray=np.sort(eigVals)   #升序  
    sortArray=sortArray[-1::-1]  #逆转，即降序  
    arraySum=sum(sortArray)  
    tmpSum=0  
    num=0  
    for i in sortArray:  
        tmpSum += i  
        num += 1
        if tmpSum >= arraySum * percentage:  
            return num

if __name__ == '__main__':
    data = np.random.randint(1,10,size = (3,5))
    fig = plt.figure()
    ax = plt.subplot(111,projection='3d')
    #ax.scatter(data[0],data[1],data[2],c='y') #绘制数据点
    ax.set_zlabel('Z') #坐标轴
    ax.set_ylabel('Y')
    ax.set_xlabel('X')
    #plt.show()
    print data
    fin = pca(data,0.9)
    mat =fin[1]
    print mat
    #ax.scatter(mat[0],mat[1],mat[2],c='y') #绘制数据点
    #plt.show()[/mw_shl_code]
应用

我们对PCA后降维后的数据最直接的应用是聚类，这里我们还是选择kmeans算法：

kmeans

具体算法可百度，我这里提下不同类型变量相异度计算方法：

• 标量也就是无方向意义的数字，也叫标度变量:
  
  • 一种很自然的想法是用两者的欧几里得距离来作为相异度，欧几里得距离的定义如下：
    
    
• 对于向量，由于它不仅有大小而且有方向，所以闵可夫斯基距离不是度量其相异度的好办法，一种流行的做法是用两个向量的余弦度量，其度量公式为：
  
  

k均值算法的计算过程非常直观：  
1、从D中随机取k个元素，作为k个簇的各自的中心。  
2、分别计算剩下的元素到k个簇中心的相异度，将这些元素分别划归到相异度最低的簇。  
3、根据聚类结果，重新计算k个簇各自的中心，计算方法是取簇中所有元素各自维度的算术平均数。
4、将D中全部元素按照新的中心重新聚类。
5、重复第4步，直到聚类结果不再变化。
6、将结果输出。

我们上代码，是基于spark MLlib的：

[mw_shl_code=scala,true]import org.apache.spark.SparkConf
import org.apache.spark.SparkContext
import org.apache.spark.mllib.clustering.KMeans

import org.apache.spark.mllib.linalg.Matrix
import org.apache.spark.mllib.linalg.Vectors
import org.apache.spark.mllib.linalg.distributed.RowMatrix

import scala.util.Random

object PCA {
    def getRandom(num: Int) = {
        (for(i <- 1 to num ) yield "%1.1f".format(new Random().nextDouble()*10).toDouble ).toArray
    }
    def main(args: Array[String]): Unit = {
        val line = 5
        val conf =new SparkConf().setAppName("PCA");
        val sc = new SparkContext(conf)
        val data =(for(i <- 1 to line ) yield Vectors.dense(getRandom(10))).toArray
        //我们产生一个5X10的矩阵
        val dataRDD = sc.parallelize(data, 2)
        //RowMatrix 分布式矩阵,RowMatrix.numRows，RowMatrix.numCols
        val mat: RowMatrix = new RowMatrix(dataRDD)
        val pc: Matrix = mat.computePrincipalComponents(3)//PCA只需要保留前3个特征
        //得到的矩阵结果
        val projected: RowMatrix = mat.multiply(pc)
        val newdateRDD = projected.rows
        // 打印出降维的数据
        val numIterations = 20//迭代的次数
        val ks:Array[Int] = Array(2,3,4,5)
           ks.foreach(cluster => {
            val model = KMeans.train(newdateRDD, cluster,numIterations,1)
            val ssd = model.computeCost(newdateRDD)
            println(" when k=" + cluster + " -> "+ ssd)
        })
        val Knum = 3  //将目标数据分成几类
        //将参数，和训练数据传入，形成模型
        val clusters1 = KMeans.train(dataRDD, Knum , numIterations)//训练原始数据
        val clusters2 = KMeans.train(newdateRDD, Knum , numIterations)//训练降维数据


        val collect = projected.rows.collect()
        println("主成分投影列矩阵:")
        collect.foreach { vector => println(vector) }
        val source = mat.rows.collect()
        println("原始投影列矩阵:")
        source.foreach { vector => println(vector) }
        /*println("原始数据中心分布:")
        for (c <- clusters1.clusterCenters) {
              println("  " + c.toString)
        }
        println("训练降维数据中心分布:")
        for (c <- clusters2.clusterCenters) {
              println("  " + c.toString)
        }*/
        for( i <- 0 to source.length - 1) {
            println("训练原始数据:" + source(i).toString + "属于" + clusters1.predict(source(i)).toString + "类")
        }
        for( i <- 0 to collect.length -1) {
            println("训练降维数据:" + collect(i).toString + "属于" + clusters2.predict(collect(i)).toString + "类")
        }        
      }
    }[/mw_shl_code]

后记

大概就这样，最近遇到一个问题，就是数据的维度太多，特征向量都不知道选什么。很是头疼~


作者：michaelgbw
來源：http://www.jianshu.com/p/16f835a94444

感谢分享