sparsesvd_padedeprecated.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, int *WarmStart, int info,  int checkgap, double *dualitygap_alliter, double *cputime_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,fmu,lambda;
	int n2=n*n,incx=1,precision_flag=0,iteration_flag=0,error_flag=0;
	int lwork,inflapack,indmax,k=*WarmStart,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],tolerance[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,eigsflag;
	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],*input_eigs[4],*output_eigs[3];
	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) || (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
	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);
	if (*WarmStart==0)
	{
		alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);
	} else {
		cblas_dcopy(n2,Fmat,incx,Fmattemp,incx);
		cblas_dcopy(n2,Xmat,incx,Xmattemp,incx);
	}
	MaxIter+=*WarmStart-1;
	cputime=start_time;
	while ((precision_flag+iteration_flag+error_flag)==0)
	{
		if (k-(*WarmStart)==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);
			if (*WarmStart==0)
			{
				alpha=0.0;cblas_dscal(n2,alpha,Xmat,incx);
				alpha=0.0;cblas_dscal(n2,alpha,Fmat,incx);
			} else {
				cblas_dcopy(n2,Fmattemp,incx,Fmat,incx);
				cblas_dcopy(n2,Xmattemp,incx,Xmat,incx);
			}				
			k=*WarmStart;
			checkgap_count=0;
			free(Fmattemp);
			free(Xmattemp);
			changedmu=1;
		}

		// eigenvalue decomposition of A+X 
		cblas_dcopy(n2,Xmat,incx,Vmat,incx);
		alpha=1.0;  
		cblas_daxpy(n2,alpha,Amat,incx,bufmata,incx);
		// do pade approximation to exp(A+X)
		cblas_dcopy(n2,bufmata,incx,Vmat,incx);
		bufmata_shift=frobnorm(bufmata,n);// simple bound on largest magnitude eigenvalue
		sprintf(tolerance,"options.tol=%.15f;",.0001/bufmata_shift); // this will be the tolerance parameter for eigs
		mexEvalString(tolerance);
		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];
		memcpy(evalue,mxGetPr(output_eigs[2]),sizeof(double));
		eigsflag=evalue[0];
		if (eigsflag>0.50) {
			mexPrintf("ERROR: Not all eigenvalues converged in calling Matlab eigs function in pade approximation.  Perhaps maxit not set high enough\n");
			error_flag=1;
		}
		mxDestroyArray(output_eigs[0]);mxDestroyArray(output_eigs[1]);
		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 (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;
			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\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
	*WarmStart=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);
}