【深圳大学机器学习】实验一：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）

取用不同数量特征值和正确率        基于此情况下人脸识别识别率曲线

数据可视化

             二维空间                       三维空间

2，（条件：基于orl数据集的前40人，每人取前三张图像作为训练数据，后七张作为识别数据）

原图

不同维度下重构效果（orl1_1.bmp）

取用不同数量特征值和正确率        基于此情况下人脸识别识别率曲线          

数据可视化

二维空间                       三维空间

3，（条件：基于Yale数据集的前15人，每人取前三张图像作为训练数据，后七张作为识别数据）

原图

不同维度下重构效果(subject01_1.bmp)

取用不同数量特征值和正确率         基于此情况下人脸识别识别率曲线

数据可视化

             二维空间                       三维空间

二，PCA,2DPCA,B2DPCA重构图像与识别率差异

条件：在ORL数据集上的重构效果（40个人，取前六张训练，后四张识别的情况下展示重构效果）

重构效果：

原图

PCA：

2DPCA:

B2DPCA

识别正确率（orl数据集，40个人，前六张训练，后四张识别）：

PCA:

最高正确率为96.3%

2DPCA:

最高正确率为97.5%

B2DPCA:

最高正确率为97.5%

尝试不同训练图像与识别图像比例后结果如下：

从表中可以知道，2DPCA和B2DPCA在ORL训练集上的最高正确率基本相同，高于PCA的正确率，同时随着训练数据的增加，PCA与它们识别正确率的差距逐渐减小。

进一步发现，在训练图像6张，识别图像四张的情况下，PCA，2DPCA,B2DPCA识别时间如下：

可以看出，2DPCA和B2DPCA虽然在正确率上差别不大，但B2DPCA却比2DPCA更加节省时间，二者所耗时间均远远小于传统PCA算法。

代码：

PCA:

%1，测试与训练图像的录入（附加计时以查看效率）
tic
reshaped_faces=[];
%训练图像读取（以orl为例）
for i=1:40
    for 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)-1
     train_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:k
    d = 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_number
        each_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_distance
                min_distance = distance;
                predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;
            end
        end
        % 获取该数据真实的标签
        real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;
        % 判断预测结果是否正确
        if predicted_label == real_label
            correct_predict_number = correct_predict_number + 1; % 如果预测正确，增加正确数
        end
    end 
    % 计算正确率
    correct_rate = correct_predict_number / test_face_number;
    % 将当前维度对应的正确率存入Y
    Y(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)];
    end
     if (t == 2)
        waitfor(scatter(projected_test_data(1, :), projected_test_data(2, :), [], color));
     else
        waitfor(scatter3(projected_test_data(1, :), projected_test_data(2, :), projected_test_data(3, :), [], color));
     end
 end

2DPCA

%1，测试与训练图像的录入
tic
reshaped_faces=[];
%训练图像读取（以orl为例）
for i=1:40
    for 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)-1
    train_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:4
    dim = 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:k
    d = 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_number
        each_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_distance
                min_distance = distance;
                predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;
            end
        end
        
        % 获取实际标签
        real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;
        
        % 判断预测结果是否正确
        if predicted_label == real_label
            correct_predict_number = correct_predict_number + 1;
        end
    end
    
    % 计算正确率
    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:40
    for 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)-1
    train_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:8
    dim_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:k
    d = 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_number
        each_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_distance
                min_distance = distance;
                predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;
            end
        end
        
        % 获取实际标签
        real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;
        
        % 判断预测结果是否正确
        if predicted_label == real_label
            correct_predict_number = correct_predict_number + 1;
        end
    end
    
    % 计算正确率
    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('人脸识别正确率与维度变化');