lnmf.cpp

人工智能模式识别中基于非负矩阵分解生成特征空间的算法

// LNMF.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "time.h"

void InvertMatrix(double *Matrix, int m, int n)
{//将矩阵转置,输入矩阵为m行\n列,输出矩阵为n行\m列
	int i,j;
	double *NewMatrix = new double[n*m];
	for(i = 0; i < m; i++)
		for(j = 0; j < n ; j++)
			NewMatrix[j*m+i] = Matrix[i*n+j];
    for(i = 0; i < n; i++)
		for(j = 0; j < m; j++)
         Matrix[i*m+j] = NewMatrix[i*m+j];
	delete []NewMatrix;
}
void main()
{
	int n,N,r;								 //n为每幅人脸图像的像素点数，N为人脸图像数目,r为分解后的维数
	n = 4;
	N = 3;
	r = 2;
	double Alpha = 0.4;
	double Beta = 0.5;
	double *B = new double[n*r];             //分解后的基向量
	double *BB = new double[n*r];
	double *H = new double[r*N];             //分解后的系数矩阵
	double *HH = new double[r*N];            //分解后的系数矩阵
	double *OriginalImage = new double[n*N]; //所有人脸图像组成的矩阵
	double *NewImage = new double[n*N];      //每次更新过程中保存BH的值
	double *Binvert = new double[r*n];
	double *Hinvert = new double[N*r];
	double *U = new double[r*r];
	double *V = new double[r*r];	

	int i, j, k,l, p, s, ii, jj, kk;
	
	/////////////////给待分解矩阵赋值////////////////////////
	int temp = 1;
	for (i=0; i<n; i++)
	{
		for (j=0; j<N; j++)
		{
			OriginalImage[i*N+j] = temp;
			temp++;
		}
	}
	
	/////////////////初始非0矩阵 B，H/////////////////////////
	for (i=0; i<n; i++)
		for (j=0; j<r; j++)
			B[i*r+j] = 0.5;
	for (j=0; j<r; j++)
		for (k=0; k<N; k++)
			H[j*N+k] = 0.5;

	//////更新B，H，直到收敛
	double temp1, temp2, temp3,temp4;
	double distance;
	int count = 0;
	do
	{
		count++;
	
		for (k=0; k<r; k++)   
		{  
			for (l=0; l<N; l++)      //更新H一行
			{
				temp1 = 0.0;
				for (i=0; i<n; i++)
				{
					temp2 = 0.0;
					for (p=0; p<r; p++)
						temp2 += (B[i*r+p]*H[p*N+l]);
					temp1 = temp1 + (OriginalImage[i*N+l]*B[i*r+k])/temp2;
				}				   
				HH[k*N+l] = sqrt(H[k*N+l]*temp1);
			}
			
			for (l=0; l<N; l++)      //更新H
				H[k*N+l] = HH[k*N+l];
			
			for (s=0; s<n; s++)      //更新B一列
			{
				temp1 = 0.0;
				for (j=0; j<N; j++)
				{
					temp2 = 0.0;
					for (l=0; l<r; l++)
						temp2 += B[s*r+l]*H[l*N+j];
					temp1 = temp1 + (OriginalImage[s*N+j]*H[k*N+j])/temp2;
				}
				
				temp3 = 0.0;
				for (j=0; j<N; j++)
					temp3 += H[k*N+j];
				
				BB[s*r+k] = B[s*r+k]*temp1/temp3;			
			}
			
			for (s=0; s<n; s++)    //更新B
				B[s*r+k] = BB[s*r+k];				
			
			temp4 = 0;
			for (s=0; s<n; s++)     //更新B一列后,对一列进行归一化
				temp4 += B[s*r+k];
			for (s=0; s<n; s++)      	
				B[s*r+k] = B[s*r+k]/temp4;				
		}
		
		//计算U
		for (ii=0; ii<n*r; ii++)
			Binvert[ii] = B[ii];

		InvertMatrix(Binvert, n, r);
	
		for (ii=0; ii<r; ii++)
		{
			for (kk=0; kk<r; kk++)
			{
				temp1 = 0.0;
				for (jj=0; jj<n; jj++)
					temp1 += Binvert[ii*n+jj]*B[jj*r+kk];
				U[ii*r+kk] = temp1;
			}
		}			
		
		//计算V
		for (ii=0; ii<r*N; ii++)
			Hinvert[ii] = H[ii];
		
		InvertMatrix(Hinvert,r,N);
		
		for (ii=0; ii<r; ii++)
		{
			for (kk=0; kk<r; kk++)
			{
				temp1 = 0.0;
				for (jj=0; jj<N; jj++)
					temp1 += H[ii*N+jj]*Hinvert[jj*r+kk];
				V[ii*r+kk] = temp1;
			}
		}
		
		//计算BH
		for (ii=0; ii<n; ii++)
		{
			for (kk=0; kk<N; kk++)
			{
				temp1 = 0.0;
				for (jj=0; jj<r; jj++)
					temp1+= B[ii*r+jj]*H[jj*N+kk];
				NewImage[ii*N+kk] = temp1;
			}
		}			
	
		//计算OriginalImage和BH之间的前差异距离
		distance = 0;
		temp1 = 0;		
		for (ii=0; ii<n; ii++)
			for (jj=0; jj<N; jj++)
				temp1 += (OriginalImage[ii*N+jj]*log(OriginalImage[ii*N+jj]/NewImage[ii*N+jj]) - OriginalImage[ii*N+jj] + NewImage[ii*N+jj]);
		
		temp2 = 0;
		for (ii=0; ii<r; ii++)
			for (jj=0; jj<r; jj++)
				temp2 += U[ii*r+jj];
		
		temp3 = 0;
		for (ii=0; ii<r; ii++)
		  	temp3 += V[ii*r+ii];
        distance = temp1+Alpha*temp2-Beta*temp3;					
	
	}while((fabs(distance)>0.0001)&&(count < 1000));
	
	double *DiffX = new double[n*N];
	for (i=0; i<n; i++)
		for (j=0; j<N; j++)
			DiffX[i*N+j] = OriginalImage[i*N+j]-NewImage[i*N+j];
	
	printf("Count:%d\ndistance:%lf\n",count ,distance);
	
	printf("Original Matrix:\n");
	
	for (i=0; i<n; i++)
	{
		for (j=0; j<N; j++)
			printf("%15.8f",OriginalImage[i*N+j]);
		printf("\n");
	}
   
	printf("B:\n");
	for (i=0; i<n; i++)
	{
		for (j=0; j<r; j++)
			printf("%15.8f",B[i*r+j]);
		printf("\n");
	}

	printf("H:\n");
	for (i=0; i<r; i++)
	{
		for (j=0; j<N; j++)
			printf("%15.8f",H[i*N+j]);
		printf("\n");
	}

	printf("DiffOriginalImage:\n");
	for (i=0; i<n; i++)
	{
		for (j=0; j<N; j++)
			printf("%12.8f",DiffX[i*N+j]);
		printf("\n");
	}
	printf("\nDistance is:%lf\n",distance);

	delete []Binvert;
	delete []Hinvert;
	delete []DiffX;
	delete []U;
	delete []V;
	delete []B;
	delete []BB;
	delete []H;
	delete []HH;
}