PCA/SVD用于人脸图像数据

148 阅读7分钟

spark mllib 提供两种相似的降低维度的模型:主成分分析(PCA)和奇异值分解(SVD)用于在无监督学习中降低数据维度。
应用场景:
探索性数据分析;
提取特征区训练其它机器学习模型;
降低大型模型在预测阶段的存储和计算需求
把大量文档缩减为一组隐含话题;
当数据维度很高时,使得学习和推广更容易(如,处理文本、声音、图像、视频等非常高维的数据时)

/**
  * Created by raini on 16-4-28. spark-1.6
  */
    val path = "file:/home/raini/data/lfw/*"
    val rdd = sc.wholeTextFiles(path)
    val files = rdd.map { case (fileName, content) =>
        fileName.replace("file:", "")
    }
    println(files.first) ///home/raini/data/lfw/Ali_Naimi/Ali_Naimi_0002.jpg
    println(files.count) //1054
/**
渲染器  文件类型     描述

AGG png 光栅图 –使用 Anti-Grain Geometry 高质量渲染引擎
PS ps eps 矢量图 – Postscript 输出
PDF pdf 矢量图– 可携带格式
SVG svg 矢量图 – 可伸缩矢量图形
Cairo png ps pdf svg … 矢量图 – Cairo图
GDK png jpg tiff … 光栅图 – gimp

    import numpy as np
    from matplotlib import *
    from matplotlib.pyplot import *
    from pylab import *
    path = "/home/raini/data/lfw/Ali_Naimi/Ali_Naimi_0002.jpg"
    ae = imread(path)
    matplotlib.use('TkAgg')
    imshow(ae)
   pylab.savefig('/home/raini/data/temp_Ali_Naimi_0002.png')
    **/
/** 提取面部图片作为向   
 //1.载入图片

import java.awt.image.BufferedImage
def loadImageFromFile(path: String): BufferedImage = {
import javax.imageio.ImageIO
import java.io.File
ImageIO.read(new File(path))
}
val aePath = “/home/raini/data/lfw/Ali_Naimi/Ali_Naimi_0002.jpg”

/** aeImage: java.awt.image.BufferedImage = BufferedImage@33ad8f4e: type = 5 ColorModel:
* #pixelBits = 24 numComponents = 3 color space = java.awt.color.ICC_ColorSpace@479b6463
* transparency = 1 has alpha = false isAlphaPre = false
* ByteInterleavedRaster: width = 250 height = 250 #numDataElements 3 dataOff[0] = 2
*/
//2.转换灰度图片并改变图片尺寸
def processImage(image: BufferedImage, width: Int, height: Int): BufferedImage = {
val bwImage = new BufferedImage(width, height, BufferedImage.TYPE_BYTE_GRAY)
val g = bwImage.getGraphics //从原始图片绘制出灰度图
g.drawImage(image, 0, 0, width, height, null) //颜色转换 尺度变化
g.dispose()
bwImage
}

    //转换测试
    val grayImage = processImage(aeImage, 100, 100) //颜色组件从 3-> 1

    /** grayImage: java.awt.image.BufferedImage = BufferedImage@c0717f8: type = 10 ColorModel:
      *  #pixelBits = 8 numComponents = 1 color space = java.awt.color.ICC_ColorSpace@25b6a0f3
      *  transparency = 1 has alpha = false isAlphaPre = false
      *  ByteInterleavedRaster: width = 100 height = 100 #numDataElements 1 dataOff[0] = 0
      */

这里写图片描述这里写图片描述
//3.存储处理过的图片
import javax.imageio.ImageIO
import java.io.File
ImageIO.write(grayImage, “jpg”, new File(“/home/raini/data/aiyaya.jpg”))

/**
from matplotlib.pyplot import *
path = "/home/raini/data/bigege.jpg"
ae = imread(path)
imshow(ae, cmap=plt.cm.gray)
pylab.savefig('/home/raini/data/temp_Naimi_0003.png')
  **/

    //4.提取真实的特征向量 作为降维模型的输入,(BufferedImage-Pixels)将通过打平二维的像素矩阵来构造一维的向量
def getPiexelsFromImage(image: BufferedImage): Array[Double] = {
  val width = image.getWidth()
  val heigh = image.getHeight()
  val pixels = Array.ofDim[Double](width * heigh)
  image.getData.getPixels(0, 0 ,width, heigh, pixels) //pixels-传入一个一维数组
}

    /**
      * 6.为了更直观,组合上面函数,将接受一个图片文件位置 和 需要处理的 宽高,返回一个包含像素的Array[DOuble]值*/
def extractPixels(path: String, width: Int, height: Int): Array[Double] = {
  val raw = loadImageFromFile(path)
  val processed = processImage(raw, width, height)
  getPiexelsFromImage(processed)
}
    /**
      * 7.将函数应用到图片路径RDD->产生包含每张图片的橡树数据RDD*/
val pixels = files.map(f => extractPixels(f, 50, 50))
println(pixels.take(10).map(_.take(15).mkString("", ",", ", ...")).mkString("\n"))//start sep end

/**
240.0,240.0,240.0,240.0,242.0,242.0,242.0,242.0,242.0,242.0, …
123.0,103.0,132.0,153.0,165.0,183.0,175.0,163.0,169.0,167.0, …
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, …
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, …
139.0,141.0,132.0,127.0,145.0,89.0,100.0,114.0,93.0,73.0, …
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, …
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, …
0.0,0.0,0.0,1.0,1.0,0.0,0.0,0.0,1.0,1.0, …
161.0,143.0,115.0,111.0,135.0,151.0,161.0,149.0,131.0,111.0, …
151.0,151.0,158.0,156.0,152.0,150.0,150.0,154.0,155.0,142.0, … */

    /**
      * 8.为每一张图片创建Mllib向量对象*/
import org.apache.spark.mllib.linalg.Vectors
val vectors = pixels.map(p => Vectors.dense(p))
vectors.setName("images-Vectors")
vectors.cache()

    /**
      * 9.在运行降维模型尤其是PCA之前,通常对输入数据---正则化*/
import org.apache.spark.mllib.feature.StandardScaler
val scaler = new StandardScaler(withMean = true, withStd = false).fit(vectors)
    /**注意:
      * 对于稠密的输入数据可以提取平均值,但是对于稀疏的数据,提取平均值将会变得稠密*/
val scaledVectors = vectors.map(v => scaler.transform(v))


    /**
      * 1.训练降维模型--------------------------------------------*/
    //模型需要向量作为输入,但不像聚类直接处理RDD[Vector],
    // PCA和SVD的计算是通过提供基于RowMatrix的方法实现(区别主要是语法不同,RowMatrix仅仅是一个RDD[Vector]的简单法封装)
import org.apache.spark.mllib.linalg.Matrix
import org.apache.spark.mllib.linalg.distributed.RowMatrix
val matrix = new RowMatrix(scaledVectors)
val k = 10
val pcK = matrix.computePrincipalComponents(k) //计算前k个主成分

// val ccss = matrix.computeColumnSummaryStatistics()//统计每一列的信息
// ccss.variance

val rows = pcK.numRows
val cols = pcK.numCols
val acts = pcK.numActives

//rows: Int = 2500 每张图片维度是50×50; 所以我们得到了前10个主成分向量,
//cols: Int = 10 每一个向量的维度都和输入图片的维度一样,可以认为这些主成分是一组包含了原始数据主要变化的隐藏特征
//acts: Int = 25000

    /**在面部识别和图像处理时,这些主成分总是被称为(特征脸),这是因为PCA和原始数据的协防差矩阵的特征值分解非常相关
      * 因为每一个主成分都和原始图像有相同维度,每一个成分本身可以看作是一张图像,这使得我们下面要做的可视化特征脸成为可能*/

    /**
      * 2.可视化特征脸*/
    //使用Breeze、numpy、matplotlib来可视化特征脸

import breeze.linalg.DenseMatrix
val pcaBreeze = new DenseMatrix(rows, cols, pcK.toArray)
//pcaBreeze.data.take(5)

    /**breeze.linalg包中提供了实用的函数把矩阵写到CSV文件中-(保存主成分)*/

import breeze.linalg.csvwrite
csvwrite(new File(“/home/raini/data/pca.csv”), pcaBreeze)

    //在IPython中加载矩阵并以图像的形式可视化主成分

/**
* pca = np.loadtxt(“/home/raini/data/pca.csv”, delimiter=”,”)
* print(pca.shape)
* //(2500, 10)
*
* —————–//绘制这10个特征脸————————
def plot_gallery(images, h, w, n_row=2, n_col=5):
“”“Helper fenction to plot ~ “””
plt.figure(figsize=(1.8 * n_col, 2.4 * n_row))
plt.subplots_adjust(bottom=0, left=.01, right=.99, top=.90, hspace=.35)
for i in range(n_row * n_col):
plt.subplot(n_row, n_col, i + 1)
plt.imshow(images[:, i].reshape((h, w)), cmap=plt.cm.gray)
plt.title(“Eigenface %d” % (i + 1), size=12)
plt.xticks(())
plt.yticks(())

plot_gallery(pca, 50, 50)
pylab.savefig('/home/raini/data/temp_2.png')    * */

这里写图片描述

   //解释特征脸
    //---------可以看到PCA模型有效的提取了反复出现的变化模式,表现了脸部图像的各种特征。
    //有些图好像选择了方向性特征(图69),有些集中表现了发型(457 10),其它似乎和面部特征更相关(眼鼻嘴)

    /**
      * 3.使用降维模型
      *   把-原始数据-投影到用主成分表示的-新的-低维(10维)空间上
      *   降维方法最终的目标是要得到数据更加压缩化的表示,并能包含原始数据之中重要的特征和变化*/
    //用矩阵乘法把 图像矩阵 和 主成分矩阵 相承 实现投影(图像矩阵是分布式的RowMatrix, Spark实现了分布式计算的multiply函数)
val projected = matrix.multiply(pcK)
println(projected.numRows(), projected.numCols())//.computePrincipalComponents(计算前k个主成分)
//(1054,10) 每幅2500维度的图像已经被转换成为一个大小为10的向量

println(projected.rows.take(5).mkString("\n"))
    //这些以向量形式表示的投影后的数据可以用来作为另一个机器学习模型的输入。例如和一些没有脸的图像产生投影数据,共同训练一个面部识别模型*/

    /**
      * 4.PCA 和 SVD模型 的关系
      *   可以使用SVD恢复出相同的主成分向量,并且应用相同的投影矩阵投射到主成分空间
      *    SVD计算产生的(右奇异值向量)等同于 之前计算得到的PCA*/
val svd = matrix.computeSVD(10, computeU = true)
println(s"矩阵(U) dimension: (${svd.U.numRows()}, ${svd.U.numCols()})")
println(s"奇异值向量(s) dimension: (${svd.s.size}, )") //${svd.s.numActives}
println(s"右奇异值向量(V)dimension: (${svd.V.numRows}, ${svd.V.numCols})")
//矩阵(U) dimension: (1054, 10)
//奇异值向量(s) dimension: (10, )//10
//右奇异值向量(V)dimension: (2500, 10) <----矩阵V与PCA完全一样(不考虑正负浮点),验证一下:

def approxEqual(array1: Array[Double], array2: Array[Double], tolerance: Double = 1e-6) = {
  val bools = array1.zip(array2).map{ case (v1, v2) =>
    if (math.abs(math.abs(v1) - math.abs(v2)) > 1e-6) false else true }
  bools.fold(true)(_ & _)
}
println(approxEqual(svd.V.toArray, pcK.toArray))
//true

    /**另一个相关性体现:
      *矩阵U和 对角矩阵S 的乘积-与-PCA中用来把原始图像数据投影到10个主成分构成的空间中的投影矩阵 相等*/
val breezeS = breeze.linalg.DenseVector(svd.s.toArray)

val projected_UmS = svd.U.rows.map{ v =>
  val breezeU = breeze.linalg.DenseVector(v.toArray)
  val multV = breezeU :* breezeS //:* 表示对向量执行对应元素和元素的乘法
  Vectors.dense(multV.data)
}
projected.rows.zip(projected_UmS).map { case (v1, v2) =>
    approxEqual(v1.toArray, v2.toArray)
}.filter(b => true).count()
//Long = 1054---------------确定投影后的每一行和SVD投影后的每一行相等

    /**
      * 5.评价降维模型-估计SVD的k值
      *   PCA/SVD都是确定性模型,这两个模型都确定可以返回多个主成分或者奇异值,因此控制模型唯一参数就是k。
      *   如聚类 增加k总是可以提高模型的表现--不确定性表现在相对误差函数值;pca/svd不确定性表现在k个成分上*/
    //当用PCA/SVD估算k值时,以一个较大的k的变化范围绘制一个奇异值图是很重要的。

val svd300 = matrix.computeSVD(300, computeU = false)//–k=300
val sMatrix = new DenseMatrix(1, 300, svd300.s.toArray)
csvwrite(new File(“/home/raini/data/s300.csv”), sMatrix)

/**
* 用python绘图
import scipy
import numpy as np
import matplotlib.pyplot as plt
s = np.loadtxt(“/home/raini/data/s300.csv”, delimiter=”,”)
print(s.shape)
plt.plot(s)
plt.savefig(‘/home/raini/data/s300_1.png’)

plt.plot(cumsum(s))
plt.yscale('log')
plt.savefig('/home/raini/data/s300_2.png') */
```
    vectors.unpersist()

这里写图片描述
这里写图片描述
//可以看到k在100左右处于拐点,表明k=100时的PCA/SVD足够解释原始数据的变化
//使用降维模型来提高另一个模型的性能,我们可以使用和那个模型相同的评价模型来帮助我们选择k值
//如可以使用AUX矩阵和交叉验证,来为分类模型选择模型参数和为降维模型选择k值,但这会耗费更高的计算资源,因为我们必须重算整个模型的训练的测试过程