sparsesvd.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(double *Amat, int n, double rho, double gapchange, int MaxIter, double *Xmat, double *Umat, double *uvec, double *Fmat, double *iter, int info, int checkgap, double *dualitygap_alliter, double *cputime_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,beta,buf,gapk;
	double dmax=0.0,fmu,lambda,tol=.01;
	int n2=n*n,incx=1,precision_flag=0,iteration_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],transa[1],transb[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 *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 checkgap_count=0,firstiter=0; // added for test variables
	
	// Start...
	if (info>=1){
		mexPrintf("DSPCA starting... Sparse eig. maximization.\n");
		mexEvalString("drawnow;");}
	// Test malloc results
	if ((Fmat==NULL) || (Vmat==NULL) || (bufmata==NULL) || (bufmatb==NULL) || (Dvec==NULL) || (workvec==NULL) || (gvec==NULL) || (hvec==NULL)){
		mexPrintf("DSPCA: memory allocation failed ... \n");
		mexEvalString("drawnow;");return;}
	// First, compute some local params
	d1=rho*rho*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)==0){ 
		// eigenvalue decomposition of A+X 
		cblas_dcopy(n2,Xmat,incx,Vmat,incx);
		alpha=1.0;cblas_daxpy(n2,alpha,Amat,incx,Vmat,incx);
		*jobz='V';*uplo='U';lwork=work_size;
		dsyev(jobz,uplo,&n,Vmat,&n,Dvec,workvec,&lwork,&inflapack); // call LAPACK (most CPU time is here) 
		// compute fmu(X) = mu*log(trace((exp(A+X)/mu)))-mu*log(n) reliably 
		indmax=idxmax(Dvec,n);dmax=Dvec[indmax];
		for (i=0;i<n;i++) {hvec[i]=exp((Dvec[i]-dmax)/mu);}
		buf=doubsum(hvec,n);fmu=dmax+mu*log(buf/n);
		// compute gradient of fmu w.r.t. X, which is the dual variable U 
		alpha=0.0;cblas_dscal(n2,alpha,bufmatb,incx);
		for (i=0;i<n;i++) {gvec[i]=hvec[i]/buf;bufmatb[i*n+i]=gvec[i];}
		alpha=1.0;beta=0.0;*transa='N';*transb='T';
		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,Vmat,n,bufmatb,n,beta,bufmata,n);
		cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,n,n,n,alpha,bufmata,n,Vmat,n,beta,Umat,n);
		// 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)||(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;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);
				alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);cblas_dscal(n2,alpha,Fmat,incx);}		
			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\n",(double)(k),dmax,gapk,left_h,left_m,left_s);
				mexEvalString("drawnow;");}
			if (firstiter==0) {firstiter=1;k--;}}
		k++;} // End of main loop...
	// 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];}
	*iter=k; // return total number of iterations
	// Free everything
	free(Vmat);
	free(bufmata);
	free(bufmatb);
	free(Dvec);
	free(workvec);
	free(gvec);
	free(hvec);
}