PCA算法实现
bunny.h
#include <iostream>
double **orAry; //初始样本集矩阵(减去均值后的新样本矩阵)
int sapNum_, dim_; //样本个数;维数;降维后维数
class BUNNY {
public :
BUNNY(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 itVtAry(void); //初始化样本向量矩阵vectAry(每一列为一特征向量)
void calValVect(double _e); //计算特征值和特征向量,特征值为valAry个对角线元素,对应的特征向量为vectAry各列
void stValVect(void); //按特征值由大到小将特征值和相应的特征向量排序
//获取模式矩阵和降维后的样本集矩阵
void lowDimtMatrixTrs(int _nwDim); //计算模式矩阵的转置矩阵并存入数组 lowDimtAryTrs(_nwDimt行,每一行存放一个特征矢量; _dimt列)
void tgtMatrix(void); //计算经降维处理后的新矩阵并存入数组 tgtAry(_nwDimt行(与原始矩阵相比,维数减少), _sapNum列)
void tgtMatrixTrs(void); //将降维后的样本集矩阵转置成sapNum_行,nwDim_列
void wtTgtMatrix(void); //将降维后的向量集输出到文本文件
//输出协方差相关矩阵
void showOrMatrix(void); //输出初始样本集矩阵
void showMeanVal(void); //输出各维度的均值
void showMeanMatrix(void); //输出减去均值后的样本矩阵
void showMeanMatrixTrs(void); //输出减去均值后的样本矩阵的转置矩阵
void showCovMatrix(void); //输出协方差矩阵
//输出特征值和特征向量相关矩阵
void sortVect(void); //输出排序后的特征向量矩阵
//输出降维后的样本向量集矩阵
void showTgtMatrix(void); //输出降维后的目标矩阵
private :
//计算协方差相关变量
int nwDim_; //样本个数;维数;降维后维数
double *meanValAry, **meanAryTrs; //各维度的均值;减去均值后的矩阵的转置矩阵
double **covAry; //协方差方阵(即为计算特征向量的初始矩阵)
//计算协方差的特征值和特征向量相关变量
double **vectAry; //和特征向量方阵
double e_;
//与模式矩阵和降维后的矩阵相关的变量
double **lowDimtAryTrs, **tgtAry; //模式矩阵的转置矩阵;降维后的样本向量集矩阵
};
bunnyFile.h
#include <iostream>
#include <malloc.h>
#include <math.h>
#include "bunny.h"
char fileName[10], nwFileName[10];
BUNNY::BUNNY() {
sapNum_ = 0; dim_ = 0; nwDim_ =0; e_ = 0.0;
orAry = NULL; tgtAry = NULL;
meanValAry = NULL; meanAryTrs = NULL;
covAry = NULL; vectAry = NULL;
lowDimtAryTrs = NULL;
}
//////////////////////////////////////////////////////
void BUNNY::clear() {
}
/////////////////////////////////////////////////////////////////
void BUNNY::rdOrMatrix(int _sapNum, int _dim) { //读取初始样本集矩阵 orAry
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(fileName,"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(ch!=EOF) {
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++;
}
std::cout << "read all the vectors" << '\n';
fclose(fp);
}
void BUNNY::meanVal() { //各维的均值数组 meanValAry
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 BUNNY::meanMatrix() { //减去均值后的样本集矩阵 orAry
for(int j=0; j<dim_; j++)
for(int i=0; i<sapNum_; i++)
*(*(orAry+i)+j) = *(*(orAry+i)+j)- *(meanValAry+j);
}
void BUNNY::meanMatrixTrs() { //减去均值后的样本集矩阵的转置矩阵 meanAryTrs
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 BUNNY::covMatrix() { //协方差矩阵 covAry
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;
}
}
std::cout << "输出协方差矩阵" << '\n';
for(i=0; i<dim_; i++) {
for(int j=0; j<dim_; j++)
std::cout << *(*(covAry+i)+j) << " ";
std::cout << '\n';
}
}
void BUNNY::itVtAry() { //初始化特征向量矩阵 vectAry
vectAry = (double **)malloc(dim_*sizeof(double));
for(int i=0; i<dim_; i++)
*(vectAry+i) = (double *)malloc(dim_*sizeof(double)); //特征向量矩阵为vectAry,每一列为响应特征值对应的特征向量;矩阵的初始状态为单位矩阵
for(i=0; i<dim_; i++)
for(int j=0; j<dim_; j++) {
if(i==j)
*(*(vectAry+i)+j) = 1;
else *(*(vectAry+i)+j) = 0;
}
}
void BUNNY::calValVect(double _e) { //计算协方差矩阵的特征值(在协方差方阵covAry自身中计算)和特征向量矩阵(vectAry)
e_ = _e;
//int p, q; //满足|*(*(valAry+p)+q)|=max|*(*(valAry+i)+j)|(i!=j)
double b, c, cn, sn; //b=(*(*(valAry+p)+p)- *(*(valAry+q)+q))/2, c=sqrt(b^2+(*(*(valAry+p)+q))^2)
//cn=(1/2+|b|/(2c))*(1/2), sn=+(*(*(valAry+p)+q)/(2c*cn))(b<=0) sn=-(*(*(valAry+p)+q)/(2c*cn))(b>0)
//但是当*(*(valAry+p)+p)=*(*(valAry+q)+q)是cn=sn=sqrt(2)/2
int count =1; //记录转换次数
double sum = 0.0; //计算非对角线元素平方和,与e_比较判断结束计算条件
for(int i=0; i<dim_; i++)
for(int j=0; j<i; j++)
if(i!=j)
sum += *(*(covAry+i)+j)* (*(*(covAry+i)+j));
int num=0;
while(num<800) { //num=800(fandisk 6475行) num=80(woman 112201行) num=80(bunny 35974行)
double max = 0.0;
int p=0, q=0; //满足|*(*(valAry+p)+q)|=max|*(*(valAry+i)+j)|(i!=j)
//找到满足|*(*(valAry+p)+q)|=max|*(*(valAry+i)+j)|(i!=j)的p,q
for(i=0; i<dim_; i++)
for(int j=0; j<dim_; j++)
if(i!=j)
if(fabs(*(*(covAry+i)+j)) > fabs(max)) {
max = *(*(covAry+i)+j);
p = i; q = j;
}
//std::cout << "p=" << p << ";q=" << q <<'\n';
//计算用于建立旋转矩阵Q的 cn,sn值
//std::cout << *(*(covAry+p)+p) << ',' << *(*(covAry+q)+q) <<'\n';
if( (double)(*(*(covAry+p)+p))== (double)(*(*(covAry+q)+q)) || (int)(*(*(covAry+p)+p))== (int)(*(*(covAry+q)+q))) {
cn = sqrt(2)/2;
sn = sqrt(2)/2;
}
else {
b = ( *(*(covAry+p)+p)- (*(*(covAry+q)+q)) )/2;
c = sqrt(b*b- *(*(covAry+p)+q)* (*(*(covAry+p)+q)));
cn = 1.0/2+ fabs(b)/(2*c);
if(b<=0)
sn = *(*(covAry+p)+q)/(2*c*cn);
else sn = -1* (*(*(covAry+p)+q)/(2*c*cn));
}
//std::cout << "cn=" << cn << ";sn=" << sn << '\n';
//计算旋转变化后的valAry方阵
double nwPP = *(*(covAry+p)+p)* (cn*cn) + *(*(covAry+q)+q)* (sn*sn) - *(*(covAry+p)+q)*sn*cn*2;
double nwQQ = *(*(covAry+p)+p)* (sn*sn) + *(*(covAry+q)+q)* (cn*cn) + *(*(covAry+p)+q)*sn*cn*2;
double nwPQ = ( *(*(covAry+p)+p) - *(*(covAry+q)+q) )*sn*cn + *(*(covAry+p)+q)*(cn*cn-sn*sn);
//if(fabs(nwPP) < e_) nwPP = 0.0;
//if(fabs(nwQQ) < e_) nwQQ = 0.0;
//if(fabs(nwPQ) < e_) nwPQ = 0.0;
double **temp; //存储两行临时值p,q行
temp = (double **)malloc(2*sizeof(double));
for(i=0; i<2; i++)
*(temp+i) = (double *)malloc(dim_*sizeof(double));
for(i=0; i<dim_; i++)
if(i!=p && i!=q) {
*(*(temp+0)+i) = *(*(covAry+p)+i)*cn - *(*(covAry+q)+i)*sn; //计算新p行元素
//if(fabs(*(*(temp+0)+i)) < e_) *(*(temp+0)+i) = 0.0;
*(*(temp+1)+i) = *(*(covAry+p)+i)*sn + *(*(covAry+q)+i)*cn; //计算新q行元素
//if(fabs(*(*(temp+1)+i)) < e_) *(*(temp+1)+i) = 0.0;
}
*(*(covAry+p)+p) = nwPP; //更新PP行
*(*(covAry+q)+q) = nwQQ; //更新QQ行
*(*(covAry+p)+q) = nwPQ; //更新PQ行
*(*(covAry+q)+p) = nwPQ; //更新QP行
for(i=0; i<dim_; i++)
if(i!=p && i!=q) {
*(*(covAry+p)+i) = *(*(temp+0)+i); //更新p行
*(*(covAry+i)+p) = *(*(temp+0)+i); //更行p列
*(*(covAry+q)+i) = *(*(temp+1)+i); //更新q行
*(*(covAry+i)+q) = *(*(temp+1)+i); //更行q列
}
delete []temp;
/*
//检验更新是否正确
std::cout << "转换" << count << "次之后的特征值矩阵" << '\n';
for(i=0; i<dim_; i++) {
for(int j=0; j<dim_; j++)
std::cout << *(*(covAry+i)+j) << " ";
std::cout << '\n';
}
*/
//计算旋转变换之后的特征向量方阵Q
double **tpVect; //存储两列临时值p,q列
tpVect = (double **)malloc(dim_*sizeof(double));
for(i=0; i<dim_; i++)
*(tpVect+i) = (double *)malloc(2*sizeof(double));
for(i=0; i<dim_; i++) {
*(*(tpVect+i)+0) = *(*(vectAry+i)+p)*cn - *(*(vectAry+i)+q)*sn; //存放新p列
//if(fabs(*(*(tpVect+i)+0)) < e_) *(*(tpVect+i)+0) = 0.0;
*(*(tpVect+i)+1) = *(*(vectAry+i)+p)*sn + *(*(vectAry+i)+q)*cn; //存放新q列
//if(fabs(*(*(tpVect+i)+1)) < e_) *(*(tpVect+i)+1) = 0.0;
}
for(i=0; i<dim_; i++) {
*(*(vectAry+i)+p) = *(*(tpVect+i)+0); //更新p列
*(*(vectAry+i)+q) = *(*(tpVect+i)+1); //更新p列
}
delete []tpVect;
/*
//检验一下更行后的特征向量矩阵
std::cout << "转换" << count << "次之后的特征向量矩阵" << '\n';
for(i=0; i<dim_; i++) {
for(int j=0; j<dim_; j++)
std::cout << *(*(vectAry+i)+j) << " ";
std::cout << '\n';
}
count++;
*/
//进行下一轮运算
num++;
sum = 0.0;
for(i=0; i<dim_; i++) //计算非对角线元素平方和,与e_比较判断结束计算条件
for(int j=0; j<i; j++)
if(i!=j)
sum += *(*(covAry+i)+j)* (*(*(covAry+i)+j));
}//end while
}
void BUNNY::stValVect() { //按特征值由大到小排序并将相应的特征向量矩阵排序 vectAry
double *sort; //特征向量排序时存储数据的临时数组
sort = (double *)malloc(dim_*sizeof(double));
for(int i=0; i<dim_-1; i++) {
int flag=i;
double max = *(*(covAry+i)+i);
for(int j=i+1; j<dim_; j++)
if(*(*(covAry+j)+j) > max) {
max = *(*(covAry+j)+j);
flag=j;
}
if(flag!=i) {
for(int k=0; k<dim_; k++) { //将相应的特征向量交换
*(sort+k)= *(*(vectAry+k)+i);
*(*(vectAry+k)+i)= *(*(vectAry+k)+flag);
*(*(vectAry+k)+flag)= *(sort+k);
}
//将相应的特征值进行交换,否则下一次比较会出错
double tmp = *(*(covAry+i)+i);
*(*(covAry+i)+i)= *(*(covAry+flag)+flag);
*(*(covAry+flag)+flag)= tmp;
}
}
}
void BUNNY::lowDimtMatrixTrs(int _nwDim) { //降维后的模式矩阵的转置矩阵 lowDimtAryTrs (直接转置即可,不需开辟新内存存储模式矩阵)
//转置矩阵的每一行为一个特征向量,列数为维数dim_
nwDim_ = _nwDim;
lowDimtAryTrs = (double **)malloc(nwDim_*sizeof(double));
for(int i=0; i<nwDim_; i++)
*(lowDimtAryTrs+i) = (double *)malloc(dim_*sizeof(double));
for(i=0; i<nwDim_; i++)
for(int j=0; j<dim_; j++)
*(*(lowDimtAryTrs+i)+j) = *(*(vectAry+j)+i);
}
void BUNNY::tgtMatrix() { //降维后的新矩阵 tgtAry (nwDim_*sapNum_)
tgtAry = (double **)malloc(nwDim_*sizeof(double));
for(int i=0; i<nwDim_; i++)
*(tgtAry+i) = (double *)malloc(sapNum_*sizeof(double));
//计算方法:模式矩阵的转置矩阵(lowDimtAryTrs(nwDim_*dim_))*减去均值后的样本集矩阵的转置矩阵(meanAryTrs(dim_*sapNum_))
for(i=0; i<nwDim_; i++)
for(int j=0; j<sapNum_; j++) {
double data = 0.0;
for(int k=0; k<dim_; k++)
data += *(*(lowDimtAryTrs+i)+k) * (*(*(meanAryTrs+k)+j));
*(*(tgtAry+i)+j) = data;
}
}
void BUNNY::tgtMatrixTrs() { //将降维后的样本集矩阵转置存放在矩阵orAry(sapNum_行dim_列)中(注意浪费了一列)
for(int i=0; i<sapNum_; i++)
for(int j=0; j<nwDim_; j++)
*(*(orAry+i)+j)= *(*(tgtAry+j)+i);
}
void BUNNY::wtTgtMatrix() {
FILE *fp;
if((fp=fopen(nwFileName,"w"))==NULL) {
printf("Cannot build file strike any key exit!");
exit(0);
}
for(int i=0; i<sapNum_; i++) {
fprintf(fp,"%s"," v ");
for(int j=0; j<nwDim_; j++)
fprintf(fp,"%.4f%s",*(*(orAry+i)+j)," ");
fprintf(fp,"%s","\n");
}
fclose(fp);
std::cout << "write all the vectors" << '\n';
}
void BUNNY::showTgtMatrix() { //输出降维后的矩阵的转置矩阵 orAry
std::cout << "============================================================" << '\n';
std::cout << "降维后矩阵的转置矩阵(每一行为一个样本向量,列数为维数,行数为样本个数) :" << '\n';
for(int i=0; i<sapNum_; i++) {
for(int j=0; j<nwDim_; j++)
std::cout << *(*(orAry+i)+j) << " ";
std::cout << '\n';
}
}
bunny.cpp
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include "bunnyFile.h"
#define er 1.0E-4
////////////////////////////////////
#include <GL/ogl.h>
#define STEP 0.01
///////////////////////////////////////////////////////////////////////////////
double angle, rx, ry, rz; //旋转设定
double bx, by, bz; //模型放大倍数
double x, y, z; //模型在相应轴上、下移动的像素数
double eyex, eyez, upx; //设置lookAt函数变换观察位置
///////////////////////////////////////////////////////////////////////////////
void myDrawing(void);
void myKeyboard(unsigned char _key, int _x, int _y);
void reshape(GLsizei _wid, GLsizei _hei);
///////////////////////////////////////////////////////////////////////////////
int main(int argc, char **argv) {
rx = ry= rz = 0;
angle = 0;
bx = by = bz =0;
x = y = z =0;
eyex = eyez =0; upx =0;
std::cout << "sapNum_ = 112201//woman.txt 6475//fandisk.txt 35947//bunny.txt" << '\n';
std::cout << "input the fileName: ";
std::cin >> fileName;
std::cout << "input the nwFileName: ";
std::cin >> nwFileName;
int sapNum;
std::cout << "input the sapNum of the file: ";
std:: cin >> sapNum;
int dim;
std::cout << "input the dim of the file: ";
std:: cin >> dim;
BUNNY bunny;
bunny.rdOrMatrix(sapNum,dim); //初始样本集矩阵 sapNum_ = 112201//woman 6475//fandisk 35947//bunny
bunny.meanVal(); //均值数组
bunny.meanMatrix(); //减去均值后的样本集矩阵
bunny.meanMatrixTrs(); //减去均值后的样本矩阵的转置矩阵
bunny.covMatrix(); //协方差矩阵
bunny.itVtAry(); //初始化特征向量矩阵
bunny.calValVect(er); //计算协方差矩阵的特征值和特征向量矩阵
bunny.stValVect(); //对特征向量排序
bunny.lowDimtMatrixTrs(2); //计算模式矩阵的转置矩阵并存入数组 lowDimtAryTrs(_nwDimt行,每一行存放一个特征矢量; _dimt列)
bunny.tgtMatrix(); //计算经降维处理后的新矩阵并存入数组 tgtAry(_nwDimt行(与原始矩阵相比,维数减少), _sapNum列)
bunny.tgtMatrixTrs(); //降维后矩阵转置
bunny.wtTgtMatrix(); //输出转置后的降维矩阵
glutInitDisplayMode(GLUT_SINGLE|GLUT_RGB);
glutInitWindowSize(500, 500);
glutInitWindowPosition(400, 0);
glutCreateWindow("Hello");
glutDisplayFunc(myDrawing);
glutKeyboardFunc(myKeyboard);
glutReshapeFunc(reshape);
glutMainLoop();
return 0;
}
///////////////////////////////////////////////////////////////////////////////
void reshape(GLsizei _wid, GLsizei _hei) {
glViewport(0, 0, _wid, _hei);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glMatrixMode(GL_MODELVIEW);
}
///////////////////////////////////////////////////////////////////////////////
// 图形绘制函数
void myDrawing(void) {
glClearColor(1.0, 1.0, 1.0f, 0.0); //设置背景清除蓝色
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //执行清除
glPushMatrix();
glLoadIdentity();
gluLookAt(eyex, 0, eyez, 0, 0, -1, upx, 1, 0); //会在modelview矩阵上叠加
glRotatef(angle,rx,ry,rz);
glTranslatef(x, y, z);
glScalef(bx,by,bz);
glScalef(0.01,0.01,0.01); //fandisk时需要另加此句(与坐标大小有关!!!)
glPointSize(1);
glBegin (GL_POINTS);
for(int i=0; i<sapNum_; i++) {
int k=0;
glColor3f(1.0,0.0,0.0);
glVertex2f(*(*(orAry+i)+k),*(*(orAry+i)+k+1));
}
glEnd();
//gluLookAt(eyex, 0, eyez, 0, 0, -1, upx, 1, 0); //会在modelview矩阵上叠加----顺序
glPopMatrix();
glFlush();
}
///////////////////////////////////////////////////////////////////////////////
// 键盘事件处理函数 <==========NEW!!!
void myKeyboard(unsigned char _key, int _x, int _y) {
switch (_key) {
case 'q': //x向右移动
{
x++;
break;
}
case 'w': //y向上移动
{
y++;
break;
}
case 'e': //z向前移动
{
z++;
break;
}
case 'a': //x向左移动
{
x--;
break;
}
case 's': //y向下移动
{
y--;
break;
}
case 'd': //z向后移动
{
z--;
break;
}
case 'b': //放大的倍数
{
bx++; by++; bz++;
break;
}
case 'v': //缩小的倍数
{
bx--; by--; bz--;
break;
}
case 'x': //绕相应轴旋转的角度
{
angle+=30;
if(angle==360) angle=0;
rx=1.0; ry=0; rz=0;
break;
}
case 'y': //绕相应轴旋转的角度
{
angle+=30;
if(angle==360) angle=0;
rx=0.0; ry=1.0; rz=0;
break;
}
case 'z': //绕相应轴旋转的角度
{
angle+=30;
if(angle==360) angle=0;
rx=0.0; ry=0; rz=1.0;
break;
}
case 'o': eyex += STEP; break; //观察点沿x 右移动
case 'i': eyex -= STEP; break; //观察点沿x 左移动
case 'l': eyez += STEP; break; //观察点沿z 上移动
case 'k': eyez -= STEP; break; //观察点沿z 下移动
case 'm': upx += STEP; break; //头向x 右偏
case 'n': upx -= STEP; break; //头向x 左偏
}
glutPostRedisplay();
}
浙公网安备 33010602011771号