BDPCA—双向主成分分析详解及编程

最新推荐文章于 2024-04-23 18:05:34 发布

置顶翻滚的石子

最新推荐文章于 2024-04-23 18:05:34 发布

阅读量1.1k

点赞数 1

分类专栏：特征提取文章标签：机器学习线性代数

本文链接：https://blog.csdn.net/qq_38962956/article/details/116083604

版权

特征提取专栏收录该内容

2 篇文章 0 订阅

订阅专栏

一、为什么提出BDPCA（why）

（一）PCA的过拟合问题

（二）2DPCA解决PCA的过拟合问题以及不足

2.1 2DPCA如何解决PCA的过度拟合问题

第一幅插图是对此篇文章的总结，建议看完文章以后再行观看，也希望大家阅读此篇文章之前先行阅读我的2DPCA文章《特征提取—二维主成分分析（2DPCA）详解及编程》。本文主要参考论文《Bidirectional PCA With Assembled Matrix Distance Metric for Image Recognition》。

一、为什么提出BDPCA（why）

（一）PCA的过拟合问题

PCA过拟合问题指的是训练样本集的归一化均方误差远小于测试样本集的归一化均方误差。使用 $x_{i}$ 和 $x\widetilde{}_{i}$ 分别表示第 $i$ 个样本向量及其重构样本向量， $\bar{x}$ 表示样本向量的平均向量，归一化均方误差可以由下式计算：

$MSE^{train}=\frac{\sum_{i=1}\left \| x^{train}_{i}-\widetilde{x}^{train}_{i} \right \|^{2}}{\sum_{i=1}\left \| x^{train}_{i}-\bar{x}^{train}_{i} \right \|^{2}}$

$MSE^{test}=\frac{\sum_{i=1}\left \| x^{test}_{i}-\widetilde{x}^{test}_{i} \right \|^{2}}{\sum_{i=1}\left \| x^{test}_{i}-\bar{x}^{test}_{i} \right \|^{2}}$

上述公式的分子代表的是投影误差（重构误差），PCA的优化目标是投影误差最小。注：在进行PCA的公式推导时，大部分博文使用的优化目标是投影后数据方差最大，如果不理解投影误差最小这一优化目标，可以百度查找PCA博文进行阅读理解，之所以使用这个概念，是因为更容易解释PCA的过拟合问题。如下图，将二维数据投影到直线，降维到一维数据，数据点与投影直线的距离就是投影误差。

PCA的过拟合问题可以从以下两点进行解释：

1、PCA及其系列算法属于无监督学习，易导致过拟合问题：

PCA系列算法的优化目标是使得图像协方差矩阵的迹最大，且方差越大，信息量越大。因此第一个主成分特征对应信息量最大，第二个主成分特征对应信息量排第二，依次类推。但是PCA系列算法是无监督学习算法，不考虑训练数据与训练标签的关系，其计算的特征向量是使得整个训练集协方差矩阵的迹最大，即使得训练集投影后的数据的整体方差最大，而不是训练集类别投影数据之间方差最大。因此前几个主成分特征虽然对应信息量大，但不一定就是最佳的分类特征，也许某些方差小的特征直接决定了分类结果，因此特征选择(Feature Selection)研究方向一直伴随特征提取方向同行。

举一个分类识别的极端例子，要根据上半身图像分类识别十个人，十个人长相极为相似，仅仅是眼睛颜色不同。训练集中十个人穿着不同颜色的上衣，眼镜和项链，而且同一个人的上衣、眼镜和项链颜色不同；测试集中十个人也是皆穿着不用颜色的上衣、眼镜和项链，且同一个人的上衣、眼镜和项链颜色不同，但是每个人的上衣、眼镜和项链颜色与训练集不同。首先使用PCA进行特征提取，显而易见，第一个主成分对应信息是上衣，第二个主成分对应信息是眼镜边框，第三个主成分对应信息是项链，第四个主成分对应信息是眼镜，后续的主成分对应着脸部其他特征。此时使用测试集进行识别，倘若使用前三个主成分，势必识别错误，而仅仅使用第四个主成分恰恰会识别成功。

2、PCA由于训练样本容量小容易导致过拟合，即样本维度远大于样本数量时，PCA在训练集上过拟合：

首先给出一个结论，在使用PCA进行数据降维时，特征维度要小于样本数量，接下来对此进行解释，如果有解释有误，还望指正：

（1）将高维数据降维到一维，是将数据投影到直线上；将高维数据降维到二维，是将数据投影到平面上；将高维数据降维到三维，是将数据投影到立体上；......。

（2）常讲，一条直线需要两点确定，一个平面需要三点确定，一个立体需要四点确定，......。这里的点我们理解为样本点。

（3）以二维数据降维为例，下图展示了六个样本点的降维过程，可以找到一条直线使得样本点的投影误差最小；如果只有两个样本点，可以找到过两个样本点的唯一条直线使得投影误差为0；如果只有一个样本点，可以找到过该样本点的无数条直线，使得投影误差为0。对于第三种情况，无数条直线都可以使得训练集投影误差为0，但到底那一条直线也使得测试集投影误差最小时，我们无从得知，因此降维过程是没有意义的。

（4）对于第二种情况，两个样本点确定一条投影直线，但是如果这两个训练样本点是噪声点，测试样本点距离投影直线较远时，那么测试集的投影误差会很大。在不确定测试样本点的情况下，投影直线完全受到训练样本点约束，降维过程同样也没有意义，因此我更加倾向于特征维度要小于样本量-1。对此进行推广，三个高维训练样本点无法降维到二维，以三个三维训练样本点为例，降维到二维需要计算一个投影平面使得投影误差最小，而一个平面需要三个点确定，因此三个训练样本点确定了唯一一个平面，使得投影误差为0，因此投影平面的计算完全受到了三个训练样本点约束；四个高维训练样本点无法降维到三维，一个投影立体需要四个点进行确定，因此四个训练样本点确定了唯一一个立体，使得投影误差为0，投影立体完全受到四个训练样本点约束。倘若测试集与训练集差距较大，测试集距离投影直线、平面、立体......较远，则不可避免出现：训练集投影误差为0，测试集投影误差较大，PCA过拟合到训练集。

上面对特征维度与样本数量的关系进行了解释，现在继续解释PCA由于样本容量小引起的过拟合问题。由于特征维度受到训练样本数量的约束，因此相当多的特征维度会被丢失掉，而某个丢失的特征维度也许恰恰是训练集和测试集共有的特征，这样导致了训练集投影误差小，测试集投影误差大。例如，对ORL人脸数据集的两个人的20张样本图像进行PCA特征提取，每人10张图像，一张图像尺寸为112×92。图像向量化后的维度为10304，而由于我们只有20张图像，只能降维到18维，因此10286维数据被丢失，极易导致过拟合问题。

（二）2DPCA解决PCA的过拟合问题以及不足

2.1 2DPCA如何解决PCA的过度拟合问题

再次注明一下，2DPCA解决的是由于PCA自身原因导致的过拟合问题。

首先，我们回顾一下2DPCA的原理，使用 $\left \{ X_{1},X_{2},...,X_{N} \right \}$ 表示N个二维图像，图像尺寸为 $m\times n$ ， $\overline{X}$ 表示平均图像，则图像总散射矩阵 $G_{t}$ 可由下式计算：

$G_{t}= \frac{1}{N}\sum_{i=1}^{N}(X_{i}-\overline{X})^{T}(X_{i}-\overline{X})$

将 $m\times n$ 的二维图像视作 $m$ 个 $1\times n$ 行向量，即：

$G_{t}$ 又可由下式表示：

$G_{t}= \frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{m}(x_{i}^{j}-\overline{x}^{j})^{T}(x^{j}_{i}-\overline{x}^{j})$

其中 $\bar{x}^{j}$ 表示均值图像 $\bar{X}$ 的第 $j$ 个行向量。

参考PCA，将二维图像的行向量作为训练样本，则定义行总散射矩阵 $S_{t}^{row}$ 如下：

$S^{row}_{t}= \frac{1}{Nm}\sum_{i=1}^{N}\sum_{j=1}^{m}(x_{i}^{j}-\overline{x})^{T}(x^{j}_{i}-\overline{x})$

其中 $\bar{x}$ 表示均值图像 $\bar{X}$ 的所有行向量的均值行向量。

比较上述两个公式，二者几乎一致，除了行总散射矩阵 $S_{t}^{row}$ 使用 $\bar{x}$ 代替 $\bar{x}^{j}$ ，乘以系数 $1/m$ 。因此，2DPCA是行PCA，将m×n的二维图像视作 $m$ 个 $1\times n$ 行向量，在行方向上进行投影。 $G_{t}$ 的尺寸为 $n\times n$ ，通过特征分解算法可以计算出前k个特征向量组成投影矩阵 $W$ ，尺寸为 $n\times k$ 。进一步可由 $Y=XW$ 计算出特征矩阵 $Y$ ，尺寸为 $m\times n$ ，行方向维度由 $n$ 降维到 $k$ 。

2DPCA的训练样本数为mN，样本维度为n，不再像PCA一样样本维度远大于训练样本数，不再是 $N\ll mn$ ，而是 $n\ll mN$ ，因此2DPCA解决了过拟合问题。

2.2 2DPCA的不足

2DPCA是行PCA，在行方向上进行投影，因此只是从列方向上提取特征，而忽略了行方向上的特征。

2DPCA提取的特征矩阵的维度较大。假设二维样本图像尺寸 $m\times n$ ，提取前 $k$ 个特征向量组成投影矩阵 $W$ ，则PCA计算的投影矩阵 $W_{PCA}$ 尺寸为 $mn\times k$ ，2DPCA计算的投影矩阵 $W_{2DPCA}$ 尺寸为 $n\times k$ ，将样本图像投影到投影矩阵计算特征矩阵 $Y$ ，则PCA计算的特征矩阵 $Y_{PCA}$ 尺寸为 $1\times k$ ，2DPCA计算的特征矩阵 $Y_{2DPCA}$ 尺寸为 $m\times k$ 。可见，2DPCA提取的特征维度远大于PCA。

二、BDPCA数学公式（What）

（一）模型

使用 $\left \{ X_{1},X_{2},...,X_{N} \right \}$ 表示N个二维图像，图像尺寸为 $m\times n$ 。

1、PCA将 $m\times n$ 的样本图像进行向量化为 $1\times mn$ 的向量 $x$ ，旨在使用向量 $x$ 计算前 $k$ 个 $mn$ 维标准正交列向量构成的投影矩阵 $W$ ( $mn\times k$ )作为投影轴，实现如下的特征提取：

$Y=\left ( x-\overline{x }\right )W$

其中为图像向量投影后得到的新特征，尺寸为 $1\times k$ 。

2、2DPCA旨在使用二维图像 $X$ 计算前 $k$ 个 $n$ 维标准正交列向量构成的投影矩阵 $W$ ( $n\times k$ )作为投影轴，实现如下的特征提取：

$Y=\left ( X-\overline{X }\right )W$

其中为二维图形投影后得到的新特征，尺寸为 $m\times k$ 维。

3、BDPCA旨在使用二维图像 $X$ 计算前 $k_{row}$ 个 $n$ 维标准正交列向量构成的投影矩阵 $W_{row}$ ( $m\times k_{row}$ )和前 $k_{col}$ 个 $m$ 维标准正交行向量构成的投影矩阵 $W_{col}$ ( $n\times k_{col}$ )作为投影轴，实现如下的特征提取：

$Y=W_{col}^{T}\left ( X-\overline{X }\right )W_{row}$

其中为二维图形投影后得到的新特征，尺寸为 $k_{col}\times k_{row}$ 维。

（二）策略（优化目标）

BDPCA与2DPCA相同，引入投影后特征的总散射，即投影后特征构造的协方差矩阵的迹，进行衡量投影轴或者投影后特征。通过最大化协方差矩阵的迹计算多个投影轴，构成投影矩阵。由2DPCA文章可知，最大化协方差矩阵的迹计算投影轴可以转换为计算图像协方差矩阵的特征向量。

BDPCA与2DPCA不同的是，2DPCA是行PCA，仅在行方向上进行投影，仅仅构造了行总散射矩阵进行特征分解，而BDPCA既在行方向上投影，又要在列方向上投影，因此要同时构造行总散射矩阵和列总散射矩阵。

在构造行总散射矩阵和列总散射矩阵时，整合了第二章第二小节三个公式，系数上使用了第三个公式 $\frac{1}{Nm}$ ；累加项使用了第二个公式，因为其可以转换为第一个公式，直接对二维图像进行操作。

$G_{t}= \frac{1}{N}\sum_{i=1}^{N}(X_{i}-\overline{X})^{T}(X_{i}-\overline{X})$

$G_{t}= \frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{m}(x_{i}^{j}-\overline{x}^{j})^{T}(x^{j}_{i}-\overline{x}^{j})$

$S^{row}_{t}= \frac{1}{Nm}\sum_{i=1}^{N}\sum_{j=1}^{m}(x_{i}^{j}-\overline{x})^{T}(x^{j}_{i}-\overline{x})$

行总散射矩阵 $S_{t}^{row}$ 由下式计算：

$S_{t}^{row}= \frac{1}{Nm}\sum_{i=1}^{N}(X_{i}-\overline{X})^{T}(X_{i}-\overline{X})$

列总散射矩阵 $S_{t}^{col}$ 由下式计算：

$S_{t}^{col}= \frac{1}{Nn}\sum_{i=1}^{N}(X_{i}-\overline{X})(X_{i}-\overline{X})^{T}$

（三）算法（如何计算投影向量和投影矩阵）

1、载入所有二维图像 $\left \{ X_{1},X_{2},...,X_{N} \right \}$ ，计算均值图像 $\overline{X}=\frac{1}{N}\sum_{i=1}^{N}X_{i}$ ;

2、对所有二维图像进行中心化处理，计算行总散射矩阵 $S_{t}^{row}= \frac{1}{Nm}\sum_{i=1}^{N}(X_{i}-\overline{X})^{T}(X_{i}-\overline{X})$ 和列总散射矩阵 $S_{t}^{col}= \frac{1}{Nn}\sum_{i=1}^{N}(X_{i}-\overline{X})(X_{i}-\overline{X})^{T}$ ;

3、计算行总散射矩阵的特征向量和特征值，选择前 $k_{row}$ 个最大特征值对应的特征向量组成投影矩阵 $W_{row}$ ；计算列总散射矩阵的特征向量和特征值，选择前 $k_{col}$ 个最大特征值对应的特征向量组成投影矩阵 $W_{col}$ 。

（四）特征提取和图像重构

假设图像矩阵为 $X$ ，尺寸为 $m\times n$ ；对应前 $k_{row}$ 个最大特征值对应的特征向量组成投影矩阵 $W_{row}$ ，尺寸为 $n\times k_{row}$ ；对应前 $k_{col}$ 个最大特征值对应的特征向量组成投影矩阵 $W_{col}$ ，尺寸为 $m\times k_{col}$ 。特征矩阵（特征图像） $Y$ 由下式计算：

$Y=W_{col}^{T}\left ( X-\overline{X }\right )W_{row}$

重构图像公式可由下式计算:

$\widetilde{X}=\overline{X}+W_{col}YW_{row}^{T}$

三、BDPCA编程（How）

使用VS2015+opencv 3.4.1进行编程，编程语言为C++，主要根据第二节的第三小节的步骤进行。我将ORL人脸数据库中每个人的10幅图像随机分为5幅训练图像和5幅测试图像，因此得到了训练数据集和测试数据集，大小都为200幅图像。使用训练数据集计算特征向量，载入测试数据集中一副图像进行特征提取和重构。

1、程序初始参数设置：

#include<iostream>
#include<opencv2/opencv.hpp>

using namespace std;
using namespace cv;

//void main()
//{

    //定义需要修改的变量
	int nrows = 112, ncols = 92;         //样本大小
	int rowEigenNum = 10;                 //行投影特征向量数目
	int colEigenNum = 10;                 //列投影特征向量数目
	int trainNum = 200;                 //训练样本数目
	int testNum = 200;                   //测试样本数目
	string trainPath = "F:\\机器学习\\DataSet\\facetrain\\"; //训练样本存储路径
	string testPath = "F:\\机器学习\\DataSet\\facetest\\";   //测试样本存储路径

2、载入所有二维图像 $\left \{ X_{1},X_{2},...,X_{N} \right \}$ ，计算均值图像 $\overline{X}=\frac{1}{N}\sum_{i=1}^{N}X_{i}$ ：

    //加载图像，计算均值
	Mat srcSample;
	string samplePath;
	Mat sampleAve = Mat::zeros(nrows, ncols, CV_32FC1);
	vector<Mat> trainedSamples;
	for (int i = 1; i <= trainNum; i++)
	{
		samplePath = trainPath + to_string(i) + ".png"; //图像路径
		srcSample = imread(samplePath, 0);                     //载入当前图像
		srcSample.convertTo(srcSample, CV_32FC1, 1, 0); 
		trainedSamples.push_back(srcSample.clone());           //存储当前样本
		sampleAve += srcSample/double(trainNum);               //计算均值图像
	}

3、对所有二维图像进行中心化处理，计算行总散射矩阵 $S_{t}^{row}= \frac{1}{Nm}\sum_{i=1}^{N}(X_{i}-\overline{X})^{T}(X_{i}-\overline{X})$ 和列总散射矩阵 $S_{t}^{col}= \frac{1}{Nn}\sum_{i=1}^{N}(X_{i}-\overline{X})(X_{i}-\overline{X})^{T}$ ：

    //计算行总散射矩阵和列总散射矩阵
	Mat rowScatter = Mat::zeros(ncols, ncols, CV_32FC1);
	Mat colScatter = Mat::zeros(nrows, nrows, CV_32FC1);
	for (int i = 0; i < trainedSamples.size(); i++)
	{
		srcSample = trainedSamples[i] - sampleAve;                     //样本中心化处理
		rowScatter += srcSample.t()*srcSample / double(trainNum*nrows);//计算行总散射矩阵
		colScatter += srcSample*srcSample.t() / double(trainNum*ncols);//计算列总散射矩阵
	}

4、计算行总散射矩阵的特征向量和特征值，选择前 $k_{row}$ 个最大特征值对应的特征向量组成投影矩阵 $W_{row}$ ；计算列总散射矩阵的特征向量和特征值，选择前 $k_{col}$ 个最大特征值对应的特征向量组成投影矩阵 $W_{col}$ ：

    //计算行方向投影矩阵和列方向投影矩阵
	Mat rowEigenVecs, rowEigenVals, colEigenVecs, colEigenVals;
	Mat tarRowVecs, tarColVecs;
	eigen(rowScatter, rowEigenVals, rowEigenVecs);             //计算行方向投影向量
	eigen(colScatter, colEigenVals, colEigenVecs);             //计算列方向投影向量
	tarRowVecs = (rowEigenVecs.rowRange(0, rowEigenNum)).t();  //保留前rowEigenNum个最大特征值对应的特征向量组成行方向投影矩阵
	tarColVecs = (colEigenVecs.rowRange(0, colEigenNum)).t();  //保留前colEigenNum个最大特征值对应的特征向量组成列方向投影矩阵

	//保存行方向投影矩阵和列方向投影矩阵
	CvMat storeMat;
	storeMat = tarRowVecs;
	cvSave("Data\\RowEigenVecs.txt",&storeMat);
	storeMat = tarColVecs;
	cvSave("Data\\ColEigenVecs.txt", &storeMat);

5、计算训练样本投影后的特征矩阵，用于后续分类模型训练；计算测试样本投影后的特征矩阵，用于输入到分类模型进行分类识别：

    //计算训练样本特征矩阵，按行存储，用于训练分类器模型
	Mat trainData(trainNum, rowEigenNum*colEigenNum, CV_32FC1);
	for (int i = 0; i < trainedSamples.size(); i++)
	{
		srcSample = trainedSamples[i] - sampleAve; 
		Mat tmpFeature = tarColVecs.t()*srcSample*tarRowVecs;
		Mat xi = trainData.row(i);
		tmpFeature.reshape(1, 1).convertTo(xi, CV_32FC1);
	}
	storeMat = trainData;
	cvSave("Data\\TrainFeature.txt", &storeMat);

	//计算测试样本特征矩阵，按行存储，用于后续分类识别
	Mat testData(testNum, rowEigenNum*colEigenNum, CV_32FC1);
	for (int i = 1; i <= testNum; i++)
	{
		samplePath = testPath + to_string(i) + ".png";
		srcSample = imread(samplePath, 0);
		srcSample.convertTo(srcSample, CV_32FC1);
		srcSample = srcSample - sampleAve; 
		Mat tmpFeature = tarColVecs.t()*srcSample*tarRowVecs;
		Mat xi = testData.row(i-1);
		tmpFeature.reshape(1, 1).convertTo(xi, CV_32FC1);
	}
	storeMat = testData;
	cvSave("Data\\TestFeature.txt", &storeMat);

6、载入测试集中的一副人脸图像用于特征提取和重构图像：

	//特征提取图像重构
	samplePath = testPath + "2.png";
	srcSample = imread(samplePath, 0);
	imshow("src", srcSample);
	srcSample.convertTo(srcSample, CV_32FC1);
	srcSample = srcSample - sampleAve;
	Mat feature = tarColVecs.t()*srcSample*tarRowVecs;              //计算特征矩阵
	Mat reconImg = tarColVecs*feature*tarRowVecs.t() + sampleAve;   //计算重构图像
	normalize(reconImg, reconImg, 0, 255, CV_MINMAX, CV_8UC1);      //将重构图像像素值缩放至0-255
	imshow("reconstruction", reconImg);
	samplePath = "Data\\"+to_string(rowEigenNum) + ".png";
	imwrite(samplePath, reconImg);

在程序中，我将行方向特征向量和列方向特征向量的数目设置为相同，程序结果输出：