PCA原理及实现

╰+哭是因爲堅強的太久メ 2022-08-05 04:23 180阅读 0赞

深度神经网路已经在语音识别，图像识别等领域取得前所未有的成功。本人在多年之前也曾接触过神经网络。本系列文章主要记录自己对深度神经网络的一些学习心得。

第五篇，谈谈PCA模型。本来PCA模型与深度学习是没有任何联系的。通常我们只是用PCA来对机器学习的数据做预处理。

本来想详细记录一下PCA的原理，但网上已经有一篇不错的文章，链接如下：

[http://hi.baidu.com/ifengzh/item/8851b6387aebefc4382ffa60][http_hi.baidu.com_ifengzh_item_8851b6387aebefc4382ffa60]

本文前面部分内容引用了这篇文章的内容。

## **一、简介** ##

PCA（Principal Components Analysis）即主成分分析，是图像处理中经常用到的降维方法，大家知道，我们在处理有关数字图像处理方面的问题时，比如经常用的图像的查询问题，在一个几万或者几百万甚至更大的数据库中查询一幅相近的图像。这时，我们通常的方法是对图像库中的图片提取响应的特征，如颜色，纹理，sift，surf，vlad等等特征，然后将其保存，建立响应的数据索引，然后对要查询的图像提取相应的特征，与数据库中的图像特征对比，找出与之最近的图片。这里，如果我们为了提高查询的准确率，通常会提取一些较为复杂的特征，如sift，surf等，一幅图像有很多个这种特征点，每个特征点又有一个相应的描述该特征点的128维的向量，设想如果一幅图像有300个这种特征点，那么该幅图像就有300\*vector（128维）个，如果我们数据库中有一百万张图片，这个存储量是相当大的，建立索引也很耗时，如果我们对每个向量进行PCA处理，将其降维为64维，是不是很节约存储空间啊？对于学习图像处理的人来说，都知道PCA是降维的，但是，很多人不知道具体的原理，为此，我写这篇文章，来详细阐述一下PCA及其具体计算过程：

## **二、PCA原理** ##

**1、原始数据：**

为了方便，我们假定数据是二维的，借助网络上的一组数据，如下：

x=\[2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2, 1,1.5, 1.1\]T  
y=\[2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9\]T

**2、计算协方差矩阵**

什么是协方差矩阵？相信看这篇文章的人都学过数理统计，一些基本的常识都知道，但是，也许你很长时间不看了，都忘差不多了，为了方便大家更好的理解，这里先简单的回顾一下数理统计的相关知识，当然如果你知道协方差矩阵的求法你可以跳过这里。

**（1）协方差矩阵：**

首先我们给你一个含有n个样本的集合，依次给出数理统计中的一些相关概念：

**均值： ![20130607154414484][]**

**标准差： ![20130607154423250][]**

**方差： ![20130607154443656][]**

既然我们都有这么多描述数据之间关系的统计量，为什么我们还要用协方差呢？我们应该注意到，标准差和方差一般是用来描述一维数据的，但现实生活我们常常遇到含有多维数据的数据集，最简单的大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解这几科成绩之间的关系，这时，我们就要用协方差，协方差就是一种用来度量两个随机变量关系的统计量，其定义为：

![20130607154452296][]

从协方差的定义上我们也可以看出一些显而易见的性质，如：

![20130607154503718][]

需要注意的是，协方差也只能处理二维问题，那维数多了自然就需要计算多个协方差，比如n维的数据集就需要计算CN2【此乃组合数基本公式】个协方差，那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义：

![20130607155311062][]

这个定义还是很容易理解的，我们可以举一个简单的三维的例子，假设数据集有三个维度\{x,y,z\}，则协方差矩阵为

![20130607155321562][]

可见，协方差矩阵是一个对称的矩阵，而且对角线是各个维度上的方差。

**（2）协方差矩阵的求法：**

**协方差矩阵计算的是不同维度之间的协方差，而不是不同样本之间的。下面我们将在matlab中用一个例子进行详细说明：**

首先，随机产生一个10\*3维的整数矩阵作为样本集，10为样本的个数，3为样本的维数。

**\[cpp\]** [view plain][] [copy][view plain]

1.  MySample = fix(rand(10,3)\*50)

![20130607155332203][]

根据公式，计算协方差需要计算均值，那是按行计算均值还是按列呢，我一开始就老是困扰这个问题。前面我们也特别强调了，协方差矩阵是计算**不同维度间**的协方差，要时刻牢记这一点。**样本矩阵的每行是一个样本，每列为一个维度**，所以我们要**按列计算均值**。为了描述方便，我们先将三个维度的数据分别赋值：

**\[cpp\]** [view plain][] [copy][view plain]

1.  dim1 = MySample(:,1); 
2.  dim2 = MySample(:,2); 
3.  dim3 = MySample(:,3); 
4.  %计算dim1与dim2，dim1与dim3，dim2与dim3的协方差： 
5.  sum( (dim1-mean(dim1)) .\*(dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) %得到 74.5333 
6.  sum( (dim1-mean(dim1)) .\* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10.0889 
7.  sum( (dim2-mean(dim2)) .\* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10\*\*\*000 
8.  %搞清楚了这个后面就容易多了，协方差矩阵的对角线就是各个维度上的方差，下面我们依次计算： 
9.  std(dim1)^2 % 得到 108.3222 
10. std(dim2)^2 % 得到 260.6222 
11. std(dim3)^2 % 得到 94.1778 
12. %这样，我们就得到了计算协方差矩阵所需要的所有数据，调用Matlab自带的cov函数进行验证： 
13. cov(MySample)

![20130607160112921][]

可以看到跟我们计算的结果是一样的，说明我们的计算是正确的。但是通常我们不用这种方法，而是用下面简化的方法进行计算：

先让样本矩阵中心化，即每一维度减去该维度的均值，然后直接用新的到的样本矩阵乘上它的转置，然后除以(N-1)即可。其实这种方法也是由前面的公式通道而来，只不过理解起来不是很直观而已。大家可以自己写个小的矩阵看一下就明白了。其Matlab代码实现如下：

**\[cpp\]** [view plain][] [copy][view plain]

1.  X = MySample –repmat(mean(MySample),10,1); %中心化样本矩阵 
2.  C = (X’\*X)./(size(X,1)-1) 
3.  %为方便对matlab不太明白的人，小小说明一下各个函数，同样，对matlab有一定基础的人直接跳过： 
4.  %B = repmat(A,m,n ) %%将矩阵 A复制 m×n块，即把 A 作为 B的元素，B由 m×n个 A平铺而成。B的维数是 \[size(A,1)\*m, (size(A,2)\*n\] 
5.  %B = mean(A)的说明： 
6.  %如果你有这样一个矩阵：A = \[1 2 3; 3 36; 4 6 8; 4 7 7\]; 
7.  %用mean(A)（默认dim=1）就会求每一列的均值 
8.  % ans = 
9.  % 3.0000 4.5000 6.0000 
10. % 用mean(A,2)就会求每一行的均值 
11. % ans = 
12. % 2.0000 
13. % 4.0000 
14. % 6.0000 
15. % 6.0000 
16. size(A,n)%% 如果在size函数的输入参数中再添加一项n，并用1或2为n赋值，则 size将返回矩阵的行数或列数。其中r=size(A,1)该语句返回的是矩阵A的行数， %c=size(A,2)该语句返回的是矩阵A的列数

上面我们简单说了一下协方差矩阵及其求法，言归正传，我们用上面简化求法，求出样本的协方差矩阵为：

![20130607155340718][]

**3、计算协方差矩阵的特征向量和特征值**

因为协方差矩阵为方阵，我们可以计算它的特征向量和特征值，如下：

**\[cpp\]** [view plain][] [copy][view plain]

1.  \[eigenvectors,eigenvalues\] = eig(cov)

![20130607155347578][]

我们可以看到这些矢量都是单位矢量，也就是它们的长度为1，这对PCA来说是很重要的。

**4、选择成分组成模式矢量**

求出协方差矩阵的特征值及特征向量之后，按照特征值由大到小进行排列，这将给出成分的重要性级别。现在，如果你喜欢，可以忽略那些重要性很小的成分，当然这会丢失一些信息，但是如果对应的特征值很小，你不会丢失很多信息。如果你已经忽略了一些成分，那么最后的数据集将有更少的维数，精确地说，如果你的原始数据是n维的，你选择了前p个主要成分，那么你现在的数据将仅有p维。现在我们要做的是组成一个模式矢量，这只是几个矢量组成的矩阵的一个有意思的名字而已，它由你保持的所有特征矢量构成，每一个特征矢量是这个矩阵的一列。

对于我们的数据集，因为有两个特征矢量，因此我们有两个选择。我们可以用两个特征矢量组成模式矢量：

![20130607160147437][]

我们也可以忽略其中较小特征值的一个特征矢量，从而得到如下模式矢量：

![20130607160200875][]

**5、得到降维后的数据**

![20130607160300656][]

其中rowFeatureVector是由模式矢量作为列组成的矩阵的转置，因此它的行就是原来的模式矢量，而且对应最大特征值的特征矢量在该矩阵的最上一行。rowdataAdjust是每一维数据减去均值后，所组成矩阵的转置，即数据项目在每一列中，每一行是一维，对我们的样本来说即是，第一行为x维上数据，第二行为y维上的数据。FinalData是最后得到的数据，数据项目在它的列中，维数沿着行。

这将给我们什么结果呢？这将仅仅给出我们选择的数据。我们的原始数据有两个轴（x和y），所以我们的原始数据按这两个轴分布。我们可以按任何两个我们喜欢的轴表示我们的数据。如果这些轴是正交的，这种表达将是最有效的，这就是特征矢量总是正交的重要性。我们已经将我们的数据从原来的xy轴表达变换为现在的单个特征矢量表达。

说明：如果要恢复原始数据，只需逆过程计算即可，即：

![20130607160322515][]

![20130607160352968][]

到此为止，相信你已经掌握了PCA的原理了。

## 三 . PCA的应用 ##

PCA及其改进算法主要应用的人脸识别领域，是人脸识别的经典算法之一。OpenCv2.4以后的版本实现了三种经典的人脸识别算法，其中就包括PCA。对openCv比较老的版本也可以调用PCA的算法去做，只是稍显复杂而已，网上有一篇博文如下：

[http://www.cognotics.com/opencv/servo\_2007\_series/part\_5/index.html][http_www.cognotics.com_opencv_servo_2007_series_part_5_index.html]

该代码运行在openCv2.1之前的版本当中，但是该代码有个重要的bug就是特征数K被设置为固定的值，而选择更小的值的时候，代码会crash。

PCA另外一个主要的用途是作为其他算法的预处理，术语叫做数据的白化。由于PCA具有压缩数据的作用，所以可以认为经过PCA处理过之后的数据是不相关的，但一般未必是独立的。实际可用的PCA算法一般不是以解析解的形式给出的，而是在线学习算法。有很多的原因决定了只能使用在线学习算法。在线学习算法主要有基于神经网络学习的算法和递归最小二乘法，相关的文献如下：

[http://wenku.baidu.com/view/c91f31c058f5f61fb73666f8.html][http_wenku.baidu.com_view_c91f31c058f5f61fb73666f8.html]

要注意的是openCv的实现不是在线学习算法。

## 四. PCA 的实现 ##

前面已经谈到了PCA的实现分为解析解和在线学习算法。解析解适合于数据量小并且数据完全已知的情况下，这里给出一种高效的解析解的实现代码。

4.1 数据结构定义及API说明如下：

**\[cpp\]** [view plain][] [copy][view plain]

1.  \#ifndef \_FCE\_\_PCA\_\_H\_\_
2.  \#define \_FCE\_\_PCA\_\_H\_\_
3.  
4.  \#define HIGH\_PRECISON
5.  
6.  \#ifdef HIGH\_PRECISON
7.  \#define real float
8.  \#else
9.  \#define real double
10. \#endif
11. 
12. 
13. \#ifdef \_cplusplus
14. \{ 
15. **extern**"C"
16. \#endif
17. 
18. 
19. **typedef****struct** \_FCE\_PCA\{ 
20. **int** count; //the number of sample
21. **int** n; // the number of features
22.  real \*covariance; 
23.  real \*mean; 
24.  real \*z; 
25. \}FCE\_PCA; 
26. 
27. 
28. FCE\_PCA \*fce\_pca\_init(**int** n); 
29. 
30. **void** fce\_pca\_push\_add(FCE\_PCA \*pca, real \*v); 
31. 
32. **int** fce\_pca\_solve\_eig(FCE\_PCA \*pca, real \*eigenvector, real \*eigenvalue); 
33. 
34. **void** fce\_pca\_free(FCE\_PCA \*pca); 
35. 
36. \#ifdef \_cplusplus
37. \} 
38. \#endif
39. 
40. \#endif

函数 fce\_pca\_push\_add 用于把每一个样本点添加到PCA模型之中，例如，一个人脸的样本数据。

函数 fce\_pca\_solve\_eig 采用雅克比迭代法快速求解对称矩阵的特征值和特征向量，其它两个函数分别用以创建PCA模型和释放PCA模型。

4.2 各函数的实现

**\[cpp\]** [view plain][] [copy][view plain]

1.  \#include "fce\_pca.h"
2.  
3.  
4.  \#define FCE\_MIN(i,j) (((i) > (j)) ? (j) : (i))
5.  \#define FCE\_MAX(i,j) (((i) > (j)) ? (i) : (j))
6.  
7.  FCE\_PCA \*fce\_pca\_init(**int** n)\{ 
8.   FCE\_PCA \*pca; 
9.   real zero = 0.0; 
10. **if**(n <= 1) 
11. **return** NULL; 
12. 
13.  pca = (FCE\_PCA\* )malloc(**sizeof**(FCE\_PCA)); 
14. **if** (pca == NULL)\{ 
15. **return** NULL; 
16.  \} 
17. 
18.  pca->n = n; 
19.  pca->z = (real\* )malloc(**sizeof**(\*pca->z) \* n); 
20. **if** (pca->z == NULL)\{ 
21.  free(pca); 
22. **return** NULL; 
23.  \} 
24. 
25.  memset(pca->z, zero, **sizeof**(\*pca->z) \* n); 
26. 
27.  pca->count=0; 
28.  pca->covariance= (real\* )malloc(**sizeof**(real) \* n \* n); 
29. **if** (pca->covariance == NULL)\{ 
30.  free(pca->z); 
31.  free(pca); 
32. **return** NULL; 
33.  \} 
34. 
35.  memset(pca->covariance, zero, **sizeof**(real) \* n \* n); 
36. 
37.  pca->mean = (real\* )malloc(**sizeof**(real) \* n); 
38. **if** (pca->mean == NULL)\{ 
39.  free(pca->covariance); 
40.  free(pca->z); 
41.  free(pca); 
42. **return** NULL; 
43.  \} 
44. 
45.  memset(pca->mean, zero, **sizeof**(real) \* n); 
46. 
47. **return** pca; 
48. \} 
49. 
50. **void** fce\_pca\_free(FCE\_PCA \*pca)\{ 
51.  free(pca->covariance); 
52.  free(pca->mean); 
53.  free(pca->z); 
54.  free(pca); 
55. \} 
56. 
57. **void** fce\_pca\_push\_add(FCE\_PCA \*pca, real \*v)\{ 
58. **int** i, j; 
59. **const****int** n = pca->n; 
60. **for**(i = 0; i < n; i++)\{ 
61.  pca->mean\[i\] += v\[i\]; 
62. **for**(j = i; j < n; j++) 
63.  pca->covariance\[j + i \* n\] += v\[i\]\*v\[j\]; 
64.  \} 
65.  pca->count++; 
66. \} 
67. 
68. **int** fce\_pca\_solve\_eig(FCE\_PCA \*pca, real \*eigenvector, real \*eigenvalue)\{ 
69. **int** i, j, pass; 
70. **int** k = 0; 
71. **const****int** n = pca->n; 
72.  real \*z = pca->z; 
73.  real zero = 0.0; 
74. 
75.  memset(eigenvector, zero, **sizeof**(real)\*n\*n); 
76. 
77. **for**(j = 0; j < n; j++)\{ 
78.  pca->mean\[j\] /= pca->count; 
79.  eigenvector\[j + j \* n\] = 1.0; 
80. **for**(i = 0; i <= j; i++)\{ 
81.  pca->covariance\[j + i \* n\] /= pca->count; 
82.  pca->covariance\[j + i \* n\] -= pca->mean\[i\] \* pca->mean\[j\]; 
83.  pca->covariance\[i + j \* n\] = pca->covariance\[j + i \* n\]; 
84.  \} 
85.  eigenvalue\[j\] = pca->covariance\[j + j\*n\]; 
86.  z\[j\] = 0; 
87.  \} 
88. 
89. **for**(pass=0; pass < 50; pass++)\{ 
90.  real sum = 0; 
91. **for**(i = 0; i < n; i++) 
92. **for**(j = i+1; j < n; j++) 
93.  sum += fabs(pca->covariance\[j + i \* n\]); 
94. 
95. **if**(sum == 0)\{ 
96. **for**(i = 0; i < n; i++)\{ 
97.  real maxvalue = -1; 
98. **for**(j = i; j < n; j++)\{ 
99. **if**(eigenvalue\[j\] > maxvalue)\{ 
100.  maxvalue = eigenvalue\[j\]; 
101.  k= j; 
102.  \} 
103.  \} 
104.  eigenvalue\[k\] = eigenvalue\[i\]; 
105.  eigenvalue\[i\] = maxvalue; 
106. **for**(j = 0; j < n; j++)\{ 
107.  real tmp = eigenvector\[k + j \* n\]; 
108.  eigenvector\[k + j \* n\] = eigenvector\[i + j \* n\]; 
109.  eigenvector\[i + j \* n\] = tmp; 
110.  \} 
111.  \} 
112. **return** pass; 
113.  \} 
114. 
115. **for**(i = 0; i < n; i++)\{ 
116. **for**(j = i + 1; j < n; j++)\{ 
117.  real covar = pca->covariance\[j + i \* n\]; 
118.  real t,c,s,tau,theta, h; 
119. 
120. **if**(pass < 3 && fabs(covar) < sum / (5\*n\*n)) 
121. **continue**; 
122. **if**(fabs(covar) <= 0.00000000001) 
123. **continue**; 
124. **if**(pass >=3 && fabs((eigenvalue\[j\]+z\[j\])/covar) > (1LL<<32) && fabs((eigenvalue\[i\]+z\[i\])/covar) > (1LL<<32))\{ 
125.  pca->covariance\[j + i \* n\]=0.0; 
126. **continue**; 
127.  \} 
128. 
129.  h = (eigenvalue\[j\] + z\[j\]) - (eigenvalue\[i\] + z\[i\]); 
130.  theta = 0.5 \* h/covar; 
131.  t = 1.0 /(fabs(theta) + sqrt(1.0 + theta \* theta)); 
132. **if**(theta < 0.0) t = -t; 
133. 
134.  c = 1.0 /sqrt(1 + t \* t); 
135.  s = t \* c; 
136.  tau = s /(1.0 + c); 
137.  z\[i\] -= t \* covar; 
138.  z\[j\] += t \* covar; 
139. 
140. \#define ROTATE(a,i,j,k,l) \{\\
141.  real g =a\[j + i\*n\];\\ 
142.  real h =a\[l + k\*n\];\\ 
143.  a\[j + i\*n\] = g - s \* (h + g \* tau);\\ 
144.  a\[l + k\*n\] = h + s \* (g - h \* tau); \} 
145. **for**(k = 0; k < n; k++) \{ 
146. **if**(k != i && k != j)\{ 
147.  ROTATE(pca->covariance,FCE\_MIN(k,i),FCE\_MAX(k,i),FCE\_MIN(k,j),FCE\_MAX(k,j)) 
148.  \} 
149.  ROTATE(eigenvector,k,i,k,j) 
150.  \} 
151.  pca->covariance\[j + i \* n\]=0.0; 
152.  \} 
153.  \} 
154. **for** (i = 0; i < n; i++) \{ 
155.  eigenvalue\[i\] += z\[i\]; 
156.  z\[i\]=0.0; 
157.  \} 
158.  \} 
159. 
160. **return** 0; 
161. \}

[http_hi.baidu.com_ifengzh_item_8851b6387aebefc4382ffa60]: http://hi.baidu.com/ifengzh/item/8851b6387aebefc4382ffa60
[20130607154414484]: /images/20220805/f252299274fe4269b4b44bd059a25d5b.png
[20130607154423250]: /images/20220805/0edc984e088c4b4dac4c1bb0718155b1.png
[20130607154443656]: /images/20220805/d08c4acd46d244dcb8a95b66e237ef2a.png
[20130607154452296]: /images/20220805/c75bdd2257734de68cab8fdf71417b68.png
[20130607154503718]: /images/20220805/ddf146c1d1f642f5a75dab660dd177be.png
[20130607155311062]: /images/20220805/30e17b6b1b0140b1812258300bfe9711.png
[20130607155321562]: /images/20220805/7d1d41c25285433a87e8b17c81ff0c65.png
[view plain]: http://blog.csdn.net/celerychen2009/article/details/9048033#
[20130607155332203]: /images/20220805/cca09e31189e4f61842732a3e61ee961.png
[20130607160112921]: /images/20220805/1ae018650e2941ce90f2ba21bb662cdf.png
[20130607155340718]: /images/20220805/e23487f01b8e4aee826321fd0c102d88.png
[20130607155347578]: /images/20220805/d3f75671a9cd4629980c246f39f4136f.png
[20130607160147437]: /images/20220805/5801270cbd034697a18f4e547c3ef88b.png
[20130607160200875]: /images/20220805/3b59e001a77a4884842453b46095c699.png
[20130607160300656]: /images/20220805/526d21e95e14427a964a49003cbdf8c0.png
[20130607160322515]: /images/20220805/cbd82243ab1e49cc9e93771ec35d502d.png
[20130607160352968]: /images/20220805/371004f68a0049deac2f97c79a9e3e4a.png
[http_www.cognotics.com_opencv_servo_2007_series_part_5_index.html]: http://www.cognotics.com/opencv/servo_2007_series/part_5/index.html
[http_wenku.baidu.com_view_c91f31c058f5f61fb73666f8.html]: http://wenku.baidu.com/view/c91f31c058f5f61fb73666f8.html