PCA的定义
主成分分析，即Principal Component Analysis（PCA），是多元统计中的重要内容，也广泛应用于机器学习和其它领域。它的主要作用是对高维数据进行降维。PCA把原先的n个特征用数目更少的k个特征取代，新特征是旧特征的线性组合，这些线性组合最大化样本方差，尽量使新的k个特征互不相关。
PCA的实现步骤
计算样本每个特征的平均值；
每个样本数据减去该特征的平均值；
对样本组成的矩阵求协方差矩阵；
对协方差矩阵做分解得到特征值和特征向量，并将特征值和特征向量重新按照降序排列；
对特征值求取累计贡献率；
通过累计贡献率选取适合的特征向量子集合；
对原始数据进行转换
注意：假设有m个样本，每个样本有n个特征，即组成m行n列的矩阵A,如果是图像的话，由于图像是二维向量，所以先要将图像转换为一个行向量。求均值是对A矩阵的每一列求均值，此处协方差矩阵应该是A.T * A而不是A*A.T,最后得到的协方差矩阵的形状应为n*n。
协方差矩阵分解
原始算法（硬算，特征数很大时计算量非常大）
numpy的svd算法来分解（使用svd分解出来的右奇异矩阵和直接求出来的特征向量矩阵是一个意思，但是svd并不是通过协方差矩阵来计算的，而是通过某种强大的算法近似求解）
scikit-learn实现的PCA
PCA实现
mat = [[-1,-1,0,2,1],[2,0,0,-1,-1],[2,0,1,1,0]]
mat = np.array(mat)
#原始算法
def origin_algorithm_pca(mat, 2):
#对列求均值
mean = np.mean(mat, 0)
#去中心化
A = mat - mean
#求方差矩阵
cov = np.dot(A.T, A)/2
#求解协方差矩阵的特征值和特征向量
s, vector = np.linalg.eig(cov)
#对特征进行压缩
compression = np.dot(A, vector[:2,:].T)
return compression
def numpy_svd(mat, 2):
#对列求均值
mean = np.mean(mat, 0)
#去中心化
A = mat - mean
#
u,s,VT = np.linalg.svd(A)
#压缩
compression = np.dot(A, VT[:2,:])
return compressionPCA实现人脸识别
将每张图像拉直为一维行向量
将所有样本存储为矩阵形式，矩阵的每一列表示图像的一个特征也是一个像素
按列求取均值，并用原始矩阵减去均值
利用A.T * A得到协方差矩阵(这一步代码中直接利用svd方法得到U,S,VT)
将原始矩阵乘以VT得到压缩后的矩阵，将VT存储起来用于测试测试样本
train_A*VT和test_A*VT，将test_A*VT中的每个样本与train_A*VT计算平方差，并按行求和后排序，得到距离最小的对应的train_A的label，就是预测结果
根据预测结果和真实标签计算准确率
# coding:utf-8
import os
from numpy import *
import numpy as np
import cv2
import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 图片矢量化
def img2vector(image):
img = cv2.imread(image, 0) # 读取图片
rows, cols = img.shape #获取图片的像素
imgVector = np.zeros((1, rows * cols))
imgVector = np.reshape(img, (1, rows * cols))#使用imgVector变量作为一个向量存储图片矢量化信息，初始值均设置为0
return imgVector
orlpath = "./ORL"
# 读入人脸库,每个人随机选择k张作为训练集,其余构成测试集
def load_orl(k):#参数K代表选择K张图片作为训练图片使用
'''
对训练数据集进行数组初始化，用0填充，每张图片尺寸都定为112*92,
现在共有40个人，每个人都选择k张，则整个训练集大小为40*k,112*92
'''
train_face = np.zeros((40 * k, 112 * 92))
train_label = np.zeros(40 * k) # [0,0,.....0](共40*k个0)
test_face = np.zeros((40 * (10 - k), 112 * 92))
test_label = np.zeros(40 * (10 - k))
# sample=random.sample(range(10),k)#每个人都有的10张照片中，随机选取k张作为训练样本(10个里面随机选取K个成为一个列表)
sample = random.permutation(10) + 1 # 随机排序1-10 (0-9）+1
for i in range(40): # 共有40个人
people_num = i + 1
for j in range(10): # 每个人都有10张照片
image = orlpath + '/s' + str(people_num) + '/' + str(sample[j]) + '.jpg'
# 读取图片并进行矢量化
img = img2vector(image)
if j < k:
# 构成训练集
train_face[i * k + j, :] = img
train_label[i * k + j] = people_num
else:
# 构成测试集
test_face[i * (10 - k) + (j - k), :] = img
test_label[i * (10 - k) + (j - k)] = people_num
return train_face, train_label, test_face, test_label
# 定义PCA算法
def PCA(data, r):#降低到r维
data = np.float32(np.mat(data))
rows, cols = np.shape(data)
# print(rows, cols)
data_mean = np.mean(data, 0) # 对列求平均值
A = data - np.tile(data_mean, (rows, 1)) # 将所有样例减去对应均值得到A
u, s, VT = np.linalg.svd(A) #利用svd求解右奇异向量即需要将原始数组映射的向量空间矩阵
V_r = VT[:, 0:r] # 按列取前r个特征向量
#将原始数据乘上新的空间向量得到降维后的矩阵
final_data = A * V_r
return final_data, data_mean, V_r
def face_recongize():
for r in range(10, 81, 10): # 最多降到40维,即选取前40个主成分（因为当k=1时，只有40维)
print("当降维到%d时" % (r))
x_value = []
y_value = []
train_face, train_label, test_face, test_label = load_orl(7) # 得到数据集
# 利用PCA算法进行训练
data_train_new, data_mean, V_r = PCA(train_face, r)
# print(data_train_new.shape)
num_train = data_train_new.shape[0] # 训练脸总数
num_test = test_face.shape[0] # 测试脸总数
temp_face = test_face - np.tile(data_mean, (num_test, 1))
data_test_new = temp_face * V_r # 得到测试脸在特征向量下的数据
data_test_new = np.array(data_test_new) # mat change to array
# print(data_test_new.shape)
data_train_new = np.array(data_train_new)
# 测试准确度
true_num = 0
for i in range(num_test):
testFace = data_test_new[i, :]
# print(testFace.shape)
diffMat = data_train_new - np.tile(testFace, (num_train, 1)) # 训练数据与测试脸之间距离
print(diffMat.shape)
sqDiffMat = diffMat ** 2
sqDistances = sqDiffMat.sum(axis=1) # 按行求和
sortedDistIndicies = sqDistances.argsort() # 对向量从小到大排序，使用的是索引值,得到一个向量
indexMin = sortedDistIndicies[0] # 距离最近的索引
if train_label[indexMin] == test_label[i]:
true_num += 1
else:
pass
accuracy = float(true_num) / num_test
x_value.append(7)
y_value.append(round(accuracy, 2))
print('当每个人选择%d张照片进行训练时，The classify accuracy is: %.2f%%' % (7, accuracy * 100))
if __name__ == '__main__':
face_recongize()PCA的优缺点
优点
仅仅需要以方差衡量信息量，不受数据集以外的因素影响。
各主成分之间正交，可消除原始数据成分间的相互影响的因素
计算方法简单，主要运算是特征值分解，易于实现。
缺点
主成分各个特征维度的含义具有一定的模糊性，不如原始样本特征的解释性强。
方差小的非主成分也可能含有对样本差异的重要信息，因降维丢弃可能对后续数据处理有影响。
SVD原理
任意形状的矩阵A经过svd分解会得到U，Σ，VT（V的转置）三个矩阵，其中U是一个M×M的方阵，被称为左奇异向量，方阵里面的向量是正交的；Σ是一个M×N的对角矩阵，除了对角线的元素其他都是0，对角线上的值称为奇异值；VT(V的转置)是一个N×N的矩阵，被称为右奇异向量，方阵里面的向量也都是正交的。
可以计算ΣΣ前x个奇异值的平方和占所有奇异值的平方和的比例，如果大于90%，我们就选这x个奇异值重构矩阵（剩余的数据代表的可能是噪声，无用数据）。
那么重构的A'为：
SVD的应用
svd实现图像压缩
import scipy.misc
from scipy.linalg import svd
import matplotlib.pyplot as plt
import numpy as np
import numpy
img = scipy.misc.face()[:,:,0]
img = np.array(img)
U,s,Vh=svd(img)
#将特征压缩到十维，并还原图像
A = numpy.dot(U[:,0:10],numpy.dot(numpy.diag(s[0:10]),Vh[0:10,:]))
plt.subplot(111,aspect='equal')
plt.title(":10")
plt.imshow(A)
plt.imsave('a10.png', A)svd实现信号去噪
import matplotlib.pylab as plt
import numpy as np
import numpy
from scipy.linalg import svd
x = np.linspace(-np.pi, np.pi, 400)
mu, sigma = 0, 0.05
ss = np.random.normal(mu, sigma, 400)
a = (np.sin(x) + ss).reshape(20, 20)
U,s,Vh = svd(a)
print s
A = numpy.dot(U[:,0:2],numpy.dot(numpy.diag(s[0:2]),Vh[0:2,:]))
aa = A.reshape((400,))
plt.subplot(221,aspect='equal')
plt.plot(x, np.sin(x))
plt.title("sin")
plt.subplot(222,aspect='equal')
plt.plot(x, np.sin(x) + ss)
plt.title("sin with noise")
plt.subplot(223,aspect='equal')
plt.title("sin remove noise x = 2")
plt.plot(x, aa)
A = numpy.dot(U[:,0:1],numpy.dot(numpy.diag(s[0:1]),Vh[0:1,:]))
aa = A.reshape((400,))
plt.subplot(224,aspect='equal')
plt.title("sin remove noise x = 1")
plt.plot(x, aa)
plt.show()第一副图是纯sin波形，第2副图是加了噪声的波形，第三副选择2个奇异值重构A结果很接近标准sin波形实现了去噪，第四副图是选了1个奇异值重构A的结果。
参考文献
http://liao.cpython.org/scipy06.html
https://www.cnblogs.com/jclian91/p/8024101.html
https://blog.csdn.net/m0_50572604/article/details/121018988
https://www.cnblogs.com/pinard/p/6239403.html
https://www.cnblogs.com/pinard/p/6251584.html
https://github.com/Gaoshiguo/PCA-Principal-Components-Analysis
