svd.cpp

一个可分解上万阶稀疏矩阵的SVD算法

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

#include "stdafx.h"

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <fstream.h>
#ifndef EPSILON
#define EPSILON 0.000000001
#endif
 
void inputtext();
void decompose();
void shrinkvector();
void outputtext(double *S, double *U, double *V);
int symeigen(double *eigen,double *AN,double *U);
extern int singdec( double *A, double *singular,
        		     double *U, double *V);
void computVS();

struct element
{
	int row;   
	int col;
	double value;
 
};
typedef struct element ELEMENT;

struct elementchain
{
	ELEMENT data;
	struct elementchain *next;
};
typedef struct elementchain ELEMENTCHAIN;


int mm;
int nn;
int precision=9;
double *finalVS;
double *A;
double *finalU;
double *finalV;
double *finalS;
int finalR;
ELEMENTCHAIN *inputmatrix;





//=========================================================================================
void main(int argc, char* argv[])
{
   inputtext();
   decompose();
   computVS();
   shrinkvector();
   outputtext(finalS,finalU,finalV);
}







//========================================================================================
int singdec( double *A, double *singular,
            double *U, double *V)
{
int i,j,jj,l,rank;
int rowA,colA,valueA;
static double sum,*vj,*eigen,*ATA,*U1;
ELEMENTCHAIN *test,*matrixA,*colI,*colJ,*tempcolI,*tempcolJ,*nextI,*nextJ;
int ci,cj;
  
vj=new double[mm];
U1=new double[nn*nn];
ATA=new double[nn*nn];
 
eigen = new double[nn];
 
for(i=0;i<nn;i++)
  {
  for(j=0;j<nn;j++)
  {
  
	ci=i;
	cj=j;
   
    colI=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
    colJ=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
	test=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
    nextI=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
    nextJ=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
    nextI=colI;
    nextJ=colJ;
 	matrixA=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
 	matrixA=inputmatrix->next;
	//把一个链表分解成两个链表；
    while (matrixA!=NULL)
	{
	   if (matrixA->data.col==ci)
		{
		   tempcolI=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
		   tempcolI->data=matrixA->data;
		   nextI->next=tempcolI;
		   nextI=tempcolI;
		}

	    if (matrixA->data.col==cj)
		{
		   tempcolJ=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
		   tempcolJ->data=matrixA->data;
	       nextJ->next=tempcolJ;
		   nextJ=tempcolJ;
	 	}
		matrixA=matrixA->next;
		test=colI;
	}
       nextI->next=NULL;
       nextJ->next=NULL;
 	   //计算ATA矩阵；
       tempcolI=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
       tempcolI=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
       tempcolI=colI->next;
       tempcolJ=colJ->next;  
 	   ATA[i*nn+j]=0.0;
  while((tempcolI!=NULL)&&(tempcolJ!=NULL))
   {
    if (tempcolI->data.row > tempcolJ->data.row)
		tempcolJ=tempcolJ->next;
	else if (tempcolI->data.row==tempcolJ->data.row)
	   {
	  	ATA[i*nn+j]+=tempcolI->data.value*tempcolJ->data.value;
	    tempcolI=tempcolI->next;   
		tempcolJ=tempcolJ->next;
	}
	else tempcolI=tempcolI->next;
   }
 
  }
  }

rank=symeigen( eigen, ATA, U );
for(i=0;i<rank;i++) singular[i]=sqrt(eigen[i]);
for(j=0;j<rank;j++)
  for(i=0;i<nn;i++) U1[i*nn+j]=U[i*nn+j]/singular[j];
for(i=0;i<mm;i++)
  {
  for(j=0;j<rank;j++)
    {
    V[i*mm+j]=0.0;
    for(l=0;l<nn;l++) V[i*mm+j]=V[i*mm+j]+A[i*nn+l]*U1[l*nn+j];
    }
  }
for(j=rank;j<mm;j++)
  {
  sum=0.0;                               /*for creat vj[]*/

  for(i=0;i<mm;i++) { vj[i]=rand(); sum=sum+vj[i]*vj[i]; }
  sum=sqrt(sum);
  for(i=0;i<nn;i++) vj[i]=vj[i]/sum;      /* vj[] is created*/
  for(jj=0;jj<j;jj++)
    {
    sum=0.0;
    for(i=0;i<mm;i++) sum=sum+vj[i]*V[i*mm+jj];
    for(i=0;i<mm;i++) vj[i]=vj[i]-sum*V[i*mm+jj];
    }
  sum=0.0;
  for(i=0;i<mm;i++) sum=sum+vj[i]*vj[i];
  sum=sqrt(sum);
  for(i=0;i<mm;i++) V[i*mm+j]=vj[j]/sum;
  }
return( rank );
}


/***************By using power and inverse power method to compute all
  eigenvalues and eigenvectos of symmetric matrix eigens. ***********/
double invpower(int init, double *fx0, double *B );
void power( int fs, double *fx0, double *A );
double lamb1,*V;





//================================================================================
int symeigen(double *eigen,double *A,double *V0)
{
int i,l,j,step,rank;
static double *B,*x0,mlamb;

V=new double[nn*nn];
x0=new double[nn];
B=new double[nn*nn];
for(i=0;i<nn;i++)
  { eigen[i]=0.0; for(j=0;j<nn;j++) V[i*nn+j]=0.0; }
rank=0;
for(step=0;step<nn;step++)
  {
  power(step,x0,A);
  if(lamb1==0.0) { eigen[step]=0.0; goto bbb; }
  mlamb=lamb1;
  for(i=0;i<nn;i++)
    {
    for(j=0;j<nn;j++) B[i*nn+j]=A[i*nn+j];
    B[i*nn+i]=A[i*nn+i]-lamb1;
    }
  lamb1=invpower( 1, x0, B );
  if( fabs(lamb1) > EPSILON ) eigen[step]=mlamb+1.0/lamb1;
  if(lamb1==0.0) eigen[step]=mlamb;
bbb:if( fabs(eigen[step]) > EPSILON ) rank=rank+1;
  for(i=0;i<nn;i++) V[i*nn+step]=x0[i];
  }
for(i=0;i<nn;i++)
   for(j=0;j<nn;j++) V0[i*nn+j]=V[i*nn+j];
return(rank);
}







//===========================================================================================
void power(int i0,double *x0,double *A)
{
static double *x1,lamb,sum,err,err1,eps;
int i,j,jj,ii;
eps=EPSILON*1000.0;
raa:ii=0;
x1=new double[nn];      
/*for creat x0[]*/
for(j=0;j<i0;j++) { for(i=0;i<nn;i++) x0[0]=rand(); }
sum=0.0;
for(i=0;i<nn;i++)
  { x0[i]=rand(); sum=sum+x0[i]*x0[i]; }
sum=sqrt(sum);
if( sum < EPSILON ) goto raa;
for(i=0;i<nn;i++) x1[i]=x0[i]/sum;      /* x0[] is created*/
lamb=0.5;
goto bbb;
aaa:ii=ii+1;
    for(j=0;j<nn;j++)
      {
      x1[j]=0.0;
      for(jj=0;jj<nn;jj++) x1[j]=x1[j]+A[j*nn+jj]*x0[jj];
      }
    lamb1=0.0;
    for(j=0;j<nn;j++) lamb1=lamb1+x0[j]*x1[j];
    if( fabs(lamb1) < EPSILON )
       {
       lamb1=0.0; err=0.0; err1=0.0;
       for(j=0;j<nn;j++) x1[j]=x0[j];
       goto bbb;
       }
    else
       { err=fabs(lamb1-lamb); err1=err/fabs(lamb1); }
    lamb=lamb1;
bbb:if(i0>0)
      {
      for(i=0;i<i0;i++)
        {
        sum=0.0;
        for(j=0;j<nn;j++) sum=sum+x1[j]*V[j*nn+i];
        for(j=0;j<nn;j++) x1[j]=x1[j]-sum*V[j*nn+i];
        }
      }
sum=0.0;
for(j=0;j<nn;j++) sum=sum+x1[j]*x1[j];
sum=sqrt(sum);
if( sum < EPSILON ) goto endd;
for(j=0;j<nn;j++) x0[j]=x1[j]/sum;
if( ii < 10 ) goto aaa;
if( ( err > eps ) || (err1 > eps) ) goto aaa;
endd:;
}






/******************************************************************======================**/
double invpower( int init, double *x0, double *B )
{
int i,j,k,i0,istep,*ik;
static double *x1,lamb,lamb1,sum,err,temp;
x1=new double[nn];
ik=new int[nn];
istep=0;
if( init == 0 )              /* for creat x0[] */
  {

raa: sum=0.0;
  for(i=0;i<nn;i++)
    { x0[i]=rand(); sum=sum+x0[i]*x0[i]; }
  sum=sqrt(sum);
  if( sum < EPSILON ) goto raa;
  for(i=0;i<nn;i++) x0[i]=x0[i]/sum;      /* x0[] is created*/
  }
for(k=0;k<nn;k++)                /*colum pivot Doolittle decompsition*/
  {
  temp=fabs( B[k*nn+k] );       ik[k]=k;
  for(i=k;i<nn;i++)
    if( fabs(B[i*nn+k]) > temp ) { temp=fabs(B[i*nn+k]); ik[k]=i; }
  i0=ik[k];
  if( i0 != k )
    {
    for(j=0;j<nn;j++)
      { temp=B[k*nn+j]; B[k*nn+j]=B[i0*nn+j]; B[i0*nn+j]=temp; }
    }
  if( fabs(B[k*nn+k]) < EPSILON ) { lamb1=0.0; goto end; }
  for(i=k+1;i<nn;i++)
    {
    B[i*nn+k]=B[i*nn+k]/B[k*nn+k];
    for(j=k+1;j<nn;j++) B[i*nn+j]=B[i*nn+j]-B[i*nn+k]*B[k*nn+j];
    }
  }                       /* decompsition is ended */
aaa:istep=istep+1;
    lamb=lamb1;
    for(j=0;j<nn;j++) x1[j]=x0[j];
    for(i=0;i<nn;i++)                             /* solut LR=Px0[] */
      { j=ik[i]; temp=x0[i]; x0[i]=x0[j]; x0[j]=temp; }
    for(i=1;i<nn;i++)
      for(j=0;j<i;j++) x0[i]=x0[i]-B[i*nn+j]*x0[j];
    for(i=nn-1;i>=0;i--)
      {
      for(j=i+1;j<nn;j++) x0[i]=x0[i]-B[i*nn+j]*x0[j];
      x0[i]=x0[i]/B[i*nn+i];
      }                                         /* end for solution */
    lamb1=0.0;
    for(j=0;j<nn;j++) lamb1=lamb1+x0[j]*x1[j];
    err=fabs(lamb1-lamb)/fabs(lamb1);
    sum=0.0;
    for(j=0;j<nn;j++) sum=sum+x0[j]*x0[j];
    sum=sqrt(sum);
    for(j=0;j<nn;j++) x0[j]=x0[j]/sum;
    if( istep > 2000 )
      {
      printf("Iterative 2000 times! Strike any key to exit.\n");
 
      exit(1);
      }
    if( err>100.0*EPSILON ) goto aaa;
end: return(lamb1) ;
}





//====================================================================================
void inputtext()
{
    //从文本中读入矩阵；分出行和列，并把矩阵存入数组；	 
    double AA;
    int i;
	int rownum,colnum;
	ELEMENTCHAIN *tempelement;
	ELEMENTCHAIN *nextelement;

	ifstream fileA("c:\\SVD\\A.txt");
	ofstream filetemp("c:\\SVD\\tempfile.txt");

	fileA>>AA;
	mm=int(AA);

	fileA>>AA;
	nn=int(AA);

	inputmatrix=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
    nextelement=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
//	tempelement=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));

	nextelement=inputmatrix;

	A=new double[mm*nn];
    i=0;
   
           for (rownum=0;rownum<mm;rownum++)
			   for (colnum=0;colnum<nn;colnum++)
			   {
		   
		              fileA>>AA;
	                  A[i]=AA;
		              i++;

					 if (AA !=0)
					  {
						 
						 tempelement=(ELEMENTCHAIN *)malloc(sizeof(ELEMENTCHAIN));
						 tempelement->data.row=rownum;
						 tempelement->data.col=colnum;
    					 tempelement->data.value=AA;
                    	 nextelement->next=tempelement;
						 nextelement=tempelement;
						 
                      }

				}

			   	  nextelement->next=NULL;
				  nextelement=inputmatrix->next;

                        while(nextelement!=NULL)
						{
                
						 filetemp<<nextelement->data.row
							     <<"   "
							     <<nextelement->data.col
								 <<"   "
								 <<nextelement->data.value
								 <<"   "
								 <<"\n";
			                nextelement=nextelement->next;


						} 
}






//==========================================================================================
void decompose()
{

 static double *U,*V,*singular; 
 U=new double[nn*nn];
 V=new double[mm*mm];
 singular=new double[nn];
 finalS=new double[nn];
 finalU=new double[nn*nn];
 finalV=new double[mm*mm];
 
 int i,j,jj,l,rank;

 system("cls"); //文件初始化；

 rank=singdec( A, singular, U, V );
 finalR=rank;
 
 for(i=0;i<nn;i++) 
  	 finalS[i]=singular[i];
  
 for(i=0;i<mm;i++)
 	 for(j=0;j<mm;j++)  
         finalV[i*mm+j]=V[i*mm+j];
  
 for(i=0;i<nn;i++)
  for(j=0;j<nn;j++) 
   finalU[i*nn+j]=U[i*nn+j];
   
}





//======================================================================================

void outputtext(double *S, double *U, double *V)
{
	//分别输出三个数组到三个文件中；
    int i,j;

	ofstream arrayS("c:\\SVD\\S.txt");
	arrayS<<"The rank="
	   <<finalR
	   <<"\n";
	arrayS.precision(precision);
	arrayS<<"The singulars:\n";
	 
    for(i=0;i<nn;i++) 
	 arrayS<<finalS[i]
	       <<"     ";
	 
    ofstream arrayV("c:\\SVD\\V.txt");
	arrayV<<"The matrix V:\n";
	arrayV.precision(precision);
    for(i=0;i<mm;i++)
	{
	 for(j=0;j<mm;j++)  
	 {
	    arrayV<<finalV[i*mm+j]
	            <<"     ";
	    if (finalV[i*mm+j]<0)
		  arrayV<<" "; 
	 }
	  arrayV<<"\n";
	}
   
   ofstream arrayU("c:\\SVD\\U.txt");
   arrayU<<"The matrix U:\n";
   arrayU.precision(precision);
   for(i=0;i<nn;i++)
   {
      for(j=0;j<nn;j++) 
	  {
		 arrayU<<finalU[i*nn+j]
		       <<"     ";
	      if (finalU[i*nn+j]<0)
			 arrayU<<" ";
	  }
         arrayU<<"\n";
   }

   printf("Complete decomposed!\n");

}






//=======================================================================================
void computVS()
{
   int i,j;
   double *tempS,*matrixS;

   ofstream arrayVS("c:\\SVD\\VS.txt");	
   
   finalVS=new double[mm*finalR];
   tempS=new double[mm*nn];
   matrixS=new double[nn*mm];

   for(i=0;i<mm;i++)
	   for(j=0;j<finalR;j++)
	   {
		      finalVS[i*finalR+j]=finalV[i*finalR+j]/finalS[j];
		      
	   }
  
}





//====================================================================================
void shrinkvector()
{
	double *xin,*xout;

	xin=new double[mm];
    xout=new double[finalR];
	ifstream fileIn("c:\\SVD\\Xin.txt");
	ofstream fileOut("c:\\SVD\\Xout.txt");
    
	for (int i=0;i<mm;i++)
		fileIn>>xin[i];
    
	for(i=0;i<finalR;i++)
	{
	   xout[i]=0;
       for(int j=0;j<nn;j++)
          xout[i]+=xin[j]*finalVS[j*nn+i];
    }

	for(i=0;i<finalR;i++)
		fileOut<<xout[i]
		       <<".......";
}