sparsesvd_pade.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_pade(double *Amat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iters, int info,  int checkgap, double *dualitygap_alliter, double *cputime_alliter)
{
	// Hard parameters
	int Nperiod=imaxf(1,info),changedmu=0;
	int work_size=3*n+n*n;
	double d1,sig1,d2,sig2,norma12,mu,Ntheo,L;
	double alpha,gapk,first_gapk0;
	double dmax=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],tolerancestr[100];
	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 *Dvec=(double *) calloc(n,sizeof(double));
	double *workvec=(double *) calloc(work_size,sizeof(double));
	double *gvec=(double *) calloc(n,sizeof(double));
	double *hvec=(double *) calloc(n,sizeof(double));
	int ideg=6,work_size2=4*n*n+ideg+1;
	int *work_in=(int *) calloc(n,sizeof(int));
	double *work_out=(double *) calloc(work_size2,sizeof(double));
	double trace,bufmata_shift=0.0,tol=.01;
	int work_size3=8*n;
	double *workvec2=(double *) calloc(work_size3,sizeof(double));
	double *numeigs_matlab=(double *) calloc(1,sizeof(double));
	double *evalue=(double *) calloc(1,sizeof(double));
	int *iwork=(int *) calloc(5*n,sizeof(int));
	mxArray *input[1],*output[1];
	double *Fmattemp=(double *) calloc(n*n,sizeof(double));
	double *Xmattemp=(double *) calloc(n*n,sizeof(double));
	int checkgap_count=0; // added for test variables
	double tolerance;
	char which[2]="LA"; // Arpack: we want largest algebraic eigs...
	int ncv=2,nconv,nummatvec,info_arpack,maxitr_arpack=1000;
	double *evector_temp=(double *) calloc(n,sizeof(double)); // only one eigenvector
	mxArray *input_eigs[4],*output_eigs[3];
	
	// Start...
	if (info>=1)
	{
		mexPrintf("DSPCA starting ... \n");
		mexEvalString("drawnow;");
	}
	// Test malloc results
	if ((Fmat==NULL) || (Vmat==NULL) || (bufmata==NULL) || (bufmatb==NULL) || (Dvec==NULL) || (workvec==NULL) || (gvec==NULL) || (hvec==NULL)||(work_in==NULL)||(work_out==NULL)||(workvec2==NULL)||(numeigs_matlab==NULL)||(evalue==NULL)||(iwork==NULL)||(Fmattemp==NULL)||(Xmattemp==NULL))
	{
		mexPrintf("DSPCA: memory allocation failed ... \n");
		mexEvalString("drawnow;");return;
	}
	input[0] = mxCreateDoubleMatrix(n,n,mxREAL); // for use in calling Matlab function expm
	mexEvalString("options.disp=0\;");			 // for use in calling Matlab function eigs
	mexEvalString("options.maxit=500\;");		 // for use in calling Matlab function eigs
	input_eigs[0] = mxCreateDoubleMatrix(n,n,mxREAL); 
	input_eigs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
	*numeigs_matlab=1.0;
	memcpy(mxGetPr(input_eigs[1]),numeigs_matlab,sizeof(double));		
	input_eigs[2]=mxCreateString("la");
	input_eigs[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;
			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;
			changedmu=1;
		}
		// eigenvalue decomposition of A+X 
		cblas_dcopy(n2,Xmat,incx,Vmat,incx);
		alpha=1.0; 
		cblas_daxpy(n2,alpha,Amat,incx,Vmat,incx);
		// do pade approximation to exp(A+X)
		symmetrize(Vmat,bufmata,n);	// symmetrize A+X so no precision problems		
		cblas_dcopy(n2,bufmata,incx,Vmat,incx);
		bufmata_shift=frobnorm(bufmata,n);// simple bound on largest magnitude eigenvalue
		//i=1;tolerance=.0001/bufmata_shift;
		//info_arpack=simarpack(bufmata,n,i,ncv,tolerance,which,maxitr_arpack,0,evalue,evector_temp,&nconv,&nummatvec);
		//if (info_arpack<0) {
			sprintf(tolerancestr,"options.tol=%.15f;",.0001/bufmata_shift); // this will be the tolerance parameter for eigs
			mexEvalString(tolerancestr);
			memcpy(mxGetPr(input_eigs[0]),bufmata,n*n*sizeof(double));
			mexCallMATLAB(3,output_eigs,4,input_eigs,"eigs");
			memcpy(evalue,mxGetPr(output_eigs[1]),sizeof(double));
		//}
		dmax=evalue[0];
		for (i=0;i<n;i++) {bufmata[i*n+i]-=dmax;}
		alpha=1.0/mu;cblas_dscal(n2,alpha,bufmata,incx);
		memcpy(mxGetPr(input[0]),bufmata,n*n*sizeof(double));
		mexCallMATLAB(1,output,1,input,"expm");
		memcpy(Umat,mxGetPr(output[0]),n*n*sizeof(double));
		mxDestroyArray(output[0]);
		trace=0.0;for (i=0;i<n;i++){trace+=Umat[i*n+i];}
		fmu=dmax+mu*log(trace)-mu*log(n);
		alpha=1.0/trace;cblas_dscal(n2,alpha,Umat,incx);
		// 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 ((changedmu==1)&&((k%checkgap==0)||(k%Nperiod==0)||(((double)(clock())/CLOCKS_PER_SEC-last_time)>=900))){
			gapk=dmax-doubdot(Amat,Umat,n2)+rho*doubasum(Umat,n2);
			dualitygap_alliter[checkgap_count]=gapk;cputime_alliter[checkgap_count]=cputime;checkgap_count++;
			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)&&(k%Nperiod==0))||(precision_flag==1)||iteration_flag==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\n",(double)(k),dmax,gapk,left_h,left_m,left_s);
				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,Dvec,workvec,&lwork,&inflapack);
	indmax=idxmax(Dvec,n);dmax=Dvec[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(Dvec);
	free(workvec);
	free(workvec2);
	free(gvec);
	free(hvec);
	free(iwork);
	free(numeigs_matlab);
	free(evalue);
	free(evector_temp);
}