pca2.cpp

对输入的高维特征向量进行pca降维后输出低维的特征向量

				{ 
					ix = (i-1)*n+kk-1; 
					v[ix] = -v[ix];
				}
			}
			while ( (kk!=m1) && (s[kk-1] < s[kk]) )
			{ 
				d = s[kk-1]; 
				s[kk-1] = s[kk]; 
				s[kk] = d;
				if ( kk < n )
					for (i=1; i<=n; i++)
					{ 
						ix = (i-1)*n+kk-1; iy = (i-1)*n+kk;
						d = v[ix]; v[ix] = v[iy]; v[iy] = d;
					}
					if ( kk < m )
						for ( i = 1; i <= m; i++ )
						{ 
							ix = (i-1)*m+kk-1; iy = (i-1)*m+kk;
							d = u[ix]; u[ix] = u[iy]; u[iy] = d;
						}
					kk = kk+1;
			}
			it = 60;
			mm = mm-1;
		}
		else
		{ 
			ks = mm;
			while ( (ks>kk) && (fabs(s[ks-1]) != 0.0 ) )
			{ 
				d = 0.0;
				if ( ks != mm ) d = d+fabs(e[ks-1]);
				if ( ks != kk+1 ) d = d+fabs(e[ks-2]);
				dd = fabs(s[ks-1]);
				if ( dd > eps*d ) ks = ks-1;
				else s[ks-1] = 0.0;
			}
			if ( ks == kk )
			{ 
				kk = kk+1;
				d = fabs(s[mm-1]);
				t = fabs(s[mm-2]);
				if ( t > d ) d = t;
				t = fabs( e[mm-2] );
				if ( t > d ) d = t;
				t = fabs( s[kk-1] );
				if ( t > d ) d = t;
				t = fabs( e[kk-1] );
				if ( t > d ) d = t;
				sm = s[mm-1]/d; sm1 = s[mm-2]/d;
				em1 = e[mm-2]/d;
				sk = s[kk-1]/d; ek = e[kk-1]/d;
				b = ((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
				c = sm*em1; c = c*c; shh = 0.0;
				if ( (b!=0.0) || (c!=0.0) )
				{ 
					shh = sqrt( b*b+c );
					if ( b < 0.0 ) shh = -shh;
					shh = c/(b+shh);
				}
				fg[0] = (sk+sm)*(sk-sm)-shh;
				fg[1] = sk*ek;
				for ( i = kk; i <= mm-1; i++ )
				{ 
					sss( fg, cs );
					if ( i != kk ) e[i-2] = fg[0];
					fg[0] = cs[0]*s[i-1]+cs[1]*e[i-1];
					e[i-1] = cs[0]*e[i-1]-cs[1]*s[i-1];
					fg[1] = cs[1]*s[i];
					s[i] = cs[0]*s[i];
					if ( (cs[0] != 1.0) || (cs[1] != 0.0) )
						for ( j = 1; j <= n; j++ )
						{ 
							ix = (j-1)*n+i-1;
							iy = (j-1)*n+i;
							d = cs[0]*v[ix]+cs[1]*v[iy];
							v[iy] = -cs[1]*v[ix]+cs[0]*v[iy];
							v[ix] = d;
						}
						sss( fg, cs );
						s[i-1] = fg[0];
						fg[0] = cs[0]*e[i-1]+cs[1]*s[i];
						s[i] = -cs[1]*e[i-1]+cs[0]*s[i];
						fg[1] = cs[1]*e[i];
						e[i] = cs[0]*e[i];
						if ( i < m )
							if ( ( cs[0] != 1.0 ) || ( cs[1] != 0.0 ) )
								for ( j = 1; j <= m; j++ )
								{ 
									ix = (j-1)*m+i-1;
									iy = (j-1)*m+i;
									d = cs[0]*u[ix]+cs[1]*u[iy];
									u[iy] = -cs[1]*u[ix]+cs[0]*u[iy];
									u[ix] = d;
								}
				}
				e[mm-2] = fg[0];
				it = it-1;
			}
			else
			{ 
				if ( ks == mm )
				{ 
					kk = kk+1;
					fg[1] = e[mm-2]; e[mm-2] = 0.0;
	                for ( ll = kk; ll <= mm-1; ll++ )
		            { 
						i = mm+kk-ll-1;
				        fg[0] = s[i-1];
					    sss( fg, cs );
						s[i-1] = fg[0];
					    if ( i != kk )
						{ 
							fg[1] = -cs[1]*e[i-2];
							e[i-2] = cs[0]*e[i-2];
						}
						if ( (cs[0]!=1.0)||(cs[1]!=0.0) )
							for ( j = 1; j <= n; j++ )
							{ 
								ix = (j-1)*n+i-1;
								iy = (j-1)*n+mm-1;
								d = cs[0]*v[ix]+cs[1]*v[iy];
								v[iy] = -cs[1]*v[ix]+cs[0]*v[iy];
								v[ix] = d;
							}
					}
				}
				else
				{ 
					kk = ks+1;
					fg[1] = e[kk-2];
					e[kk-2] = 0.0;
					for ( i = kk; i <= mm; i++ )
					{ 
						fg[0] = s[i-1];
						sss( fg, cs );
						s[i-1] = fg[0];
						fg[1] = -cs[1]*e[i-1];
						e[i-1] = cs[0]*e[i-1];
						if ( (cs[0]!=1.0)||(cs[1]!=0.0) )
							for (j=1; j<=m; j++)
							{ 
								ix = (j-1)*m+i-1;
								iy = (j-1)*m+kk-2;
								d = cs[0]*u[ix]+cs[1]*u[iy];
								u[iy] = -cs[1]*u[ix]+cs[0]*u[iy];
								u[ix] = d;
							}
					}
				}
			}
		}
    }

    return( 1 );
}

# define EPS  0.000001

/*
  PCA2: Perform PCA using SVD.
  data     --- MxN matrix of input data ( M dimensions, N trials )
  signals  --- MxN matrix of projected data 
  PC       --- each column is a PC
  V        --- Mx1 matrix of variances
  row = M dimensions, col = N trials 
*/
int pca2( double * data, int row, int col,
	      double * signals, double * PC, double * V )
{
	int x, y, i, ka, rvalue;
	double * mean;
	double * u;
    double * d;
    double * v;
	double * sd;
	double * data1;
	double sqrt_col_1;
	int col1, row1;

	if ( row <= 1 || col <= 1 ) return( -1 );

	/* subtract off the mean for each dimension */
	mean = ( double * )malloc( sizeof( double )*row );
	data1 = ( double * ) malloc( sizeof( double )*row*col );
	for ( y = 0; y < row; y++ ) /* calculate mean */
	{
		mean[y] = 0;
		for ( x = 0; x < col; x++ )
			mean[y] += data[y*col+x];
	}
	for ( y = 0; y < row; y++ ) mean[y] = mean[y]/col;
	for ( y = 0; y < row; y++ ) /* subtract mean */
		for ( x = 0; x < col; x++ )
		{
			data1[y*col+x] = data[y*col+x]-mean[y];
		}
	free( mean );

	/* construct the matrix Y: Y = data' / sqrt(col-1); */

	/* construct the matrix Y: Y = data' / sqrt(col-1); */
	sqrt_col_1 = sqrt(float(col-1));
	row1 = col;
	col1 = row;
	sd = (double *)malloc( sizeof(double)*col1*row1 );
	u = (double *)malloc( sizeof(double)*row1*row1 );
	v = (double *)malloc( sizeof(double)*col1*col1 );
	for ( y = 0; y < row1; y++ )
		for ( x = 0; x < col1; x++ )
		{
			sd[y*col1+x] = data1[x*row1+y]/sqrt_col_1;
		}

	/* SVD does it all: [u, S, PC] = svd( Y ); */
	if ( row1 >= col1 ) 
		ka = row1+1;
	else 
		ka = col1+1;
	rvalue = svd( sd, row1, col1, u, v, EPS, ka ); /* svd decomposition */
	d = (double *)malloc( sizeof(double) * col1 );
	for ( i = 0; i < col1; i++ ) d[i] = 0;
	if ( row1 <= col1 )
	{
		for ( i = 0; i < row1; i++ )
			d[i] = sd[i*col1+i];
	}
	else
	{
		for ( i = 0; i < col1; i++ )
			d[i] = sd[i*col1+i];
	}

	/* calculate the variances: S = diag( S ); V = S .* S; */
	for ( x = 0; x < col1; x++ )
		V[x] = d[x] * d[x];
	/* get PC */
	for ( y = 0; y < col1; y++ )
		for ( x = 0; x < col1; x++ )
		{
			PC[y*col1+x] = v[x*col1+y];
		}

	/* project the original data: signals = PC' * data; */
	for ( y = 0; y < row; y++ )
		for ( x = 0; x < col; x++ )
		{
			signals[y*col+x] = 0;
			for ( i = 0; i < row; i++ )
				signals[y*col+x] += PC[i*row+y] * data1[i*col+x];
		}

	free( u );
	free( v );
	free( sd );
	free( d );
	free( data1 );

	return( 1 );
}

/*
   main function
*/
int main()
{
	double data[36*200] = { 0.1657 , 0.153734 , 0.142081 , 0.251764 , 0.377804 , 0.258022 , 0.114693 ,
		0.160907 , 0.174727 , 0.379836 , 0.159222, 0.0441471 ,0.0572539 ,0.0531755 ,
		0.109691 ,0.0702412 ,  0.12878 ,  0.21951, 0.0708107 , 0.150538 , 0.132882 ,
		0.21048 ,   0.2584 ,  0.19579 , 0.149602 ,0.0662608, 0.0720444, 0.0500046 ,
		0.0444386 ,0.0304997, 0.0534205 , 0.107401 , 0.107563 , 0.106128 ,0.0401567 ,
		0.0257917 , 0.126627 ,0.0542111 ,0.0583784 , 0.070953 , 0.173271 , 0.134469 ,
		0.0205659 ,0.0184031, 0.0646504 , 0.173839,  0.134699 , 0.141953,  0.123935 ,
		0.403615 , 0.338064, 0.0276296, 0.0195569 ,  0.11133, 0.00880792, 0.0100133 ,
		0.0885434 , 0.137133 , 0.282997,  0.136031, 0.0254509 ,0.00751139, 0.00465978 ,
		0.045524 ,0.0286458 , 0.185939 , 0.268302 , 0.403615  ,0.324212 ,0.0296117 ,
		0.0138353 , 0.017288 , 0.381674 , 0.133414, 0.0790253 ,0.0565466 , 0.183495 ,
		0.145601, 0.0182931, 0.0311231 , 0.170597 , 0.381674 , 0.173484 ,0.0229865 ,
		0.0130349, 0.0138423, 0.0341162,  0.011718, 0.0232139 , 0.307974,  0.230878 ,
		0.0566228, 0.0683206 , 0.106407 , 0.381674 , 0.148563 , 0.012566, 0.0254223 ,
		0.0656758 , 0.381674, 0.0916176 ,0.0220906, 0.0289558 , 0.068351, 0.0376567 ,
		0.00745692, 0.0191836 , 0.152901 , 0.407509 , 0.134734 ,0.00585401 ,0.00321429 ,
		0.00112352, 0.0166332, 0.0511024 ,0.00981551,  0.325877,  0.407509 , 0.155979 ,
		0.0387866, 0.00342308 ,0.00670642, 0.0232931, 0.0658894, 0.0574715,  0.407509 ,
		0.333105 , 0.065247, 0.0129569 ,0.00714528, 0.00874974 , 0.014964 ,0.0248751 ,
		0.00934936 ,  0.13611 , 0.340628, 0.0638973, 0.0247342, 0.00464725 ,0.0141411 ,
		0.0204052, 0.0359512, 0.0501483 , 0.199702 , 0.329576 , 0.149673 ,0.0621065 ,
		0.00780398 ,0.0329433, 0.0877676 , 0.103574,  0.135157 , 0.323298 , 0.329576 ,
		0.193228 ,0.0686342, 0.0385592, 0.0669586 , 0.155242,  0.129693 , 0.112795 ,
		0.224888 , 0.138876, 0.0651972 ,0.0293058 ,0.0109269 , 0.121833 , 0.159987 ,
		0.133214 , 0.161032 , 0.127249, 0.0744482 , 0.108871 ,0.0535444 ,0.0549595 ,
		0.237275 , 0.329576 , 0.265336 , 0.207154, 0.0535686 , 0.405886 ,  0.24916 ,
		0.0881366 ,0.0454212, 0.0714913, 0.0844257, 0.0698007 , 0.067213 , 0.347346 ,
		0.405886  , 0.13592, 0.0523571, 0.0275502 ,0.0450864, 0.0273337, 0.0287424 ,
		0.0380657 , 0.405886,  0.126306 , 0.136926 , 0.108214, 0.0857887 , 0.211221 ,
		0.216013 , 0.185343 , 0.110809 ,0.0516402 ,0.0770488 ,0.0496148, 0.0539304 ,
		0.0426737 , 0.096447, 0.0444659 , 0.044463, 0.0325991 ,0.0508148,  0.427374 ,
		0.130947, 0.0249754 , 0.023727, 0.0638681, 0.0653918, 0.0569512 ,0.0480317, 
		0.392653,  0.427374 , 0.288818 ,0.0484743, 0.0226702, 0.0942696 ,0.0977035 ,
		0.14445, 0.0623437 , 0.282399 ,0.0568553 ,0.0127045 ,0.00721119, 0.0168199 ,
		0.0400718, 0.0513284 , 0.091616 ,0.0470645, 0.0589789 ,0.0814543, 0.0306087 ,
		0.00732747 , 0.022916 ,0.0489248 , 0.144673 , 0.368764 ,0.0725378, 0.0448737 };
	int row = 36, col = 7;
	double signals[36*200], PC[1296], V[36];
	int x, y;

	pca2( data, row, col, signals, PC, V );

	double sum = 0, subsum = 0;
	int l;
	for ( y = 0; y < row; y++ )
		sum += V[y];
	for ( y = 0; y < row; y++ ){
		subsum += V[y];
		if(subsum/sum >= 0.98){
			l = y;
			break;
		}
	}

	printf( "l: %d ", l  );
	double subsignals[36*200];
	printf( "Project to Principal Component: \n" );
	for ( y = 0; y < l; y++ )
	{
		for ( x = 0; x < col; x++ )
		{ 
			subsignals[y*col+x] = signals[y*col+x];
			printf( "%7.4f ", subsignals[y*col+x] );
		}
		printf( "\n" );
	}
	printf( "Eigen Values (in deceasing order): \n" );
	for ( y = 0; y < row; y++ )
		printf( "%7.4f ", V[y] );
	printf( "\n" );
	printf( "Eigen Vectors: \n" );
	for ( y = 0; y < row; y++ )
	{
		for ( x = 0; x < row; x++ )
		{
			printf( "%7.4f ", PC[y*row+x] );
		}
		printf( "\n" );
	}

	return( 1 );
}