求解协方差矩阵算法

covHead.h

#include <iostream>

class COV {
public :
	COV(void);
	void clear(void);
public :
	void rdOrMatrix(int _sapNum, int _dim); 
	//_sapNum为样本个数, _dimt为维数; 创建样本集矩阵并存入数组 orAry(_sapNum行, _dimt列)
	void meanVal(void); //计算各维度的均值, 并存入数组 meanValAry(一维数组, 维数为_dimt)
	void meanMatrix(void); //计算各维减去各维均值后的新矩阵并存入数组 orAry(_sapNum行, _dimt列)
	void meanMatrixTrs(void); //计算新矩阵的转置矩阵并存入数组 meanAryTrs(_dimt行,_sapNum列)
	void covMatrix(void); //计算协方差矩阵并存入数组 covAry(_dimt * _dimt方阵)

	void showOrMatrix(void); //输出初始样本集矩阵
	void showMeanVal(void); //输出各维度的均值
	void showMeanMatrix(void); //输出减去均值后的样本矩阵
	void showMeanMatrixTrs(void); //输出减去均值后的样本矩阵的转置矩阵
	void showCovMatrix(void); //输出协方差矩阵
private :
	int sapNum_, dim_; //样本个数;维数;降维后维数
	double **orAry; //初始样本集矩阵(减去均值后的新样本矩阵)
	double *meanValAry, **meanAryTrs; //各维度的均值;减去均值后的矩阵的转置矩阵
	double **covAry; //协方差方阵
};

 covFile.h

#include <iostream>
#include "covHead.h"
#include <malloc.h>
#include <math.h>

//////////////////////////////////////////////////////////////
COV::COV() {
	sapNum_ = 0;	dim_ = 0;
	orAry = NULL;
	meanValAry = NULL;	meanAryTrs = NULL;
	covAry = NULL;
}

//////////////////////////////////////////////////////
void COV::clear() {
	
}

/////////////////////////////////////////////////////////////////
void COV::rdOrMatrix(int _sapNum, int _dim) { //读取初始样本集矩阵
	sapNum_ = _sapNum;	dim_ = _dim;
	std::cout << "维数为: " << dim_ << "样本个数为: " << sapNum_ << '\n';
	orAry = (double **)malloc(sapNum_*sizeof(double));
	for(int i=0; i<sapNum_; i++)
		*(orAry+i) = (double *)malloc(dim_*sizeof(double)); //每一行存放一个样本向量,列数为维数,行数位总样本向量个数

	std::cout << "read the original vectors" << '\n';
	FILE *fp;
	char ch;
	if((fp=fopen("covariance.txt","rt"))==NULL) { //covariance文件中每一行代表一个样本向量,列数为维数                            
		printf("Cannot open file strike any key exit!");
		exit(0);
	}

	ch=fgetc(fp);
	std::string str=""; //这里str是类string的一个对象
	int k=0;
	while(k<sapNum_) {
		for(int i=0; i<dim_; i++) {
			while(ch==' ' || ch=='v' || ch=='\n')
				ch=fgetc(fp);
			while(ch!=' ' && ch!='\n' && ch!='v') {
				str+=ch;
				ch=fgetc(fp);
			}
			double aa=atof(str.c_str()); //将字符串转换成double型
			*(*(orAry+k)+i) = aa; //文本中的一行为矩阵valAry中的一行
			str="";
		}
		ch=fgetc(fp);
		k++;
	}
}

void COV::showOrMatrix() {
	std::cout << "============================================================" << '\n';
	std::cout << "初始样本矩阵 :" << '\n';
	for(int i=0; i<sapNum_; i++) {
		for(int j=0; j<dim_; j++)
			std::cout << *(*(orAry+i)+j) << "    ";
		std::cout << '\n';
	}			
}


/////////////////////////////////////////////////////////////////
void COV::meanVal() { //各维的均值
	meanValAry = (double *)malloc(dim_*sizeof(double));

	double sum=0.0;
	for(int i=0; i<dim_; i++) {
		sum = 0.0;
		for(int j=0; j<sapNum_; j++)
			sum += *(*(orAry+j)+i);
		*(meanValAry+i) = sum/sapNum_;
	}
}

void COV::showMeanVal() {
	std::cout << "====================================" << '\n';
	std::cout << "各维的均值 :" << '\n';
	for(int i=0; i<dim_; i++)
		std::cout << *(meanValAry+i) << "   ";
	std::cout << '\n';
}


/////////////////////////////////////////////////////////////////
void COV::meanMatrix() { //减去均值后的样本集矩阵
	for(int j=0; j<dim_; j++)
		for(int i=0; i<sapNum_; i++)
			*(*(orAry+i)+j) = *(*(orAry+i)+j)- *(meanValAry+j);
}

void COV::showMeanMatrix() {
	std::cout << "============================================================" << '\n';
	std::cout << "减去均值后的样本矩阵 :" << '\n';
	for(int i=0; i<sapNum_; i++) {
		for(int j=0; j<dim_; j++)
			std::cout << *(*(orAry+i)+j) << "    ";
		std::cout << '\n';
	}			
}


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

void COV::showMeanMatrixTrs() {
	std::cout << "============================================================" << '\n';
	std::cout << "减去均值后的样本矩阵的转置矩阵 :" << '\n';
	for(int i=0; i<dim_; i++) {
		for(int j=0; j<sapNum_; j++)
			std::cout << *(*(meanAryTrs+i)+j) << "    ";
		std::cout << '\n';
	}			
}


///////////////////////////////////////////////////////////////////
void COV::covMatrix() { //协方差矩阵
	covAry = (double **)malloc(dim_*sizeof(double));
		for(int i=0; i<dim_; i++)
			*(covAry+i) = (double *)malloc(dim_*sizeof(double));

	//计算协方差矩阵代码:减去均值的样本矩阵的转置矩阵 * 减去均值的样本矩阵本身
	for(i=0; i<dim_; i++) {
		for(int j=0; j<dim_; j++) {
			double data = 0.0;
			for(int k=0; k<sapNum_; k++)
				data +=  (*(*(meanAryTrs+i)+k))* (*(*(orAry+k)+j));
			data /= sapNum_-1;
			*(*(covAry+i)+j) = data;
		}
	}
}

void COV::showCovMatrix() {
	std::cout << "============================================================" << '\n';
	std::cout << "初始样本集矩阵的协方差矩阵 :" << '\n';
	for(int i=0; i<dim_; i++) {
		for(int j=0; j<dim_; j++)
			std::cout << *(*(covAry+i)+j) << "    ";
		std::cout << '\n';
	}			
}

 covariance.cpp

#include <iostream>
#include "covFile.h"

int main() {
	COV cov;
	//pca.PCA();
	cov.rdOrMatrix(10,3); //初始样本集矩阵
	cov.showOrMatrix();
	cov.meanVal(); //均值数组
	cov.showMeanVal();
	cov.meanMatrix(); //减去均值后的样本集矩阵
	cov.showMeanMatrix();
	cov.meanMatrixTrs(); //减去均值后的样本矩阵的转置矩阵
	cov.showMeanMatrixTrs();
	cov.covMatrix(); //协方差矩阵
	cov.showCovMatrix();
	
	return 0;
}