图像处理——主成分分析

图像处理——主成分分析
1 引⾔
1.1 维度灾难
分类为例：如最近邻分类⽅法（基本思想：以最近的格⼦投票分类）
问题：当数据维度增⼤，分类空间爆炸增长。

如图1所⽰，
图1 维度增加⽰意图
1.2 解决⽅法
缓解维度遭难的⼀个重⽤途径是降维。

降维是通过某种数学变换，将原始⾼维属性空间转换为⼀个低维“⼦空间”，在这个⼦空间中样本
密度⼤幅度提⾼，距离计算也变得更为容易。

1.3 降维的可⾏性
数据样本虽然是⾼维的，但与我们关⼼的也许仅是某个低维分布，因⽽可以对数据进⾏有效的降维。

1.4 主成分分析
1.4.1 概念
主成分分析(Principal Component Analysys,简称PCV)，是将具有相关性的⾼维指标,通过线性变换，化为少数相互独⽴的低维综合指标的
⼀种多元统计⽅法。

⼜称为主分量分析。

使⽤PCV降维时应该满⾜：
（1）使降维后样本的⽅差尽可能⼤。

（2）使降维后数据的均⽅误差尽可能⼩。

1.4.2 基本计算步骤
（1）计算给定样本{x n},n=1,2,...,N的均值mean_x和⽅差S；
（2）计算S的特征向量与特征值，X = UΛU T；
（3）将特征值从⼩到⼤排列，前M个特征值λ1，λ1，...,λM所对应的特征向量u1,u2,...,u M构成投影矩阵。

2 矩阵特征值及特征向量求解
2.1 定义
A为n阶矩阵，若数λ和n维⾮0列向量x满⾜Ax=λx，那么数λ称为A的特征值，x称为A的对应于特征值λ的特征向量。

式Ax=λx也可写成( A-λE)x=0，并且|λE-A|叫做A 的特征多项式。

|λE-A| = 0，称为A的特征⽅程（⼀个齐次线性⽅程组），求解特征值的过程其实就是求解特征⽅程的解。

2.2 特征值分解
特征值分解是将⼀个矩阵分解成下⾯的形式：
其中Q是这个矩阵A的特征向量组成的矩阵。

Σ = diag(λ1, λ2, ..., λn)是⼀个对⾓阵，每⼀个对⾓线上的元素就是⼀个特征值。

3 奇异值分解（SVD）
特征值分解只适⽤于⽅阵。

⽽在现实的世界中，我们看到的⼤部分矩阵都不是⽅阵，我们怎样才能像描述特征值⼀样描述这样⼀般矩阵的重要特征呢？奇异值分解（Singular Value Decomposition, 简称SVD）可⽤来解决这个问题。

SVD⽐特征值分解的使⽤范围更⼴，是⼀个适⽤于任⼀矩阵分解的⽅法。

假设X是⼀个m * n的矩阵，那么得到的U是⼀个m * m的⽅阵（⾥⾯的向量是正交的，U⾥⾯的向量称为左奇异向量）,
Σ是⼀个m * n的实数对⾓矩阵（对⾓线以外的元素都是0，对⾓线上的元素称为奇异值）,
V T(V的转置)是⼀个n * n的矩阵，⾥⾯的向量也是正交的，V⾥⾯的向量称为右奇异向量，
从图⽚来反映⼏个相乘的矩阵的⼤⼩可得下⾯的图⽚,
将⼀个矩阵X的转置 X T * X，将会得到 X T X 是⼀个⽅阵，我们⽤这个⽅阵求特征值可以得到：
(A T A)v i = λi v i
　这⾥得到的V，就是我们上⾯的右奇异向量。

这⾥的σi就是就是上⾯说的奇异值λi，u i就是上⾯说的左奇异向量。

常见的做法是将奇异值由⼤⽽⼩排列。

如此Σ便能由M唯⼀确定了。

奇异值σ跟特征值类似，在矩阵Λ中也是从⼤到⼩排列，⽽且σ的减少特别的快，在很多情况下，前10%甚⾄1%的奇异值的和就占了全部的奇异值之和的99%以上了。

也就是说，我们也可以⽤前r⼤的奇异值来近似描述矩阵，这⾥定义⼀下部分奇异值分解：
r是⼀个远⼩于m、n的数，这样矩阵的乘法看起来像是下⾯的样⼦：
右边的三个矩阵相乘的结果将会是⼀个接近于X的矩阵，在这⼉，r越接近于n，则相乘的结果越接近于A。

⽽这三个矩阵的⾯积之和（在存储观点来说，矩阵⾯积越⼩，存储量就越⼩）要远远⼩于原始的矩阵A，我们如果想要压缩空间来表⽰原矩阵A，我们存下这⾥的三个矩阵：U、Λ、V就好了。

4 特征分解与奇异值分解
4.1 低秩近似（Low-rank Approximation）
4.2 关联
5 pca(matrix)函数代码
1 def pca(X):
2 """ 主成分分析(PCA)
3 输⼊: X, （每⾏存储训练数据（⼀张图像）的扁平数组）
4 返回: 投影矩阵(V,重要维度优先)、⽅差(S)和均值(mean_X)
5 """
6
7 # 获得维度：num_data,dim = 图⽚数量 * 图⽚⼤⼩ ------>图⽚⼤⼩ = 长*宽
8 num_data,dim = X.shape
9
10 # 中⼼数据
11 mean_X = X.mean(axis=0) # 每列求平均值得到 1 * dim
12 X = X - mean_X # X.shape = num_data * dim
13
14 if dim > num_data:
15 # PCA - 紧凑的使⽤技巧
16 M = dot(X,X.T) # 协⽅差矩阵
17 e,EV = linalg.eigh(M) # e,EV = M的特征值，特征向量
18 tmp = dot(X.T,EV).T # 使⽤紧凑技巧？
19 V = tmp[::-1] # 为了得到最后的特征向量
20 S = sqrt(e)[::-1] # 反向，因为特征值是递增的
21 for i in range(V.shape[1]): # V.shape[1] = ？？？
22 V[:,i] /= S
23 else:
24 # PCA - 使⽤奇异值分解(SVD)
25 U,S,V = linalg.svd(X)
26 V = V[:num_data] # 只有返回前num_data个才有意义但是我认为此处不应该为V[:num_data]⽽应该为 V[:dim] ------>因为num_data >= dim
27
28 # 返回n the projection matrix, the variance and the mean
29 return V,S,mean_X
6 简单实例说明PCV实现的过程
6.1 列数⼤于⾏数时
1# 初始化矩阵X （列数⼤于⾏数时）
2 X = [[1 0 1 0 0 0]
3 [0 1 0 0 0 0]
4 [1 1 0 0 0 0]
5 [1 0 0 1 1 0]
6 [0 0 0 1 0 1]]
7#计算每列的平均值
8 mean_X =X.mean(axis=0)= [0.6 0.4 0.2 0.4 0.2 0.2]
9# 计算中⼼数据
10 X =X - mean_X = [[ 0.4 -0.4 0.8 -0.4 -0.2 -0.2]
11 [-0.6 0.6 -0.2 -0.4 -0.2 -0.2]
12 [ 0.4 0.6 -0.2 -0.4 -0.2 -0.2]
13 [ 0.4 -0.4 -0.2 0.6 0.8 -0.2]
14 [-0.6 -0.4 -0.2 0.6 -0.2 0.8]]
15# 计算协⽅差矩阵
16 M =dot(X, X.T)=[[ 1.20000000e+00 -4.00000000e-01 5.55111512e-17 -2.00000000e-01
17 -6.00000000e-01]
18 [-4.00000000e-01 1.00000000e+00 4.00000000e-01 -8.00000000e-01
19 -2.00000000e-01]
20 [ 5.55111512e-17 4.00000000e-01 8.00000000e-01 -4.00000000e-01
21 -8.00000000e-01]
22 [-2.00000000e-01 -8.00000000e-01 -4.00000000e-01 1.40000000e+00
23 0.00000000e+00]
24 [-6.00000000e-01 -2.00000000e-01 -8.00000000e-01 0.00000000e+00
25 1.60000000e+00]]
26# 计算协⽅差的特征值和特征向量
27 e, EV = numpy.linalg.eigh(M)
28 e = [-6.04093951e-16 2.66097365e-01 1.17478886e+00 2.00000000e+00
29 2.55911378e+00]
30 EV = [[ 4.47213595e-01 -1.38203855e-01 6.92423635e-01 5.00000000e-01
31 -2.26824169e-01]
32 [ 4.47213595e-01 -6.14411448e-01 -1.73966345e-01 -5.00000000e-01
33 -3.77139608e-01]
34 [ 4.47213595e-01 6.98980014e-01 -3.02870626e-01 4.38587603e-16
35 -4.68717745e-01]
36 [ 4.47213595e-01 -2.11286151e-01 -5.40988322e-01 5.00000000e-01
37 4.61183041e-01]
38 [ 4.47213595e-01 2.64921441e-01 3.25401658e-01 -5.00000000e-01
39 6.11498480e-01]]
40#⾏之间反转
41 tmp = dot(X.T, EV).T =[[-3.33066907e-16 -3.33066907e-16 4.16333634e-16 -3.33066907e-16
42 -3.88578059e-16 1.66533454e-16]
43 [ 3.49490007e-01 8.45685660e-02 -1.38203855e-01 5.36352895e-02
44 -2.11286151e-01 2.64921441e-01]
45 [-1.51435313e-01 -4.76836971e-01 6.92423635e-01 -2.15586664e-01
46 -5.40988322e-01 3.25401658e-01]
47 [ 1.00000000e+00 -5.00000000e-01 5.00000000e-01 -7.21644966e-16
48 5.00000000e-01 -5.00000000e-01]
49 [-2.34358872e-01 -8.45857353e-01 -2.26824169e-01 1.07268152e+00
50 4.61183041e-01 6.11498480e-01]]
51#默认返回的是按特征值升序，上下调转⼀下
52 V =tmp[::-1] = [[-2.34358872e-01 -8.45857353e-01 -2.26824169e-01 1.07268152e+00
53 4.61183041e-01 6.11498480e-01]
54 [ 1.00000000e+00 -5.00000000e-01 5.00000000e-01 -7.21644966e-16
55 5.00000000e-01 -5.00000000e-01]
56 [-1.51435313e-01 -4.76836971e-01 6.92423635e-01 -2.15586664e-01
57 -5.40988322e-01 3.25401658e-01]
58 [ 3.49490007e-01 8.45685660e-02 -1.38203855e-01 5.36352895e-02
59 -2.11286151e-01 2.64921441e-01]
60 [-3.33066907e-16 -3.33066907e-16 4.16333634e-16 -3.33066907e-16
61 -3.88578059e-16 1.66533454e-16]]
62#sigma = 特征值开根号
63#使特征值降序排序
64 S =sqrt(e)[::-1]= [1.59972303 1.41421356 1.08387677 0.51584626 nan]
65#特征向量归⼀化
66for i in range(V.shape[1]):
67 V[:, i] /= S
68 -----------------------------------------------------------------------------
69 V = [[-1.46499655e-01 -5.28752375e-01 -1.41789650e-01 6.70542025e-01
70 2.88289305e-01 3.82252720e-01]
71 [ 7.07106781e-01 -3.53553391e-01 3.53553391e-01 -5.10280049e-16
72 3.53553391e-01 -3.53553391e-01]
73 [-1.39716356e-01 -4.39936516e-01 6.38839814e-01 -1.98903298e-01
74 -4.99123458e-01 3.00220160e-01]
75 [ 6.77508074e-01 1.63941415e-01 -2.67916753e-01 1.03975338e-01
76 -4.09591321e-01 5.13566659e-01]
77 [ nan nan nan nan
78 nan nan]]
79 S = [1.59972303 1.41421356 1.08387677 0.51584626 nan]
6.2 列数不⼤于⾏数时
1 1 #初始化数组X
2 2 X = [[ 1 3]
3 3 [
4 7]
4 4 [ 3 2]
5 5 [ 1 23]]
6 6 num_data, dim = X.shape
7 7 X.shape= (4, 2)
8 8 U, S, V = linalg.svd(X)
9 9 U.shape, S.shape, V.shape=(4, 4) (2,) (2, 2)
10 10 V = [[ 0.10462808 0.99451142]
11 11 [ 0.99451142 -0.10462808]]
12 12 V = V[:dim] #V = V[:dim] or V =V[:num_data] ??? 注意：此时num_data > dim，所以我认为应该是V[:dim]
13 13 V[dim-1] =V[ 1 ]= [ 0.99451142 -0.10462808]
14 14 #计算结果
15 15 U = [[ 0.12635702 0.149642 -0.4064109 -0.89245244]
16 16 [ 0.30196809 0.71358512 -0.49810453 0.38923441]
17 17 [ 0.09422707 0.60994996 0.75324035 -0.22740114]
18 18 [ 0.94019702 -0.31042648 0.13910799 0.0177182 ]]
19 19 S = [24.43997403 4.54836997]
20 20 V = [[ 0.10462808 0.99451142]
21 21 [ 0.99451142 -0.10462808]]
7 环境配置
7.1 需要准备的安装包
下载python(x,y)2.7.x: /get/Programming/Other-Programming-Files/Python-x-y.shtml
下载PVC包：https:///willard-yuan/pcv-book-code
7.2 安装Python(x, y)
在Windows下，推荐安装Python(x,y)2.7.x。

Python(x,y)2.7.x是⼀个库安装包，除了Python⾃⾝外，还包含了第三⽅库，下⾯是安装Python(x, y)时的界⾯：
选择Full模式，⾃动安装的包见下图，
从上⾯第⼆幅图可以看出，pythonxy不仅包含了SciPy、NumPy、PyLab、OpenCV、MatplotLib,还包含了机器学习库scikits-learn。

为避免出现运⾏实例时出现的依赖问题，译者建议将上⾯的库全部选上，也就是选择“full”(译者也是⽤的全部安装的⽅式进⾏后⾯的实验的)。

安装完成后，为验证安装是否正确，可以在Python shell⾥确认⼀下OpenCV是否已安装来进⾏验证，在Python Shell⾥输⼊下⾯命令：
1from cv2 import__version__
2__version__
⾃动安装的pyopencv有问题。

需要你把opencv卸载后重新安装就可以了。

重新安装完成后，再次输⼊上⾯命令，就可以看到OpenCV的版本信息，则说明python(x,y)已安装正确。

7.2 安装PCV库
7.2.1 安装twine
cd 到对应python的scripts⽬录中，执⾏以下命令
1 pip install twine
或者
2 python
3 -m pip install twine
7.2.2 PVC的安装
PCV库是⼀个第三⽅库。

安装PVC之前必须安装twine，否则安装会出错。

假设你已从上⾯⽹址上下载了中译版源码，从Windows cmd终端
进⼊PCV所在⽬录：
1 cd PCV
2 python setup.py install
运⾏上⾯命令，即可完成PCV库的安装。

为了验证PCV库是否安装成功，在运⾏上⾯命令后，可以打开Python⾃带的Shell，在Shell输⼊：
1import PCV
如果未报错，则表明你已成功安装了该PCV库。

8 图像处理案案例
8.1 图像处理测试代码
1# -*- coding: utf-8 -*-
2import pickle
3from PIL import Image
4from numpy import *
5from pylab import *
6from PCV.tools import imtools, pca
7
8# 获取图像列表和他们的尺⼨
9 imlist = imtools.get_imlist('../data/fontimages/a_thumbs') # fontimages.zip is part of the book data set
10 im = array(Image.open(imlist[0])) # open one image to get the size
11 m, n = im.shape[:2] # get the size of the images
12 imnbr = len(imlist) # get the number of images
13print"The number of images is %d" % imnbr
14
15# Create matrix to store all flattened images
16 immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], 'f')
17
18# PCA降维
19 V, S, immean = pca.pca(immatrix)
20
21# 保存均值和主成分
22 f = open('../data/fontimages/font_pca_modes.pkl', 'wb')
23 pickle.dump(immean, f)
24 pickle.dump(V, f)
25 f.close()
26
27# Show the images (mean and 7 first modes)
28# This gives figure 1-8 (p15) in the book.
29 figure()
30 gray()
31 subplot(2, 4, 1)
32 axis('off')
33 imshow(immean.reshape(m, n))
34for i in range(7):
35 subplot(2, 4, i+2)
36 imshow(V[i].reshape(m, n))
37 axis('off')
38 show()
8.2 ⽂件⽬录相对位置
data与执⾏⽂件的相对位置：
待处理图⽚.jpg的位置
8.3 执⾏结果
参考资料：/pcvwithpython/installation.html。