a.cpp

基于非负矩阵分解（NMF）的人脸特征提取算法

#include<stdio.h>
#include<conio.h>
#include<stdlib.h>
#include<math.h>
#include<windows.h>

#define maxiter 80000

//-------------------------------数据类型定义-----------------------
//矩阵结构定义

#define m 10304
#define n 1			//测试时为1，训练时为7
#define k 5

#define k1 1
#define k2 1000000

double A[m][n];//原始矩阵
double W[m][k];//左乘矩阵
double H[k][n];//右乘矩阵

double H_[k][1];

//-------------------------------随机矩阵生成函数-------------------

void matrixrand(int tag)
{
#define RANDMAX 1000000000//可以运行前设置，也可=A中数据的最大值+1~……
	int i,j;
	if(tag==0)//tag为0时为左乘矩阵的初始化
		for(i=0;i<m;i++)
			for(j=0;j<k;j++)
				W[i][j]=double(rand()%RANDMAX*1.0/RANDMAX);
	else
		if(tag==1)//tag为1时为右乘矩阵的初始化
			for(i=0;i<k;i++)
				for(j=0;j<n;j++)
					H[i][j]=double(rand()%RANDMAX*1.0/RANDMAX);
#undef RANDMAX
}

//-------------------------------信息读入函数-----------------------

bool inforead(char *filename)
{
	int i,j;
	FILE *fp;
	fp=fopen(filename,"r");
	if(fp==NULL)
		return FALSE;
	for(i=0;i<n;i++)
		for(j=0;j<m;j++)
			fscanf(fp,"%lf",&A[j][i]);
	fclose(fp);
	return TRUE;
}

//-------------------------------正规化函数-------------------------

void normalization(void )
{
	int i,j;
	double max,min;
	for(i=0;i<m;i++)
	{
		max=min=A[i][0];
		for(j=0;j<n;j++)
		{
			if(A[i][j]<min)
				min=A[i][j];
			else
				if(max<A[i][j])
					max=A[i][j];
		}
		for(j=0;j<n;j++)
			A[i][j]=(A[i][j]-min)*1.0/(max-min);
	}
}

//-------------------------------归一化H----------------------------

void NormalizeH(void )
{
	int i,j;
	double sum;
	//H归一化
	for(i=0;i<k;i++)
	{
		sum=0;
		for(j=0;j<n;j++)
			sum+=H[i][j];
		for(j=0;j<n;j++)
			H[i][j]/=sum;
	}
}

//-------------------------------归一化W----------------------------

void NormalizeW(void )
{
	int i,j;
	double sum;
	//W归一化
	for(i=0;i<k;i++)
	{
		sum=0;
		for(j=0;j<m;j++)
			sum+=W[j][i];
		for(j=0;j<m;j++)
			W[j][i]/=sum;
	}
}

//-------------------------------计算H的迭代公式--------------------

void updateH(void )
{
	int i,j,con;
	double WT[k][m];//存储W的转置
	double Temp_1[k][n];//存储WT与A的积
	double Temp_2[k][k];//存储WT与W的积
	double Temp_3[k][n];//存储Temp_2与H的积+1e-6
	double Temp_4[k][n];//存储H.*Temp_1
	//转置操作W
	for(i=0;i<k;i++)
		for(j=0;j<m;j++)
			WT[i][j]=W[j][i];
	//求WT与A的积
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
		{
			Temp_1[i][j]=0;
			for(con=0;con<m;con++)
				Temp_1[i][j]+=WT[i][con]*A[con][j];
		}
	//求WT与W的积
	for(i=0;i<k;i++)
		for(j=0;j<k;j++)
		{
			Temp_2[i][j]=0;
			for(con=0;con<m;con++)
				Temp_2[i][j]+=WT[i][con]*W[con][j];
		}
	//求Temp_2与H的积
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
		{
			Temp_3[i][j]=0;
			for(con=0;con<k;con++)
				Temp_3[i][j]+=Temp_2[i][con]*H[con][j];
			Temp_3[i][j]+=1e-9;
		}
	//求H.*Temp_1
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
			Temp_4[i][j]=H[i][j]*Temp_1[i][j];
	//求一趟循环后的H
	for(i=0;i<k;i++)
		for(j=0;j<n;j++)
			H[i][j]=Temp_4[i][j]/Temp_3[i][j];
}

//-------------------------------计算W的迭代公式--------------------

void updateW(void )
{
	int i,j,con;
	double HT[n][k];//存储H的转置
	double Temp_01[m][k];//存储A与HT的积
	double Temp_02[k][k];//存储H与HT的积
	double Temp_03[m][k];//存储W与Temp_02的积+1e-6
	double Temp_04[m][k];//存储W.*Temp_01
	//转置操作H
	for(i=0;i<n;i++)
		for(j=0;j<k;j++)
			HT[i][j]=H[j][i];
	//求A与HT的积
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
		{
			Temp_01[i][j]=0;
			for(con=0;con<n;con++)
				Temp_01[i][j]+=A[i][con]*HT[con][j];
		}
	//求H与HT的积
	for(i=0;i<k;i++)
		for(j=0;j<k;j++)
		{
			Temp_02[i][j]=0;
			for(con=0;con<n;con++)
				Temp_02[i][j]+=H[i][con]*HT[con][j];
		}
	//求W与Temp_02的积
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
		{
			Temp_03[i][j]=0;
			for(con=0;con<k;con++)
				Temp_03[i][j]+=W[i][con]*Temp_02[con][j];
			Temp_03[i][j]+=1e-9;
		}
	//求W.*Temp_01
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
			Temp_04[i][j]=W[i][j]*Temp_01[i][j];
	//求一趟循环后的W
	for(i=0;i<m;i++)
		for(j=0;j<k;j++)
			W[i][j]=Temp_04[i][j]/Temp_03[i][j];
}

//-------------------------------原始算法函数-----------------------

void origNMF(void )
{
	int control=0;
	while(control<maxiter)
	{
		updateH();
		NormalizeH();//归一化H
		updateW();
//		NormalizeW();//归一化W
		control++;
	}
}

//-------------------------------判优函数---------------------------

double g_estimate(void )
{
	int i,con;
	double gestimate=0;
	double AWH[m][n];//存储W*H所生成的矩阵
	for(i=0;i<m;i++)
	{
		AWH[i][0]=0;
		for(con=0;con<k;con++)
			AWH[i][0]+=W[i][con]*H[con][0];
	}
	for(i=0;i<m;i++)
		gestimate+=k1*pow((A[i][0]-AWH[i][0]),2);
	for(i=0;i<k;i++)
		gestimate+=k2*pow((H[i][0]-H_[i][0]),2);
	return gestimate;
}

//-------------------------------测试函数---------------------------

bool test(void )
{
	int min;
	double tempmin;
	double Cgestimate[10];
	/////////////////////////////////////////////////////////////////

	char Wname[10][10]={"W1.txt","W2.txt","W3.txt","W4.txt","W5.txt","W6.txt","W7.txt","W8.txt","W9.txt","W10.txt"};
	char hname[10][10]={"h1.txt","h2.txt","h3.txt","h4.txt","h5.txt","h6.txt","h7.txt","h8.txt","h9.txt","h10.txt"};

	/////////////////////////////////////////////////////////////////
	FILE *fp,*fpresult,*fptemp;
	if(!(fp=fopen("测试.txt","r")))
		return false;
	if(!(fpresult=fopen("结果.txt","w")))
		return false;

	for(int i=0;i<40;i++)
	{
		for(int j=0;j<m;j++)
			fscanf(fp,"%lf",&A[j][0]);
		int ii,jj;

		//计算对应的目标函数
		for(int testtemp=0;testtemp<10;testtemp++)
		{
			fptemp=fopen(Wname[testtemp],"r");
			for(ii=0;ii<m;ii++)
				for(jj=0;jj<k;jj++)
					fscanf(fptemp,"%lf",&W[ii][jj]);
			fclose(fptemp);

			fptemp=fopen(hname[testtemp],"r");
			for(ii=0;ii<k;ii++)
				fscanf(fptemp,"%lf",&H_[ii][0]);
			fclose(fptemp);

			matrixrand(1);//初始化H_矩阵

			for(j=0;j<200;j++)
				updateH();

			Cgestimate[testtemp]=g_estimate();
		}
		
		//找出Cgestimate中最小的值
		tempmin=Cgestimate[0];
		min=0;
		for(int t=1;t<10;t++)
			if(Cgestimate[t]<tempmin)
			{
				min=t;
				tempmin=Cgestimate[t];
			}

		fprintf(fpresult,"%d\n",min);
	}
	fclose(fp);
	fclose(fpresult);
	return true;
}

//-------------------------------主函数-----------------------------

void main()
{
//                  测试

	
	test();

	
/*//					训练
	bool bo;
	bo=inforead("训练2.txt");		//具体训练时需要更改该处的名称
	if(bo==FALSE)
		printf("文件打开失败！！！\n");
	//矩阵的初始化
	matrixrand(0);
	matrixrand(1);

	origNMF();

	int i,j;
	FILE *fp;
	fp=fopen("W2.txt","w");			//具体训练时需要更改该处的名称
	for(i=0;i<m;i++)
	{
		for(j=0;j<k;j++)
			fprintf(fp,"%lf\t",W[i][j]);
		fprintf(fp,"\n");
	}
	fclose(fp);

	double Havg;
	fp=fopen("h2.txt","w");			//具体训练时需要更改该处的名称
	for(i=0;i<k;i++)
	{
		Havg=0;
		for(j=0;j<n;j++)
			Havg+=(H[i][j]*1.0/n);
		fprintf(fp,"%lf\n",Havg);
	}
	fclose(fp);
*/
	printf("OK!\n");
	getch();
}