高维数据样本集的协方差矩阵的求解及算法实现

来源：互联网发布：施里芬计划知乎编辑：程序博客网时间：2024/04/27 16:45

步骤：

1、获取样本集矩阵（本文中行表示样本，列表示维度）

void PCA::orMatrix(int _sapNum, int _dimt) { //初始样本集矩阵sapNum_ = _sapNum;dimt_ = _dimt;std::cout << "维数为: " << dimt_ << "样本个数为: " << sapNum_ << '\n';orAry = (double **)malloc(sapNum_*sizeof(double));for(int i=1; i<=sapNum_; i++)*(orAry+i) = (double *)malloc(dimt_*sizeof(double));std::cout << "input the original vectors" << '\n';for(i=1; i<=sapNum_; i++)for(int j=1; j<=dimt_; j++)std::cin >> *(*(orAry+i)+j);}

输入数据如下

即为 10*3 阶矩阵，其中维数为3（即3个变量），索取样本数为10（即10个样本向量）

2、计算每一维的均值（即每一列的均值）

void PCA::meanVal() { //各维的均值meanValAry = (double *)malloc(dimt_*sizeof(double));double sum=0.0;for(int i=1; i<=dimt_; i++) {sum = 0.0;for(int j=1; j<=sapNum_; j++)sum += *(*(orAry+j)+i);*(meanValAry+i) = sum/sapNum_;}}

结果为：第一列 16.1 ；第二列 24.2 ；第三列 19.8

3、计算出每一维减去相应均值后的新样本集矩阵

void PCA::meanMatrix() { //减去均值后的样本集矩阵meanAry = (double **)malloc(sapNum_*sizeof(double));for(int i=1; i<=sapNum_; i++)*(meanAry+i) = (double *)malloc(dimt_*sizeof(double));for(i=1; i<=dimt_; i++)for(int j=1; j<=sapNum_; j++)*(*(meanAry+j)+i) = *(*(orAry+j)+i)- *(meanValAry+i);}

方法一：

4、按公式 =（新样本集矩阵的第 i 列的转置 * 新样本的第 j 列）/（样本数（10）-1）计算

for(i=1; i<=dimt_; i++)for(int j=1; j<=dimt_; j++) {double data = 0.0;for(int k=1; k<=sapNum_; k++) data += *(*(meanAry+k)+i)* (*(*(meanAry+k)+j));data /= sapNum_-1;*(*(covAry+i)+j) = data;

最后求得协方差矩阵（3*3 阶方阵）

方法二：

4、求新样本集的转置矩阵

void PCA::meanMatrixTrs() { //减去均值后的样本集矩阵的转置矩阵meanAryTrs = (double **)malloc(dimt_*sizeof(double));for(int i=1; i<=dimt_; i++)*(meanAryTrs+i) = (double *)malloc(sapNum_*sizeof(double));for(i=1; i<=dimt_; i++)for(int j=1; j<=sapNum_; j++)*(*(meanAryTrs+i)+j) = *(*(meanAry+j)+i);}

5、（新样本集的转置矩阵*新样本集矩阵本身）/（样本数-1）即得协方差矩阵

for(i=1; i<=dimt_; i++) {for(int j=1; j<=dimt_; j++) {double data = 0.0;for(int k=1; k<=sapNum_; k++)data += *(*(meanAryTrs+i)+k)* (*(*(meanAry+k)+j));data /= sapNum_-1;*(*(covAry+i)+j) = data;}}

结果与方法一相同