纯手写,有助于深入理解PCA
'''
数据描述:
X 大小:m×n ---n个m维样本
每一列是一个样本
'''
import os
import glob
import cv2
import numpy as np
# 设置文件路径
os.getcwd()
os.chdir('C:/Users/h/Desktop/ML/ORL')
os.getcwd()
# 文件名称
path=[r'.\s'+str(i)+'\*.bmp' for i in range(1,41)]
# 读取图像
X = []
for j in range(len(path)+1):
images = glob.glob(path[j])
for img in images:
img = cv2.imread(img, 0)
temp = np.resize(img, (img.shape[0] * img.shape[1], 1))
X.append(temp.T)
X = np.array(X).squeeze().T
# 定义显示图像函数
def imshow(img):
cv2.imshow('imshow',img)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 分开训练集 测试集
train_data=X[:,[x for x in list(range(X.shape[1])) if x not in list(range(X.shape[1],10))]]
test_data=X[:,list(range(X.shape[1],10))]
# 求均值
mean1=np.mean(train_data,axis=1)
mean = np.expand_dims(mean1,axis=1)
# 减去均值
Z=X-mean
# 快速PCA
# 计算散布矩阵 Z.T*Z
R=np.dot(Z.T,Z)
#计算散布矩阵R的特征值、特征向量
a,b = np.linalg.eig(R)
# 计算原本特征向量
b_real = np.dot(Z,b)
# 把特征值从大到小排序的序号
a_sort=np.argsort(-a)
# 根据特征值排序结果,取前100特征值对应的特征向量
P = b_real.T[a_sort[:100]]
P = P.T
# P.T*Z得到投影a
A=np.dot(P.T,Z)
# 重建 mean+P*A
X_new =np.dot(P,A)
new = X_new + mean