主成分分析

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

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

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

概率论基础

线性代数基础

本文主要内容如下:

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

原理解析

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

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

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

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

假设有MN列向量数据$X\in R^{N\times M}$,第$i$个字段数据的方差计算如下:

$\mu_{i}$表示第$i$个字段数据的均值

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

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

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

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

原始数据$X$有$N$维,表示有$N$个字段,计算字段之间的协方差矩阵$C\in R^{N\times N}$,对角元素表示方差,非对角元素表示协方差

其中$c_{ii}$表示第$i$个字段的随机变量,$c_{ij}$表示第$i$个和第$j$个字段的随机变量

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

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

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

其中$\lambda_{1}, \lambda_{2}, …, \lambda_{n}$为$C$的特征值

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

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

  • 如何执行降维操作?

前$k$个特征向量(转置成行向量)组成的矩阵$P\in R^{k\times N}$,所以

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

需不需要数据去量纲?

参考:

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

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

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

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

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

如何重建原始数据

对投影数据进行逆操作

需要每个字段加上均值

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

如何确定$k$值

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

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

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

操作步骤

假设有$M$条$N$维列向量数据

  1. 将原始数据按列排成$N$行$M$列矩阵$X\in R^{N\times M}$
  2. 对$X$进行标准化操作(每行减去均值后再除以标准差
  3. 求协方差矩阵$C=\frac {1}{M}XX^{T}$
  4. 求协方差矩阵特征值及对应的特征向量
  5. 将特征向量(行向量形式)按对应特征值大小从上到下排列成矩阵,取前$k$行组成矩阵$P\in R^{k\times N}$
  6. $Y=PX = R^{k\times N}\cdot R^{N\times M}=R^{k\times 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\times N$的实矩阵或复矩阵
  • full_matrices默认为True,表示返回值$u$和$v$大小为$M\times M$和$N\times N$;否则形状为$M\times K$和$K\times 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作为预处理操作,因为在无监督降维过程中,一些关键结构信息会因为方差小而被处理掉

坚持原创技术分享,您的支持将鼓励我继续创作!