主成分分析

主成分分析(princial component analysis,简称PCA)是一种无监督的数据降维操作,它通过最大化方差方法来寻找低维空间,能够有效减轻计算量的同时保证处理数据有效性

主要参考文章PCA数学原理,里面做了生动的数学原理分析

相关线性代数、几何学以及统计学数学内容参考:

概率论基础

线性代数基础

本文主要内容如下:

  1. 原理解析
  2. 如何重建原始数据
  3. 如何确定k
  4. 操作步骤
  5. numpy实现

原理解析

  • 为什么要进行数据降维?
  1. 减少内存或磁盘存储
  2. 加快算法运行
  3. 数据可视化(将多维数据将为2维或3维,实现数据可视化
  • 为什么能够进行数据降维?

原始数据之间存在结构相关性,某些数据的丢失不影响后续对数据的处理

  • 数据降维过程应该注意哪些?
  1. 保留尽可能多的原始数据内部的结构信息
  2. 处理后的数据不存在相关性
  • 如何保证能够得到尽可能多的原始数据内部的结构信息?

原始数据内部的结构信息不会在投影后重叠在一起,低维空间数据尽可能分散,也即使说投影后每个字段数据的方差越大越好

假设有MN列向量数据XRN×M,第i个字段数据的方差计算如下:

Var(X,j)=1Mj=1M(Xijμi)2

μi表示第i个字段数据的均值

μi=1Mj=1MXij

  • 如何保证处理后的数据不存在相关性

投影后不同字段之间的协方差0,表示不存在相关性,数据完全独立

i个和第l个字段数据之间的协方差计算如下:

Cov(Xi,Xl)=1Mj=1MXi,jXl,j

  • 如何计算方差和协方差?

原始数据XN维,表示有N个字段,计算字段之间的协方差矩阵CRN×N,对角元素表示方差,非对角元素表示协方差

C=[c11c12c1Nc21c22c2NcN1cN2cNN]=1MXXT

其中cii表示第i个字段的随机变量,cij表示第i个和第j个字段的随机变量

目标是协方差值为0,也就是协方差矩阵对角化

  • 如何计算协方差矩阵对角化?

协方差矩阵C是实对称矩阵,所以存在正交矩阵Q,能得到

Q1CQ=QTCQ=Λ=diag(λ1,λ2,...,λn)

其中λ1,λ2,...,λnC的特征值

正交矩阵Q由特征值对应的特征向量正交化且单位化后组成,Q的第j列特征向量对应特征值λj

所以计算协方差矩阵的特征值和特征向量,将特征值(方差)从大到小排列,取前k个对应的特征向量(列向量),就是低维空间的基

  • 如何执行降维操作?

k个特征向量(转置成行向量)组成的矩阵PRk×N,所以

Y=PXRk×M

Y表示数据X投影到矩阵P为基的低维空间数据

需不需要数据去量纲?

参考:

PCA填坑篇——使用PCA到底需不需要数据去量纲?

PCA要对数据进行预处理的原因

需要对原始数据去量纲,也就是让属性成为单纯一个数,去除属性单位带来的影响(比如cmkg

需要统一单位,所以需要对原始数据进行标准化操作,原始数据减去均值后,再除以各字段的标准差

注意:在图像操作中,每个字段取值均为[0-255],可以不进行去量纲操作,减去均值即可

如何重建原始数据

对投影数据进行逆操作

Xapprox=PTY=RN×kRk×M=RN×M

需要每个字段加上均值

(Xapprox)i=(Xapprox)i+μi

还需要对每个字段乘以字段标准差

(Xapprox)i=(Xapprox)iσi

如何确定k

降维目的是保证能够大幅降低维度,同时能够最大限度保持原始数据内部结构信息,通过投影均方误差来评价投影结果

1M||XXapprox||2

计算投影均方误差和数据集方差(已减去均值)比值,如果要保持99%方差

1M||XXapprox||21M||X||20.01

可以利用特征值来加速计算

1i=1Kλii=1Nλi0.01i=1Kλii=1Nλi0.99

操作步骤

假设有MN维列向量数据

  1. 将原始数据按列排成NM列矩阵XRN×M
  2. X进行标准化操作(每行减去均值后再除以标准差
  3. 求协方差矩阵C=1MXXT
  4. 求协方差矩阵特征值及对应的特征向量
  5. 将特征向量(行向量形式)按对应特征值大小从上到下排列成矩阵,取前k行组成矩阵PRk×N
  6. Y=PX=Rk×NRN×M=Rk×M即为降维到k维后的数据

numpy实现

特征值和特征向量计算

函数numpy.linalg.eig用于计算协方差矩阵的特征值和特征向量

numpy.linalg.eig(a)

输入参数a是方阵

返回两个数组:eigenvalue、eigenvector

1
2
3
4
5
6
7
>>> w, v = np.linalg.eig(np.diag((1,2,3)))
>>> w
array([1., 2., 3.])
>>> v # 列向量
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])

函数numpy.linalg.svd用于奇异值分解

numpy.linalg.svd(a, full_matrices=1, compute_uv=1)

  • a是大小为M×N的实矩阵或复矩阵
  • full_matrices默认为True,表示返回值uv大小为M×MN×N;否则形状为M×KK×N,其中K=min(M,N)
  • compute_uv默认为True,表示是否计算u/v

返回值是3个数组u/s/v,其中s表示特征值并且按从大到小排列,u表示其对应的特征列向量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
>>> X
array([[3, 0, 0],
[0, 1, 0],
[0, 0, 2]])
>>> u,s,v = np.linalg.svd(X)

>>> u # 列向量
array([[1., 0., 0.],
[0., 0., 1.],
[0., 1., 0.]])
>>> s
array([3., 2., 1.])
>>> v
array([[1., 0., 0.],
[0., 0., 1.],
[0., 1., 0.]])

scipy.linalg.svd同样实现了奇异值分解

注意:除以标准差时需要注意除数不能为0

pca实现

完整实现代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
# -*- coding: utf-8 -*-

# @Time : 19-6-5 上午11:36
# @Author : zj

import numpy as np
import matplotlib.pyplot as plt
import cv2
import os
import time
import warnings

warnings.filterwarnings('ignore')


def load_cifar_data(num, data_path):
"""
加载cifar数据
"""
test_dir = os.path.join(data_path, 'test', '0')
file_list = os.listdir(test_dir)

data_list = []
for i in range(num):
img_path = os.path.join(test_dir, file_list[i])
img = cv2.imread(img_path)
if img is not None:
data_list.append(img.reshape(-1))

return np.array(data_list).astype(np.uint8)


def load_mnist_data(num, data_path):
"""
加载mnist数据, 轮流提取每个类别的图像
"""
test_dir = os.path.join(data_path, 'test')

data_list = []
for i in range(num):
img_path = os.path.join(test_dir, str(i % 10), '%d.png' % (i))
img = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
if img is not None:
data_list.append(img.reshape(-1))

return np.array(data_list).astype(np.uint8)


def pca(X, ratio=0.99, **kwargs):
"""
pca降维
:param X: 大小为NxM,其中M是个数,N是维度,每个字段已是零均值
:param ratio: 表示投影均方误差和方差比值,默认为0.99,保持99%的方差
:param kwargs: 字典参数,如果指定了k值,则直接计算
:return: 降维后数据
"""
N, M = X.shape[:2]
C = X.dot(X.T) / M
u, s, v = np.linalg.svd(C)
# u, s, v = linalg.svd(C)

k = 1
if 'k' in kwargs:
k = kwargs['k']
else:
while k < N:
s_k = np.sum(s[:k])
s_N = np.sum(s)
if (s_k * 1.0 / s_N) >= ratio:
break
k += 1
p = u.transpose()[:k]
y = p.dot(X)

return y, p


def make_image(data_array, N, W, H, is_gray=False):
if is_gray:
img = np.zeros((N * H, N * W))
else:
img = np.zeros((N * H, N * W, 3))
for i in range(N):
for j in range(N):
if is_gray:
img[i * H:(i + 1) * H, j * W:(j + 1) * W] = data_array[i * N + j].reshape(H, W)
else:
img[i * H:(i + 1) * H, j * W:(j + 1) * W, :] = data_array[i * N + j].reshape(H, W, 3)
return img.astype(np.uint8)


if __name__ == '__main__':
mnist_path = '/home/zj/data/decompress_mnist/'
cifar_path = '/home/zj/data/decompress_cifar_10/'
N = 7
W = 32
H = 32
isGray = False

# data_array = load_mnist_data(N * N, mnist_path)
data_array = load_cifar_data(N * N, cifar_path)
img1 = make_image(data_array, N, W, H, is_gray=isGray)

if N == 1:
mu = 0
else:
mu = np.mean(data_array, axis=0, keepdims=True)
data_array = data_array - mu

print(data_array.shape)
start = time.time()
y, p = pca(data_array.T, ratio=0.99)
end = time.time()
print('pca need time: %f' % (end - start))
print(y.shape)
X_reduced = (p.T.dot(y)).T
X_reduced += mu

img2 = make_image(X_reduced, N, W, H, is_gray=isGray)

plt.figure(1)
if isGray:
plt.imshow(img1, cmap='gray')
else:
plt.imshow(img1)
plt.axis('off')

plt.figure(2)
if isGray:
plt.imshow(img2, cmap='gray')
else:
plt.imshow(img2)
plt.axis('off')

plt.show()

进行mnist图像降维,保留99%的方差,图像从784维降至43

降维(49, 784)(49, 43)耗时0.088249

进行cifar图像降维,保留99%的方差,图像从3072维降至49

降维(49, 3072)(49, 41)耗时9.106328秒

小结

参考:为什么PCA不被推荐用来避免过拟合?

PCA是通用的数据降维操作,但在卷积神经网络处理中不推荐使用PCA作为预处理操作,因为在无监督降维过程中,一些关键结构信息会因为方差小而被处理掉