【机器学习】 PCA(Principal Component Analysis)——主成分分析
PCA方法原理
首先,PCA是做什么的。一般采集到的原始数据中,数据量很大,这导致数据不容易分析,但是有时候数据中潜在的关键信息其实很少,这就需要PCA将数据简化,去掉冗余的没有用的信息,留下关键的有用信息,便于分析数据。
举例说明一下PCA的作用:
假设有上图中这样一个运动系统,忽略弹簧的质量和所有摩擦力,给小球一个x轴方向的力之后,小球会沿着x轴不停的往复运动。现在我们只能观察到这个小球在运动,它为什么运动,沿着什么方向做什么样的运动我们事先都不清楚。所以我们在
A
、
B
、
C
A、B、C
A、B、C三个位置放了三个
120
H
z
120Hz
120Hz的相机拍摄小球的运动。因为不知道小球做什么样的运动,自然事先也没有坐标轴的概念,所以三个相机的位置是随机的。三个摄像机的拍摄结果如下,这就是我们采集到的数据:
上帝视角里,我们知道小球的运动就是沿着x轴的简单的往复运动,但是我们采集到的数据因为有冗余信息且关键信息不突出,我们无法高效率的分析它并得到结论。所以,我们需要PCA。PCA能去掉上面数据中冗余的数据,并突出关键数据,甚至可以根据情况进行降维和数据压缩(保留主成分的数据),而且不会损失太多的信息。
接下来逐步推导PCA的过程:
- 每个样本是 m m m维的列向量。在上面的例子中,每个相机拍摄的数据有 x , y x,y x,y两个轴的信息,三个相机也就是每个数据是一个 6 6 6维的列向量,如果相机拍摄了10分钟,得到的数据一共有 10 ∗ 60 ∗ 120 10*60*120 10∗60∗120,也就是 72000 72000 72000组。将这些数据按列排列,得到 6 ∗ 72000 6*72000 6∗72000的数据矩阵,泛化开是 m ∗ n m*n m∗n的矩阵, m m m是数据为维度, n n n是数据的数量。我们设它为 X X X。
- 这里出现了PCA的第一个假设,线性假设。我们每一个数据可以用一组基向量的线性组合表示:
上图中就是一组最简单的基向量。比如在上面小球的例子中,所有的数据可以用 x x x轴和 y y y轴的单位向量这组基向量的线性组合来表示。
P X = Y PX=Y PX=Y, X X X是上面提到的数据集, Y Y Y是通过线性变换 P X = Y PX=Y PX=Y得到的数据集的新的表示形式。 X X X和 Y Y Y都是 m ∗ n m*n m∗n的矩阵, P P P是一个 m ∗ m m*m m∗m的矩阵。
上图是 P X = Y PX=Y PX=Y拆开写的计算细节。计算得到的 Y Y Y每一列的表示如下:
p i p_i pi是 1 ∗ m 1*m 1∗m的向量, x i x_i xi是 m ∗ 1 m*1 m∗1的向量, p i ∗ x i p_i*x_i pi∗xi是一个数。 x i x_i xi能提取出来,所以 p 1 到 p m p_1到p_m p1到pm是一组新的基向量,数据集 Y Y Y中的每一个都可以用这组基向量的线性组合表示。理解这点很重要,因为后面我们知道 P P P的每一行和 X X X相乘得到的 1 ∗ n 1*n 1∗n的行向量代表数据不同重要程度的成分,而这里也解释了为什么可以只保留主成分的数据达到降维的效果。 - 数据中的噪声没有办法绝对标度,一般用SNR(signal-to-noise ratio)表示:
SNR越大,数据越干净。SNR越小,越是有噪声的数据。
A、B是数据的两个评价标准,从 1 1 1到 n n n代表数据集中的 n n n个样本。也就是前面说一个样本可以用一个 m m m维列向量表示, n n n个样本组成 m ∗ n m*n m∗n的矩阵,A、B就是这个矩阵的两个行向量。这里先约定A、B的均值都是0。
A、B各自的方差如上。
A、B的协方差如上。
重要概念:协方差是两个变量之间线性关系的程度,协方差为正数,越大,代表正相关程度越大。协方差为负数,越小,代表负相关程度越大。协方差为零,代表不相关。协方差绝对值越大,代表相关程度越大,也就是冗余程度越大。
将上面的评价标准A、B写成如a、b的行向量的形式,则a和b的协方差恒等于a与b转置相乘的结果除以n,如下:
我们采集到的数据集 X X X可以按行向量表示为如下形式:
则 X X X的协方差矩阵有如下关系:
可知协方差矩阵对角线上是方差,对角线之外是协方差。协方差反应的是数据中的中的噪声和冗余。这里出现PCA的第二个假设,协方差矩阵对角线上的值越大,对应的评价标准越重要。协方差非对角线的值越大,这两个评价标准越冗余。- 还记得前面所说的 P X = Y PX=Y PX=Y, Y Y Y是原始数据集 X X X的新的表现形式。那么可以通过 P C x = C y PC_x=C_y PCx=Cy,用 C y C_y Cy作为 C x C_x Cx的新的表现形式。那么我们的 C y C_y Cy有什么期望呢,肯定是最小化噪声和冗余,最大化关键信号啊。那么 C y C_y Cy就应该一个对角矩阵(也就是 Y Y Y是非相关的,为什么,因为 C y 恒 等 于 1 n Y Y T C_y恒等于\dfrac{1}{n}YY^T Cy恒等于n1YYT,只有 Y Y Y是非相关的时候, C y C_y Cy才是对角矩阵),且 Y Y Y按行算每个维度应该由其对应的方差大小,由大到小排序。
- 让 C y C_y Cy对角化的方法有很多,PCA选择了最简单的方法—— P P P是一个正交矩阵。这也是PCA的第三个假设。
- 为什么
P
P
P是正交矩阵,
C
y
C_y
Cy就能对角化呢。
上面这段推导没有什么好解释的,其中 ( A B ) T = B T A T (AB)^T=B^TA^T (AB)T=BTAT,这是矩阵转置的性质。
C x C_x Cx是什么样的矩阵, C x C_x Cx是一个实对称矩阵。实对称矩阵有什么重要性质呢:
上面我们说假设 P P P是正交矩阵。满足 A A T = E AA^T=E AAT=E的矩阵 A A A是正交矩阵,同时正交矩阵有性质 A T A^T AT也是正交矩阵,且 A T = A − 1 A^T=A^{-1} AT=A−1。
那么 C y = P ( C x ) P T C_y=P(C_x)P^T Cy=P(Cx)PT可以写成 C y = P T ( C x ) P C_y=P^T(C_x)P Cy=PT(Cx)P,然后再写成 C y = P − 1 C x P C_y=P^{-1}C_xP Cy=P−1CxP。所以当 P P P是某一正交矩阵时, C y C_y Cy是对角矩阵。 - 关键步骤,
P
P
P怎么求啊。线性代数里有!
所以怎么求主成分 P P P呢。
一、 X X X每行减去这一行的均值。
二、用 C x = 1 n X X T C_x=\dfrac{1}{n}XX^T Cx=n1XXT计算 X X X的协方差矩阵 C x C_x Cx。
三、计算 C x C_x Cx的所有特征值和对应的特征向量。
四、将特征向量单位化,按特征值由大到小顺序排序,组成主成分 P P P。
注意特征向量原本是一个
n
∗
1
n * 1
n∗1的列向量,但是在这里我们按行保存。
P
P
P的第一行是最重要的成分的基向量,第二行是第二重要的成分的基向量,以此类推。所以只用
P
P
P的第一行乘
X
X
X得到的
1
∗
n
1*n
1∗n的
Y
Y
Y就是只保留主成分的降维后的数据。
10. 从别人文章里copy过来一个例子,结合上面的解释理解一下。
PCA在人脸识别中的应用
先下载yalefaces数据集
因为数据集还没有处理,接下来我们用Haar+AdaBoost+cascade的方法对有数据集做人脸检测
详见【机器学习】Haar特征+AdaBoost+cascade实现人脸检测
import cv2
import glob
filename = "/home/zhangchen/yalefaces/*"
dataset_list = glob.glob(filename)
def detect(filename):
face_cscade = cv2.CascadeClassifier(
"./cascades/haarcascade_frontalface_default.xml")
save_name = filename.split(".")[0][-2:]+"_"+filename.split(".")[-1]+".png"
gif = cv2.VideoCapture(filename)
# 因为yalefaces数据集中数据的格式是.gif,所以要用cv2.VideoCapture()获取一张图片,如果用cv2.imread()会报错
_, img = gif.read()
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
face = face_cscade.detectMultiScale(img_gray, 1.3, 5)
for (x, y, w, h) in face:
save_img = cv2.resize(img_gray[y:y+h, x:x+w], (200, 200))
cv2.imwrite("/home/zhangchen/my_yalefaces/{}".format(save_name), save_img)
break
# 因为知道一张图里只有一个脸,所以保存之后break
for each in dataset_list:
detect(each)
处理后得到数据集如下:
一张灰度图就是一个数据,把
200
∗
200
200*200
200∗200的灰度图变成
40000
∗
1
40000*1
40000∗1的列向量,再把所有的列向量按列排列,就得到PCA中的
X
X
X。
yalefaces数据集中一共有164张图,15个人,我随机每个人拿出两张作为测试用,一共拿出30张图,还剩134张图。所以
X
X
X是一个
40000
∗
134
40000*134
40000∗134的矩阵。
接下来,计算每一行的平均值,再每一行的每一个元素分别减去对应的行的平均值。这些平均值是
40000
40000
40000个数值,恢复成
200
∗
200
200*200
200∗200的矩阵,得到平均脸。
import numpy as np
import cv2
dataset_dir = "/home/zhangchen/my_yalefaces/*"
dataset_list = glob.glob(dataset_dir)
def face_average(path_list):
X = np.zeros((40000, 1))
for each in path_list:
img = cv2.imread(each, 0)
x_list = np.split(img, 200, axis=1)
# opencv读取的图片是hxwxc的,这里按列切开,所以axis选择1
x = np.concatenate(x_list, axis=0)
X = np.concatenate((X, x), axis=1)
X = np.split(X, [1], axis=1)[1]
# 到这里得到了X
means = np.mean(X, axis=1)
# 计算每一行的平均值,这里得到了每一行的均值
means_list = np.split(means, 200, axis=0)
for i, each in enumerate(means_list):
means_list[i] = each[:, np.newaxis]
mean_face = np.concatenate(means_list, axis=1).astype(np.uint8)
while True:
cv2.imshow("mean_face", mean_face)
if cv2.waitKey(1000 // 12) & 0xff == ord("q"):
break
cv2.destroyAllWindows()
face_average(dataset_list)
可视化后得到的平均脸如下:
求得 X X X后,根据公式 C x = 1 n X X T C_x=\dfrac{1}{n}XX^T Cx=n1XXT计算 X X X的协方差矩阵 C x C_x Cx, X X X是一个 40000 ∗ 134 40000*134 40000∗134的矩阵,那么 C X C_X CX是 40000 ∗ 40000 40000 * 40000 40000∗40000的矩阵。然后计算 C x C_x Cx的所有特征值和对应的特征向量,每个特征向量是一个 40000 ∗ 1 40000 * 1 40000∗1的列向量,将特征向量单位化,按特征值由大到小顺序排序,按行存放,组成主成分 P P P。
但是大多数情况 X X T XX^T XXT的维度太大,比如我们这里就是一个 40000 ∗ 40000 40000*40000 40000∗40000的矩阵,你可以用python算一下,计算它的特征向量会死机。所以这个时候利用数学的小技巧:先求 X T X X^TX XTX的特征向量矩阵 V V V,其中 V V V的每一列是一个特征向量, V V V是一个 n ∗ n n*n n∗n的矩阵,然后从 V V V中取出前 r r r个最大特征值对应的特征向量,所以 V V V就变成了 n ∗ r n*r n∗r矩阵,然后 C x ‘ = X V C_x^`=XV Cx‘=XV, C x ‘ C_x^` Cx‘的转置就是计算出的 P P P,它是一个 r ∗ m r*m r∗m的矩阵,每一行代表不同重要成分。
每一个特征向量数值变化到0到255,再恢复成 200 ∗ 200 200 * 200 200∗200的矩阵就是一个个特征脸。
//为了方便在上面的代码上进行补充
import numpy as np
import cv2
dataset_dir = "/home/zhangchen/my_yalefaces/*"
dataset_list = glob.glob(dataset_dir)
def face_average(path_list):
X = np.zeros((40000, 1))
for each in path_list:
img = cv2.imread(each, 0)
x_list = np.split(img, 200, axis=1)
# opencv读取的图片是hxwxc的,这里按列切开,所以axis选择1
x = np.concatenate(x_list, axis=0)
X = np.concatenate((X, x), axis=1)
X = np.split(X, [1], axis=1)[1]
# 到这里得到了X
means = np.mean(X, axis=1)
# 计算每一行的平均值,这里得到了每一行的均值
means_list = np.split(means, 200, axis=0)
for i, each in enumerate(means_list):
means_list[i] = each[:, np.newaxis]
mean_face = np.concatenate(means_list, axis=1).astype(np.uint8)
while True:
cv2.imshow("mean_face", mean_face)
if cv2.waitKey(1000 // 12) & 0xff == ord("q"):
break
cv2.destroyAllWindows()
# 补充部分
X_list = np.split(X, 40000, axis=0)
# X按行分开,分成40000行
each_mean = np.split(means, 40000, axis=0)
# mean按行分开,分成40000行
for i in range(len(X_list)):
X_list[i] = X_list[i] - each_mean[i]
# 按行做减
X_new = np.concatenate(X_list, axis=0)
# 按行拼接回40000×134的矩阵,这是减过均值的X
C = np.mat(np.dot(X_new.T, X_new))
value, vector = np.linalg.eigh(C)
# 注意np.linalg.eigh是针对对称矩阵的,np.linalg.eig是针对非对称矩阵的
# value是特征值,vector是特征向量,注意vector的每一列是一个特征向量
V = []
# 这里存特征向量
for each in range(40):
# 我们这里保留40个特征向量
index = value.tolist().index(max(value.tolist()))
V.append(vector[index])
value[index] = -float("inf")
V = np.squeeze(np.array(V)).swapaxes(0, 1)
Cx = np.dot(X_new, V)
# 这个矩阵就是P
vector_list = np.split(Cx, 40, axis=1)
for i, each in enumerate(vector_list):
face = np.concatenate(np.split(each, 200, axis=0), axis=1).astype(np.uint8)
while True:
cv2.imshow("face", face)
if cv2.waitKey(1000 // 12) & 0xff == ord("q"):
break
cv2.destroyAllWindows()
cv2.imwrite("/home/zhangchen/face_{}.png".format(i), face)
# 补充部分
face_average(dataset_list)
得到特征脸如下:
一共生成了四十个特征脸。
上面说过,
P
X
=
Y
PX=Y
PX=Y,
p
1
到
p
m
p_1到p_m
p1到pm是一组新的基向量,数据集
Y
Y
Y中的每一个都可以用这组基向量的线性组合表示。所以人脸数据集
X
X
X中的每一个,都可以用这些特征脸的线性组合表示。
关键问题来了:那来一张新脸怎么做人脸识别呢。
特征矩阵
P
P
P是
40
∗
40000
40*40000
40∗40000 的矩阵,每一张人脸是
40000
∗
1
40000*1
40000∗1的向量,特征矩阵和人脸向量做矩阵乘法是一个
40
∗
1
40*1
40∗1的向量。将数据集中每一张人脸和特征矩阵做矩阵乘法,得到134个这样的向量,待识别的人脸也和特征矩阵做矩阵乘法得到一个这样的向量。分别将数据集中134个这样的向量和待识别人脸得到的这样的向量求欧式距离,我们认为最小的那个数据集中的人脸和待识别人脸是同一个人。
# 为了方便在上面的代码上进行补充
import numpy as np
import cv2
dataset_dir = "/home/zhangchen/my_yalefaces/*"
dataset_list = glob.glob(dataset_dir)
def face_average(path_list):
X = np.zeros((40000, 1))
for each in path_list:
img = cv2.imread(each, 0)
x_list = np.split(img, 200, axis=1)
# opencv读取的图片是hxwxc的,这里按列切开,所以axis选择1
x = np.concatenate(x_list, axis=0)
X = np.concatenate((X, x), axis=1)
X = np.split(X, [1], axis=1)[1]
# 到这里得到了X
means = np.mean(X, axis=1)
# 计算每一行的平均值,这里得到了每一行的均值
means_list = np.split(means, 200, axis=0)
for i, each in enumerate(means_list):
means_list[i] = each[:, np.newaxis]
mean_face = np.concatenate(means_list, axis=1).astype(np.uint8)
while True:
cv2.imshow("mean_face", mean_face)
if cv2.waitKey(1000 // 12) & 0xff == ord("q"):
break
cv2.destroyAllWindows()
X_list = np.split(X, 40000, axis=0)
# X按行分开,分成40000行
each_mean = np.split(means, 40000, axis=0)
# mean按行分开,分成40000行
for i in range(len(X_list)):
X_list[i] = X_list[i] - each_mean[i]
# 按行做减
X_new = np.concatenate(X_list, axis=0)
# 按行拼接回40000×134的矩阵,这是减过均值的X
C = np.mat(np.dot(X_new.T, X_new))
value, vector = np.linalg.eigh(C)
# 注意np.linalg.eigh是针对对称矩阵的,np.linalg.eig是针对非对称矩阵的
# value是特征值,vector是特征向量,注意vector的每一列是一个特征向量
V = []
# 这里存特征向量
for each in range(40):
# 我们这里保留40个特征向量
index = value.tolist().index(max(value.tolist()))
V.append(vector[index])
value[index] = -float("inf")
V = np.squeeze(np.array(V)).swapaxes(0, 1)
Cx = np.dot(X_new, V)
# 这个矩阵就是P
vector_list = np.split(Cx, 40, axis=1)
for i, each in enumerate(vector_list):
face = np.concatenate(np.split(each, 200, axis=0), axis=1).astype(np.uint8)
while True:
cv2.imshow("face", face)
if cv2.waitKey(1000 // 12) & 0xff == ord("q"):
break
cv2.destroyAllWindows()
cv2.imwrite("/home/zhangchen/face_{}.png".format(i), face)
# 补充部分
database = []
# 存放数据集中每个人脸的线性变换向量
for each in np.split(X, 134, axis=1):
database.append(np.dot(Cx.swapaxes(0, 1), each))
test_img = cv2.imread("/home/zhangchen/yalefacestest/10_noglasses.png", 0)
# 加载一张测试图
test_list = np.split(test_img, 200, axis=1)
test = np.dot(Cx.swapaxes(0, 1), np.concatenate(test_list, axis=0))
norm = []
# 存放计算欧式距离的结果
# 然后计算test与database中每一个元素的欧式距离,最小的就是对应同一个人
for each in database:
norm.append(np.linalg.norm(each - test))
which_person = cv2.imread(dataset_list[norm.index(min(norm))], 0)
while True:
cv2.imshow("face", which_person)
if cv2.waitKey(1000 // 12) & 0xff == ord("q"):
break
cv2.destroyAllWindows()
# 补充部分
face_average(dataset_list)
待识别人脸:
识别到人脸:
识别正确!
待识别人脸:
识别到人脸:
识别错误!
待识别人脸:
识别到人脸:
识别正确!
可见效果还是不错的。
结语
如果您有修改意见或问题,欢迎留言或者通过邮箱和我联系。
手打很辛苦,如果我的文章对您有帮助,转载请注明出处。