sparsesvd_partialeig_matlab.c

关于有直接稀疏PCA的方法

/* Main function
sparse_rank_one(double *Amat, int n, double rho, double tol, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, int WarmStart, int info) 
SPARSERANKONE finds a sparse rank-one approximation to 
a given symmetric matrix A, by solving the SDP
min_U lambda_max(A+X) : X = X', abs(X(i,j)) <= rho, 1<=i,j<= n
and its dual:
max_X Tr(UA) - rho sum_ij |U_ij| : U=U', U \succeq 0, Tr(U)=1

***	inputs: ***
A			nxn symmetric matrix (left unchanged)
n			problem size
rho			non-negative scalar 
gapchange	required change in gap from first gap (default: 1e-4) 
MaxIter		maximum number of iterations
info		controls verbosity: 0 silent, n>0 frequency of progress report
WarmStart	0 if cold start, k0 if WarmStart (total number of iterations in previous run)
F			Average gradient (for warm start, Fmat is updated)

***	outputs: ***
X			symmetric matrix that solves the above SDP 
U			dual variable, solves the dual SDP 
u			largest eigenvector of U 
F			Average gradient
k			number of iterations run

This code implements Nesterov's smooth minimization algorithm. 
See: Y. Nesterov "Smooth Minimization of NonSmooth Functions", CORE DP 2003/12. 

Last Modified: A. d'Aspremont, Laurent El Ghaoui, Ronny Luss July 2006.
http://www.carva.org/alexandre.daspremont
*/

#include "sparsesvd.h"

void sparse_rank_one_partialeig_matlab(double *Amat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iters, int info, int numeigs, int addeigs, int checkgap, double perceigs, int check_for_more_eigs, double *dualitygap_alliter, double *cputime_alliter, double *perceigs_alliter)
{
	// Hard parameters
	int Nperiod=imaxf(1,info),dspca_finished=0;
	int work_size=3*n+n*n,changedmu=0;
	// Working variables
	double d1,sig1,d2,sig2,norma12,mu,Ntheo,L;
	double alpha,gapk;
	double dmax=0.0,fmu,lambda;
	int n2=n*n,incx=1,precision_flag=0,iteration_flag=0,error_flag=0;
	int lwork,inflapack,indmax,k=0,i,j;
	double cputime,last_time=(double)clock();double start_time=(double)clock();int left_h=0,left_m=0,left_s=0;
	char jobz[1],uplo[1];
	double *Vmat=(double *) calloc(n*n,sizeof(double));
	double *bufmata=(double *) calloc(n*n,sizeof(double));
	double *bufmatb=(double *) calloc(n*n,sizeof(double));
	double *workvec=(double *) calloc(work_size,sizeof(double));
	double *hvec=(double *) calloc(n,sizeof(double));
	int work_size3=8*n;
	double *workvec2=(double *) calloc(work_size3,sizeof(double));
	int *iwork=(int *) calloc(5*n,sizeof(int));
	double *numeigs_matlab=(double *) calloc(1,sizeof(double));
	double *evector_temp=(double *) calloc(n*n,sizeof(double));
	double *evalue=(double *) calloc(n*n,sizeof(double));
	double *evector_store=(double *) calloc(n*n,sizeof(double));
	double eigcut,tolweight=.75,tol=.01;
	int *count=(int *) calloc(n,sizeof(int));
	mxArray *output[3],*input[4];
	double *Fmattemp=(double *) calloc(n*n,sizeof(double));
	double *Xmattemp=(double *) calloc(n*n,sizeof(double));
	int checkgap_count=0; // added for test variables

	// Start...
	if (info>=1)
	{
		mexPrintf("DSPCA starting ... \n");
		mexEvalString("drawnow;");
	}
	// Test malloc results
	if ((Fmat==NULL) || (Vmat==NULL) || (bufmata==NULL) || (bufmatb==NULL)|| (workvec==NULL) || (hvec==NULL) ||(evector_temp==NULL)||(evector_store==NULL)||(evalue==NULL)||(iwork==NULL)||(workvec2==NULL)||(numeigs_matlab==NULL)||(Fmattemp==NULL)||(Xmattemp==NULL))
	{
		mexPrintf("DSPCA: memory allocation failed ... \n");
		mexEvalString("drawnow;");return;
	}
	eigcut=(1-tolweight)*(tol/10)/(rho*pow(n,1.5)); // scale delta (tol/10) to get eig threshold
	tol=tolweight*tol; // scale of .5 for partial eig precision
	mexEvalString("options.disp=0\;");           // for use in calling Matlab function eigs	
	mexEvalString("options.maxit=500\;");		 // for use in calling Matlab function eigs	
	input[0] = mxCreateDoubleMatrix(n,n,mxREAL); // for use in calling Matlab function eigs
	input[1] = mxCreateDoubleMatrix(1,1,mxREAL);
	input[2]=mxCreateString("la");
	input[3]=mexGetVariable("caller","options");

	// First, compute some local params
	d1=rho*rho*n*n/2.0;sig1=1.0;d2=log(n);sig2=0.5;norma12=1.0;mu=tol/(2.0*d2);
	Ntheo=(4.0*norma12*sqrt(d1*d2/(sig1*sig2)))/tol;Ntheo=ceil(Ntheo);
	L=(d2*norma12*norma12)/(2.0*sig2*tol);
	alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);
	
	cputime=start_time;

	while ((precision_flag+iteration_flag+error_flag)==0)
	{
		if (k==1 && changedmu==0) {  // after 1st iteration and when algorithm hasn't been restarted, adjust tol to be a percentage change in original gap
			gapk=dmax-doubdot(Amat,Umat,n2)+rho*doubasum(Umat,n2);
			tol=gapchange*gapk;
			eigcut=(1-tolweight)*(tol/10)/(rho*n); // scale delta (tol/10) to get eig threshold
			tol=tolweight*tol; // scale of .5 for partial eig precision
			mu=tol/(2.0*d2);
			L=(d2*norma12*norma12)/(2.0*sig2*tol);
			alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);
			alpha=0.0;cblas_dscal(n2,alpha,Fmat,incx);
			k=0;
			checkgap_count=0;			
			free(Fmattemp);
			free(Xmattemp);
			changedmu=1;
		}
		*count=0;
		// eigenvalue decomposition of A+X 
		cblas_dcopy(n2,Xmat,incx,Vmat,incx);
		alpha=1.0; 
		cblas_daxpy(n2,alpha,Amat,incx,Vmat,incx);
		symmetrize(Vmat,bufmata,n);	// symmetrize A+X so no precision problems		
		// do partial eigenvalue approximation to exp(A+X)
		cblas_dcopy(n2,bufmata,incx,Vmat,incx);		
		*numeigs_matlab=1.0*numeigs;
		fmu=partial_eig_matlab(n,k,mu,eigcut,bufmata,bufmatb,numeigs_matlab,
			evector_temp,evector_store,evalue,input,output,
			hvec,Vmat,Umat,workvec,count,dmax,addeigs,perceigs,check_for_more_eigs);
		numeigs=(int)(*numeigs_matlab);
		dmax=bufmata[0];
		error_flag=(int)bufmata[1];
		// update gradient's weighted average 
		alpha=((double)(k)+1)/2.0;
		cblas_daxpy(n2,alpha,Umat,incx,Fmat,incx);
		// find a projection of X-Gmu/L on feasible set 
		cblas_dcopy(n2,Xmat,incx,bufmata,incx);
		alpha=-1/L;
		cblas_daxpy(n2,alpha,Umat,incx,bufmata,incx);
		// project again
		alpha=-(sig1/L);
		cblas_dcopy(n2,Fmat,incx,bufmatb,incx);
		cblas_dscal(n2,alpha,bufmatb,incx);
		// update X
		lambda=2.0/((double)(k)+3);
		for (j=0;j<n;j++){
			for (i=0;i<n;i++){
				Xmat[j*n+i]=lambda*dsignf(bufmatb[j*n+i])*dminif(rho,dabsf(bufmatb[j*n+i]))+(1-lambda)*dsignf(bufmata[j*n+i])*dminif(rho,dabsf(bufmata[j*n+i]));}}
		// check convergence and gap periodically
		cputime=((double)clock()-start_time)/CLOCKS_PER_SEC;
		if (k%checkgap==0) {  // check gap more often than printing info
			gapk=dmax-doubdot(Amat,Umat,n2)+rho*doubasum(Umat,n2);
			dualitygap_alliter[checkgap_count]=gapk;
			cputime_alliter[checkgap_count]=cputime;
			perceigs_alliter[checkgap_count]=100.0*numeigs/n;
			checkgap_count++;
			if (gapk<=tol) dspca_finished=1;
		}
		if ((changedmu==1)&&((dspca_finished==1)||(k%Nperiod==0)||(((double)(clock())/CLOCKS_PER_SEC-last_time)>=900)))
		{
			gapk=dmax-doubdot(Amat,Umat,n2)+rho*doubasum(Umat,n2);
			last_time=(double)(clock())/CLOCKS_PER_SEC;
			if (gapk<=tol) precision_flag=1;
			if (k>=MaxIter) iteration_flag=1;
			// report iteration, gap and time left
			if (info>=1){
				left_h=(int)floor(cputime/3600);left_m=(int)floor(cputime/60-left_h*60);left_s=(int)floor(cputime-left_h*3600-left_m*60);
				mexPrintf("Iter: %.3e   Obj: %.4e    Gap: %.4e   CPU Time: %2dh %2dm %2ds	 %% Eigs Used: %.2f\n",(double)(k),dmax,gapk,left_h,left_m,left_s,100.0*numeigs/n);
				mexEvalString("drawnow;");
			}
		}
		k++;
	}
	// set dual variable and output vector
	// eigenvalue decomposition of A+X 
	alpha=0.0;
	cblas_dscal(n2,alpha,Vmat,incx);
	cblas_dcopy(n2,Umat,incx,Vmat,incx);
	*jobz='V';*uplo='U';lwork=work_size;
	dsyev(jobz,uplo,&n,Vmat,&n,hvec,workvec,&lwork,&inflapack);
	indmax=idxmax(hvec,n);dmax=hvec[indmax]; 
	for (i=0;i<n;i++) {uvec[i]=Vmat[(indmax)*n+i];}
	// Return total number of iterations
	*iters=k;
	// Free everything
	free(Vmat);
	free(bufmata);
	free(bufmatb);
	free(workvec);
	free(workvec2);
	free(hvec);
	free(iwork);
	free(numeigs_matlab);
	free(evector_store);
	free(evector_temp);
	free(evalue);
	free(count);
}