【深圳大学机器学习】实验一：PCA算法
实验目的
1、实现PCA算法的人脸重构，即用20,40,60,80,...,160个投影来重构图像的效果。
2、实现PCA算法的人脸识别，给出不少于三个人脸数据库上 10,20,30,...,160维的人脸识别识别率，并打印出识别率变化曲线图。
3、用PCA用来进行人脸图像降维，实现3个不同数据集多个子集（或多个类）的二维和三维空间实现数据的可视化，不同类用不同颜色表示。
4、比较PCA、2DPCA、B2PCA在各个人脸数据集上的最高识别率，比较重构图像的差异。
实验步骤：
一，基于PCA算法的人脸重构与人脸识别，以及数据可视化
1，（条件：基于AR数据集的前120人，每人取前三张图像作为训练数据，后七张作为识别数据）
原图
不同维度下的重构效果（AR001-1.tif）
| 维度 | 正确率 |
| 20 | 62.70% |
| 40 | 68.52% |
| 60 | 71.33% |
| 80 | 74.27% |
| 100 | 74.53% |
| 120 | 75.20% |
| 140 | 75.41% |
| 160 | 75.89% |
取用不同数量特征值和正确率 基于此情况下人脸识别识别率曲线
数据可视化
二维空间 三维空间
2，（条件：基于orl数据集的前40人，每人取前三张图像作为训练数据，后七张作为识别数据）
原图
不同维度下重构效果（orl1_1.bmp）
| 维度 | 正确率 |
| 20 | 78.92% |
| 40 | 83.57% |
| 60 | 85.00% |
| 80 | 85.71% |
| 100 | 86.42% |
| 120 | 86.07% |
| 140 | 86.07% |
| 160 | 86.07% |
取用不同数量特征值和正确率 基于此情况下人脸识别识别率曲线
数据可视化
二维空间 三维空间
3，（条件：基于Yale数据集的前15人，每人取前三张图像作为训练数据，后七张作为识别数据）
原图
不同维度下重构效果(subject01_1.bmp)
| 维度 | 正确率 |
| 20 | 87.73% |
| 40 | 86.77% |
| 60 | 86.77% |
| 80 | 86.77% |
| 100 | 86.77% |
| 120 | 86.77% |
| 140 | 86.77% |
| 160 | 86.77% |
取用不同数量特征值和正确率 基于此情况下人脸识别识别率曲线
数据可视化
二维空间 三维空间
二，PCA,2DPCA,B2DPCA重构图像与识别率差异
条件：在ORL数据集上的重构效果（40个人，取前六张训练，后四张识别的情况下展示重构效果）
重构效果：
原图
PCA：
2DPCA:
B2DPCA
识别正确率（orl数据集，40个人，前六张训练，后四张识别）：
PCA:
最高正确率为96.3%
2DPCA:
最高正确率为97.5%
B2DPCA:
最高正确率为97.5%
尝试不同训练图像与识别图像比例后结果如下：
| 训练情况 | 训练3，识别7 | 训练4，识别6 | 训练5，识别5 | 训练6，识别4 | 训练7，识别3 |
| PCA | 86.4% | 89.6% | 90.0% | 96.3% | 96.7% |
| 2DPCA | 93.2% | 93.3% | 93.5% | 97.5% | 97.5% |
| B2DPCA | 93.2% | 93.3% | 93.5% | 97.5% | 97.5% |
从表中可以知道，2DPCA和B2DPCA在ORL训练集上的最高正确率基本相同，高于PCA的正确率，同时随着训练数据的增加，PCA与它们识别正确率的差距逐渐减小。
进一步发现，在训练图像6张，识别图像四张的情况下，PCA，2DPCA,B2DPCA识别时间如下：
| PCA | 2DPCA | B2DPCA |
| 4.09s | 1.59s | 1.22s |
可以看出，2DPCA和B2DPCA虽然在正确率上差别不大，但B2DPCA却比2DPCA更加节省时间，二者所耗时间均远远小于传统PCA算法。
代码：
PCA:
%1，测试与训练图像的录入（附加计时以查看效率）
tic
reshaped_faces=[];
%训练图像读取（以orl为例）
for i=1:40for j=1:10 a=imread(strcat('C:\Users\18023\Desktop\粤海小清华\homework\计导文件\人脸数据库\ORL56_46\orl',num2str(i),'_',num2str(j),'.bmp')); b = reshape(a,2576,1); %将每一张人脸拉成列向量b=double(b); %uint-8转换为double类型，避免人脸展示时全灰的影响 reshaped_faces=[reshaped_faces, b]; %储存每一列转化后的列向量到矩阵reshaped_faces中end
end%用每个人前六张作为训练图，后四张作为测试图
%初始化测试集集矩阵
test_data_index = [];
%初始化训练集矩阵
train_data_index = [];
for k=0:(size(reshaped_faces,2)/j)-1train_data_index= [train_data_index j*k+1:j*k+3];%建立一个用于查找所需训练图像的目录test_data_index= [test_data_index j*k+4:j*(k+1)];%同理，是测试图像的目录
end
train_data = reshaped_faces(:,train_data_index);%根据index读取矩阵 reshaped_faces的数据添加到矩阵train_data中
test_data = reshaped_faces(:, test_data_index);%根据index1读取矩阵 reshaped_faces的数据添加到矩阵test_data中% 2-特征向量矩阵（all_eigen_face)求得过程
% 求平均脸
mean_face = mean(train_data,2);%用mean函数求训练集矩阵的平均值,2是对每一行求平均，最后得到一个和原先尺寸一样的一个人脸，每个像素是平均颜色
% 中心化
normalizeddata= (train_data - mean_face);%得到中心化人脸（此处减去一列等于traindata每一列都减去这一列）
%PCA特征值分解的矩阵为去中心化的矩阵乘其转置矩阵
cov= normalizeddata * normalizeddata';
%特征值分解
[eigen_vectors, eigenvalues] = eig(cov);
% 从对角矩阵获取特征值
eigen_values = diag(eigenvalues);
% 对特征值按索引进行从大到小排序
[sorted_eigen_values, index] = sort(eigen_values, 'descend'); %此处建立index获取特征值的位置，以便寻找对应的特征向量
%按照特征值的大小排列特征向量
sorted_eigen_vectors = eigen_vectors(:, index);
all_eigen_faces = sorted_eigen_vectors;%2重构过程
single_face = normalizeddata(:,1);%代码中取第一张脸进行重构
index = 1;
figure;
for dim=[20,40,60,80,100,120,140,160]eigen_faces = all_eigen_faces(:,1:dim);%降维与再次重构rebuild_face = eigen_faces * (eigen_faces' * single_face) + mean_face;subplot(2,4, index); index = index + 1;fig = imshow(mat2gray(reshape(rebuild_face, [56, 46])));title(sprintf("维度=%d", dim)); if (dim == 160)waitfor(fig);end
end% 3.人脸识别
index = 1;
k = length(10:10:160); %选取适当间隔
Y = zeros(1, k);
% 遍历每个维度
for i = 1:kd = 10 * i; % 通过i确定当前维度% 取出相应数量特征向量eigen_faces = all_eigen_faces(:, 1:d);% 测试、训练数据降维projected_train_data = eigen_faces' * (train_data - mean_face);projected_test_data = eigen_faces' * (test_data - mean_face);% 测试前准备过程correct_predict_number = 0;test_face_number = size(projected_test_data, 2);% 读取测试数据与初始化for each_test_face_index = 1:test_face_numbereach_test_face = projected_test_data(:, each_test_face_index);% 初始化最小距离min_distance = inf;predicted_label = -1;%在欧氏距离下与每一个降维后的数据进行比较for each_train_face_index = 1:size(projected_train_data, 2)% 计算距离distance = norm(each_test_face - projected_train_data(:, each_train_face_index));% 如果当前距离小于之前的最小距离，则更新最小距离和标签（1NN，取最近的一个）if distance < min_distancemin_distance = distance;predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;endend% 获取该数据真实的标签real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;% 判断预测结果是否正确if predicted_label == real_labelcorrect_predict_number = correct_predict_number + 1; % 如果预测正确，增加正确数endend % 计算正确率correct_rate = correct_predict_number / test_face_number;% 将当前维度对应的正确率存入YY(i) = correct_rate; % 使用索引i来存储正确率，确保Y与d一致% 输出当前实验结果fprintf("维度=%d，总测试样本：%d，正确数:%d，正确率：%1f\n", d, test_face_number, correct_predict_number, correct_rate);
end
toc% 4，制作折线图使正确率随维度变化可视化
figure;
plot(10:10:160, Y, '-o'); % x轴为维度数，y轴为正确率
xlabel('维度数');
ylabel('正确率');
title('人脸识别正确率与维度变化');
waitfor(fig);%5，散点图绘制
for t=[2 3]eigen_faces=all_eigen_faces(:,1:t);projected_test_data=eigen_faces'*(test_data - mean_face);color=[];for x=1:size(test_data,2)color=[color floor((x-1)/19*2)];endif (t == 2)waitfor(scatter(projected_test_data(1, :), projected_test_data(2, :), [], color));elsewaitfor(scatter3(projected_test_data(1, :), projected_test_data(2, :), projected_test_data(3, :), [], color));endend2DPCA
%1，测试与训练图像的录入
tic
reshaped_faces=[];
%训练图像读取（以orl为例）
for i=1:40for j=1:10%a=imread(strcat('C:\Users\18023\Desktop\粤海小清华\homework\计导文件\人脸数据库\ORL56_46\orl',num2str(i),'_',num2str(j),'.bmp')); b = reshape(a,2576,1); b=double(b); reshaped_faces=[reshaped_faces, b]; end
end%用每个人前六张作为训练图，后四张作为测试图
%定义测试集集矩阵
test_data_index = [];
%定义训练集矩阵
train_data_index = [];
for k=0:(size(reshaped_faces,2)/j)-1train_data_index= [train_data_index j*k+1:j*k+6];%建立一个用于查找所需训练图像的目录test_data_index= [test_data_index j*k+7:j*(k+1)];%同理，是测试图像的目录
end
train_data = reshaped_faces(:,train_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵train_data中
test_data = reshaped_faces(:, test_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵test_data中
mean_face = mean(train_data,2);%2，2DPCA矩阵求取
%对于训练矩阵重新还原成二维叠加的状态
train_matrices = cell(size(train_data, 2), 1);
for i = 1:size(train_data, 2)train_matrices{i} = reshape(train_data(:, i), [56, 46]);
end
%所有训练图像矩阵的均值
mean_face_matrix = zeros(56, 46);
for i = 1:length(train_matrices)mean_face_matrix = mean_face_matrix + train_matrices{i};
end
mean_face_matrix = mean_face_matrix / length(train_matrices);% 初始化2DPCA协方差矩阵
cov_2dpca = zeros(46, 46);
% 计算2DPCA的协方差矩阵
for i = 1:length(train_matrices)centered_matrix = train_matrices{i} - mean_face_matrix;%每一个二维数据进行去中心化cov_contribution = centered_matrix' * centered_matrix;%求取其协方差矩阵cov_2dpca = cov_2dpca + cov_contribution;%每一个二维数据的协方差矩阵，都加到最后求解的矩阵上
end
cov_2dpca = cov_2dpca / length(train_matrices);%除以所有参与的二维矩阵数，求得均值
%cov_2dpca即为最终进行协方差分解的矩阵
%特征值分解与按照特征值从大到小排列
[eigen_vectors, eigenvalues] = eig(cov_2dpca);
eigen_values = diag(eigenvalues);
[sorted_eigen_values, index] = sort(eigen_values, 'descend');
sorted_eigen_vectors = eigen_vectors(:, index);
all_eigen_faces = sorted_eigen_vectors;figure
% 3，计算重构过程，展示第一张训练图像的重构
for idx = 1:4dim = idx * 10;eigen_faces = all_eigen_faces(:, 1:dim); % 选择前dim个特征向量% 计算投影系数test_image = train_matrices{1}; % 获取第一张训练图像centered_image = test_image - mean_face_matrix; % 中心化图像projection = centered_image * eigen_faces; % 投影到特征空间% 重构图像rebuild_face = mean_face_matrix + (projection * eigen_faces'); % 重构% 显示重构结果subplot(1, 4, idx);imshow(mat2gray(rebuild_face)); % 显示重构图像title(sprintf('维度=%d', dim));
end% 4.人脸识别
index = 1;
k = length(10:10:40);
Y = zeros(1, k);
% 将测试图像转换为矩阵形式
test_matrices = cell(size(test_data, 2), 1);
for i = 1:size(test_data, 2)test_matrices{i} = reshape(test_data(:, i), [56, 46]);
end
% 遍历每个维度
for i = 1:kd = 10* i;% 取出相应数量特征向量projection_vectors = all_eigen_faces(:, 1:d);% 计算训练和测试图像的特征train_features = zeros(size(train_data, 2), 56*d);test_features = zeros(size(test_data, 2), 56*d);% 对每张训练图像进行投影for j = 1:length(train_matrices)% 计算投影projected = train_matrices{j} * projection_vectors; % 56×d% 将投影结果展平为特征向量train_features(j, :) = projected(:)';end% 对每张测试图像进行投影for j = 1:length(test_matrices)projected = test_matrices{j} * projection_vectors;test_features(j, :) = projected(:)';end% 测试前准备过程correct_predict_number = 0;test_face_number = size(test_features, 1);% 遍历每一个待测试人脸for each_test_face_index = 1:test_face_numbereach_test_face = test_features(each_test_face_index, :)'; % 读取测试用的脸特征% 初始化最小距离min_distance = inf;predicted_label = -1;% 对比每一张训练脸for each_train_face_index = 1:size(train_features, 1)% 计算距离distance = norm(each_test_face - train_features(each_train_face_index, :)');% 如果当前距离小于之前的最小距离，则更新最小距离和标签if distance < min_distancemin_distance = distance;predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;endend% 获取实际标签real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;% 判断预测结果是否正确if predicted_label == real_labelcorrect_predict_number = correct_predict_number + 1;endend% 计算正确率correct_rate = correct_predict_number / test_face_number;Y(i) = correct_rate;fprintf("维度=%d，总测试样本：%d，正确数:%d，正确率：%1f\n", d, test_face_number, correct_predict_number, correct_rate);
end
toc% 制作折线图使正确率随维度变化可视化
figure;
plot(10:10:40, Y, '-o'); % x轴为维度数，y轴为正确率
xlabel('维度数');
ylabel('正确率');
title('人脸识别正确率与维度变化');B2DPCA
%1，测试与训练图像的录入
tic
reshaped_faces=[];
for i=1:40for j=1:10 a=imread(strcat('C:\Users\18023\Desktop\粤海小清华\homework\计导文件\人脸数据库\ORL56_46\orl',num2str(i),'_',num2str(j),'.bmp')); b = reshape(a,2576,1); b=double(b); reshaped_faces=[reshaped_faces, b]; end
end%用每个人前六张作为训练图，后四张作为测试图
%定义测试集集矩阵
test_data_index = [];
%定义训练集矩阵
train_data_index = [];
for k=0:(size(reshaped_faces,2)/j)-1train_data_index= [train_data_index j*k+1:j*k+6];%建立一个用于查找所需训练图像的目录test_data_index= [test_data_index j*k+7:j*(k+1)];%同理，是测试图像的目录
end
train_data = reshaped_faces(:,train_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵train_data中
test_data = reshaped_faces(:, test_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵test_data中
mean_face = mean(train_data,2);%2，协方差矩阵求取
%对于训练矩阵重新还原成二维叠加的状态
train_matrices = cell(size(train_data, 2), 1);
for i = 1:size(train_data, 2)train_matrices{i} = reshape(train_data(:, i), [56, 46]);
end%求所有训练图像矩阵的均值
mean_face_matrix = zeros(56, 46);
for i = 1:length(train_matrices)mean_face_matrix = mean_face_matrix + train_matrices{i};
end
mean_face_matrix = mean_face_matrix / length(train_matrices);% 行和列方向特征值分解
%行方向
cov_2dpca_col = zeros(46, 46);
for i = 1:length(train_matrices)centered_matrix = train_matrices{i} - mean_face_matrix;%求得每一个二维数据去中心化的结果cov_contribution = centered_matrix' * centered_matrix;%行方向协方差矩阵求解，用转置乘以原矩阵cov_2dpca_col = cov_2dpca_col + cov_contribution;
end
cov_2dpca_col = cov_2dpca_col / length(train_matrices);
%列方向
cov_2dpca_row = zeros(56, 56);
for i = 1:length(train_matrices)centered_matrix = train_matrices{i} - mean_face_matrix;%求得每一个二维数据去中心化的结果cov_contribution = centered_matrix * centered_matrix';%列方向协方差矩阵求解，用原矩阵乘以其转置cov_2dpca_row = cov_2dpca_row + cov_contribution;
end
cov_2dpca_row = cov_2dpca_row / length(train_matrices);% 列方向特征向量求取
[eigen_vectors_col, eigenvalues_col] = eig(cov_2dpca_col);
eigen_values_col = diag(eigenvalues_col);
[sorted_eigen_values_col, index_col] = sort(eigen_values_col, 'descend');
sorted_eigen_vectors_col = eigen_vectors_col(:, index_col);% 行方向特征向量求取
[eigen_vectors_row, eigenvalues_row] = eig(cov_2dpca_row);
eigen_values_row = diag(eigenvalues_row);
[sorted_eigen_values_row, index_row] = sort(eigen_values_row, 'descend');
sorted_eigen_vectors_row = eigen_vectors_row(:, index_row);%3，人脸重构查看效果
figure;
for idx = 1:8dim_col = idx * 5; % 列方向降维dim_row = idx * 5; % 行方向降维% 行和列方向特征向量选取eigen_faces_col = sorted_eigen_vectors_col(:, 1:dim_col);eigen_faces_row = sorted_eigen_vectors_row(:, 1:dim_row);% 获取第一张训练图像并重构test_image = train_matrices{1};centered_image = test_image - mean_face_matrix;% B2DPCA投影和重构projection = eigen_faces_row' * centered_image * eigen_faces_col;rebuild_face = mean_face_matrix + eigen_faces_row * projection * eigen_faces_col';% 显示重构结果subplot(2, 4, idx);imshow(mat2gray(rebuild_face));title(sprintf('维度=%d×%d', dim_row, dim_col));
end%4,人脸识别部分
index = 1;
k = length(10:10:40);
Y = zeros(1, k);% 将测试图像转换为矩阵形式
test_matrices = cell(size(test_data, 2), 1);
for i = 1:size(test_data, 2)test_matrices{i} = reshape(test_data(:, i), [56, 46]);
end% 遍历每个维度
for i = 1:kd = 10 * i; % 每个方向的维度% 取出特征向量projection_vectors_col = sorted_eigen_vectors_col(:, 1:d);projection_vectors_row = sorted_eigen_vectors_row(:, 1:d);% 计算训练和测试图像的特征train_features = zeros(size(train_data, 2), d * d);test_features = zeros(size(test_data, 2), d * d);% 对每张训练图像进行双向投影for j = 1:length(train_matrices)centered_matrix = train_matrices{j} - mean_face_matrix;projected = projection_vectors_row' * centered_matrix * projection_vectors_col;train_features(j, :) = projected(:)';end% 对每张测试图像进行双向投影for j = 1:length(test_matrices)centered_matrix = test_matrices{j} - mean_face_matrix;projected = projection_vectors_row' * centered_matrix * projection_vectors_col;test_features(j, :) = projected(:)';end% 测试前准备过程correct_predict_number = 0;test_face_number = size(test_features, 1);% 遍历每一个待测试人脸for each_test_face_index = 1:test_face_numbereach_test_face = test_features(each_test_face_index, :)'; % 读取测试用的脸特征% 初始化最小距离min_distance = inf;predicted_label = -1;% 对比每一张训练脸for each_train_face_index = 1:size(train_features, 1)% 计算距离distance = norm(each_test_face - train_features(each_train_face_index, :)');% 如果当前距离小于之前的最小距离，则更新最小距离和标签if distance < min_distancemin_distance = distance;predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;endend% 获取实际标签real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;% 判断预测结果是否正确if predicted_label == real_labelcorrect_predict_number = correct_predict_number + 1;endend% 计算正确率correct_rate = correct_predict_number / test_face_number;Y(i) = correct_rate;fprintf("维度=%d，总测试样本：%d，正确数:%d，正确率：%1f\n", d, test_face_number, correct_predict_number, correct_rate);
end
toc% 制作折线图使正确率随维度变化可视化
figure;
plot(10:10:40, Y, '-o'); % x轴为维度数，y轴为正确率
xlabel('维度数');
ylabel('正确率');
title('人脸识别正确率与维度变化');