纯手写，有助于深入理解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