sparsesvd_partialeig.c

关于有直接稀疏PCA的方法

/* 
Finds a sparse rank-one approximation to a given symmetric matrix A, by solving the SDP

						   min_X lambda_max(A+X) : X = X', abs(X(i,j)) <= rho, 1<=i,j<= n
and its dual:

						   max_U 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			solves the primal SDP 
U			dual variable, solves the dual SDP 
u			largest eigenvector of U 
F			Average gradient

This code implements Nesterov's smooth minimization algorithm. See: Y. Nesterov "Smooth Minimization of NonSmooth Functions."
Here, the gradient is only only computed aproximately. See A. d'Aspremont "Smooth optimization with approximate gradient."

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

#include "sparsesvd.h"

void sparse_rank_one_partialeig(double *Amat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iter, 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);
	int work_size=3*n+n*n;
	// Working variables
	double d1,sig1,d2,sig2,norma12,mu,Ntheo,L;
	double alpha,gapk,dmax=0.0,fmu,lambda;
	int n2=n*n,incx=1,precision_flag=0,iteration_flag=0;
	int lwork=work_size,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]="V",uplo[1]="U";
	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 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_store=(double *) calloc(n*n,sizeof(double));
	double *evector_temp=(double *) calloc(n*n,sizeof(double)); // TODO: one of these matrices can go.
	double *eig=(double *) calloc(n,sizeof(double));
	double eigcut,tolweight=.75,tol=.01,numeigstemp=numeigs;
	int *count=(int *) calloc(n,sizeof(int)); // Records the distribution of eig computations
	int checkgap_count=0,firstiter=0; // added for test variables
	int arcount;

	// Start...
	if (info>=1){mexPrintf("DSPCA starting... Sparse eig. maximization using approximate gradient.\n");mexEvalString("drawnow;");}
	// Test malloc results
	if ((Fmat==NULL) || (Vmat==NULL) || (bufmata==NULL) || (bufmatb==NULL)|| (Dvec==NULL) || (workvec==NULL) || (gvec==NULL) || (hvec==NULL) ||(evector_temp==NULL)||(evector_store==NULL)||(iwork==NULL)||(workvec2==NULL)||(numeigs_matlab==NULL)||(eig==NULL)){
		mexPrintf("DSPCA: memory allocation failed ... \n");mexEvalString("drawnow;");return;}
	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

	// First, compute some local params
	norma12=1.0;d1=rho*rho*n*n/2.0;sig1=1.0;d2=log(n);sig2=0.5;mu=tol/(2.0*d2); // TODO: can we get a less conservative d1?
	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)==0){		
		// eigenvalue decomposition of A+X 
		cblas_dcopy(n2,Xmat,incx,Vmat,incx);
		alpha=1.0; 
		cblas_daxpy(n2,alpha,Amat,incx,Vmat,incx);
		cblas_dcopy(n2,Vmat,incx,bufmata,incx);		
		// do partial eigenvalue approximation to exp(A+X)
		*numeigs_matlab=1.0*numeigs;*count=0;arcount=0;
		fmu=partial_eig(n,k,mu,eigcut,bufmata,bufmatb,numeigs_matlab,evector_temp,evector_store,eig,Dvec,gvec,hvec,Vmat,Umat,workvec,count,addeigs,perceigs,check_for_more_eigs,&arcount); 
		numeigs=(int)(*numeigs_matlab);dmax=bufmata[0];
		// 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.0/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.0);
		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)||(k%Nperiod==0)||(((double)(clock())/CLOCKS_PER_SEC-last_time)>=900)){
			gapk=dmax-doubdot(Amat,Umat,n2)+rho*doubasum(Umat,n2);
			if (firstiter==1) {dualitygap_alliter[checkgap_count]=gapk;cputime_alliter[checkgap_count]=cputime;perceigs_alliter[checkgap_count]=100.0*numeigs/n;checkgap_count++;}
			if (firstiter==0){// If first iteration, reset precision targets
				tol=gapk*gapchange;norma12=1.0;d1=rho*rho*n*n/2.0;sig1=1.0;d2=log(n);sig2=0.5;mu=tol/(2.0*d2);
				L=(d2*norma12*norma12)/(2.0*sig2*tol);eigcut=(1-tolweight)*(tol/10)/(rho*n);
				alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);cblas_dscal(n2,alpha,Fmat,incx);numeigs=numeigstemp;}		
			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)&&(firstiter==1))||(precision_flag+iteration_flag>0)){
				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;");}
			if (firstiter==0) {firstiter=1;k--;}}
		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);
	dsyev(jobz,uplo,&n,Vmat,&n,Dvec,workvec,&lwork,&inflapack); // TODO: switch to ARPACK here...
	indmax=idxmax(Dvec,n);dmax=Dvec[indmax]; 
	for (i=0;i<n;i++) {uvec[i]=Vmat[(indmax)*n+i];}
	*iter=k; // return total number of iterations
	// Free everything
	free(Vmat);
	free(bufmata);
	free(bufmatb);
	free(Dvec);
	free(workvec);
	free(workvec2);
	free(gvec);
	free(hvec);
	free(iwork);
	free(numeigs_matlab);
	free(evector_store);
	free(evector_temp);
	free(eig);
	free(count);
}