logreglasso.cpp

统计学习软件包

#include "LogregLASSO.h"
#include <math.h>
#include <fstream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iomanip.h>
#include <time.h>
#include "my_mat.h"


LogregLASSO::LogregLASSO(const My_Matrix* AA, const int* Yv, double c_val, 
			 double constraint, int* i_index, double* i_beta)
{

  int i,j,k;
  const_val = c_val;

  A = AA;
  n = A->Nrows();
  Y = new double[n];
  Projected = new double[n];
  rv = new double[1]; 
  for(i=0; i<n; i++)
      Y[i] = Yv[i];

  kc = Acols = A->Ncols();

  alpha = new double[kc];
  old_index = new int[kc];
  index = new int[kc];
  init_index = new int[kc];
  init_beta = new double[kc];

  for(i=0; i<kc; i++){
      init_index[i] = i_index[i];
      init_beta[i] = i_beta[i];
  }
  mu = constraint;

}
LogregLASSO::~LogregLASSO()
{
    delete [] Y;
    delete [] alpha; 
    delete [] Projected;
    delete [] init_index;
    delete [] init_beta;
    delete [] rv;  
    delete [] index;
    delete [] old_index;
}
 
void LogregLASSO::SetInitVals(int* iind, double* ibeta, double constraint){
    int i;
    for(i=0; i<kc; i++){
	 init_index[i] = iind[i];
	 init_beta[i] = ibeta[i];
    }
    mu = constraint;
}


void LogregLASSO::init_LASSO( double* h, double* r, double* v, double* beta,  double* betaDagger, double* Theta, int* c){
    int i;
    for(i=0; i<kc; i++){
	beta[i] = 0.0;
	betaDagger[i] = 0.0;	
	h[i] = 0.0;  
	v[i] = 0.0;
    }
    for(i=0; i<n; i++){
	r[i] = 0.0;
    }
    (*c) =0;
    for(i=0; i<kc; i++){
	old_index[i] = init_index[i];
	index[i] = init_index[i];
	beta[i] = init_beta[i];
	if(index[i] == 1){
	    if(beta[i]/ fabs(beta[i] ) <0.0)
		Theta[i] = -1;
	    else
		Theta[i] = 1;

	    (*c) +=1;
	}
	else{
	    Theta[i] = -1.;
	    beta[i] =0.0;
	}
    }

}

void LogregLASSO::SUBlasso(double* Solvec, double* YTd)
{
  int  i,j;
  
      
  double* beta = new double[kc];
  double* betaDagger = new double[kc];
  double* Theta = new double[kc];
  double* h = new double[kc];
  double* r = new double[n];
  double* v = new double[kc];
 
  
  int cn;
  int* c = &cn;
  init_LASSO(h, r, v, beta,  betaDagger,  Theta, c);
  
  proceed = 1;
  
  SUBlassoOPT_rob(YTd, h, beta, Theta,   c);
  SUBlassoCheckSignFeasible(YTd,h, beta, Theta,  c, betaDagger );
  SUBlassoCheckViolation( YTd,h, beta, Theta,  c, betaDagger,r,v);
  
  while(!proceed){
      cout << " >>>>>>>>>>>>>>> RESTART <<<<<<<<<<<<<<<<<<<<<<<<<" <<endl;
      proceed = 1;
      // restart with c = 1;
       for(i=0; i<kc; i++){
	   init_index[i] = 0;
	   init_beta[i] = 0.;
       }
       int ri; 
       ri =  rand()%(A->Ncols()-1);
       init_index[ri] = 1;
       init_beta[ri] = 1e-7;


       init_LASSO(h, r, v, beta,  betaDagger,  Theta, c);
      
       SUBlassoOPT_rob(YTd, h, beta, Theta,   c); 
       SUBlassoCheckSignFeasible(YTd,h, beta, Theta,  c, betaDagger );
       SUBlassoCheckViolation( YTd,h, beta, Theta, c, betaDagger,r,v);
  }
  
  
  for(i=0; i<kc; i++){
    Solvec[i] = betaDagger[i];
  }

  rv_numb = *c;
  
  delete [] rv;
  rv = new double[*c]; 
  j = 0;
  for(i=0; i<kc; i++){
    if(index[i] == 1){
      rv[j] = i;
      j++;
    }
  }

  
  delete [] beta;
  delete [] betaDagger;
  delete [] Theta;
  delete [] h;
  delete [] r;
  delete [] v;
 

 
}
void LogregLASSO::multA_Beta_index(double* beta, double* r){
    int i,j;
    double sum;
    for(i=0; i<n; i++){
	sum = 0.;
	for(j=0; j<kc; j++){
	    if(index[j] == 1){
		sum += A->el(i,j)*beta[j];
	    }
     
	}
	r[i] = sum;
     }
    
}
void LogregLASSO::SUBlassoCheckViolation(double* YTd, double* h, double* beta, 
				  double* Theta, 
				  int *c, double* betaDagger, double* r, 
				  double* v )
{
  int i,j;
  double sum;
  // cout << "CheckViolation ..." << endl;

  multA_Beta_index(betaDagger,r);
  for(i=0; i<n; i++)
      r[i] = -(1.0/(1.0+exp(-r[i])) - YTd[i]);
      


    for(j=0; j<kc; j++){
	sum = 0.;
	for(i=0; i<n; i++)
	    sum += A->el(i,j) * r[i];
	v[j] = sum;
    }



  double max = -1.;
  int indMax; 

  for(i=0; i<kc; i++){   
    if(index[i] == 0)
	continue;
    if(fabs( v[i]) > max){
	max = fabs( v[i]);
    }
  }
  
  for(i=0; i<kc; i++){
    v[i] /= max;
  }

  max = -1;
 
  for(i=0; i<kc; i++){
      if(index[i] == 1){
	  continue;
      }
      if(fabs(v[i]) > max){
	  max = fabs(v[i]);
	  indMax = i;  
      }
  }

 double err_tol = 1.0 + 1e-8;
/*
 cout << "max viol: ";

 if(max-1.0 > 0.0)
     cout << max-1.0 <<endl;
 else
     cout << 0 << endl;
*/
 if(max >err_tol){
   // ----------- B2 -----
   
   index[indMax] = 1;
   (*c) +=1;

   betaDagger[indMax]= 0.0;
   Theta[indMax] = v[indMax] /fabs(v[indMax]);
   
   for(i=0; i<kc; i++){
       beta[i] = betaDagger[i];       
   }
   
   if(proceed){
       SUBlassoOPT_rob(YTd, h, beta, Theta,   c);
       SUBlassoCheckSignFeasible(YTd,h, beta, Theta,   c, betaDagger );
       
       for(i=0; i<kc; i++){    
	   betaDagger[i] = beta[i] +h[i];
       }
       for(i=0; i<kc; i++) 
	   Theta[i] = betaDagger[i] / fabs( betaDagger[i]);
       
       SUBlassoCheckViolation(YTd,h, beta, Theta,  c, betaDagger,r,v );
   }
   
 }
}
void LogregLASSO::SUBlassoCheckSignFeasible(double* YTd,double* h, double* beta, 
				     double* Theta,   
				     int* c, double* betaDagger )
{
    int isFeas = 1;
    int indK,i,j;
    double sum;
    
    for(i=0; i<kc; i++){    
	betaDagger[i] = beta[i] +h[i];
    }
 
    for(i=0; i<kc; i++){
	if(index[i] == 1){
	    if( (betaDagger[i] / fabs( betaDagger[i])) <0 && Theta[i] >0){
		isFeas = 0;
		break;
	    }
	    if( (betaDagger[i] / fabs( betaDagger[i])) >0 && Theta[i] <0){
		isFeas = 0;
		break;
	    } 
	}
    }
 
 
  if(isFeas == 0){
  
    // -----------SUBlassoA1 -------------
    double gamma = 2.;   
    indK = 0;  
    for(i=0; i<kc; i++){
      if(index[i] == 0)
	continue;
 
      if(-beta[i]/h[i] > 0. && -beta[i]/h[i] < gamma){
	  gamma = -beta[i]/h[i];
	  indK = i;
      }
    }
   
    
    if(gamma >1.){
	cout<<"***** roundoff errors in gamma ********* :" << gamma<< endl;
	gamma = 1e-10;
	indK = index[0];
	proceed = 0;
	isFeas = 0;	
    }
    
    
    for(i=0; i<kc; i++){
	beta[i] += gamma*h[i];
    }

    beta[indK] = 0.0;
    
    if(proceed){    
	// ------- SUBlassoA2 ----------------
	Theta[indK] = -Theta[indK];
	SUBlassoOPT_rob(YTd, h, beta, Theta,  c);
	// ------- check A2-1      ----------------
	int isFeas = 1;
	for(i=0; i<kc; i++)
	    betaDagger[i] = beta[i] +h[i];
	
	for(i=0; i<kc; i++){
	    
	    if(index[i] == 1){
		if( (betaDagger[i] / fabs( betaDagger[i])) <0  && Theta[i] ==1){
		    isFeas = 0;
		    break;
		}
		if( (betaDagger[i] / fabs( betaDagger[i])) >0 && Theta[i] ==-1){
		    isFeas = 0;
		    break;
		}
		
	    }
	}
	if(isFeas == 0) { 
	    // ------------A2-2 ------
	    
	    index[indK] = 0;
	    betaDagger[indK] = 0.0;
	    h[indK] = 0.0;
	    Theta[indK] = -1;
	    (*c) -= 1;
	    // cout << "# RVs: " <<*c<< endl;
	    
	    SUBlassoOPT_rob(YTd, h, beta, Theta,   c);
	    SUBlassoCheckSignFeasible(YTd,h, beta, Theta,   c, betaDagger );
	}
	
    }
  }  
}


void LogregLASSO::SUBlassoOPT_rob(double* YTd, double* h, double* beta, double* Theta,   int *cc){

  int i,ii,j,k,l,iter;
  int c = (*cc);
  double sum, res;
  int start = 0;
  int kk = 0;

  ii = 0;
  My_Matrix* K_red = new My_Matrix();
  K_red->ReDimension(n,c); 
  for(i=0; i<n; i++){
    ii = 0;
    for(j=0; j<kc; j++){
      if(index[j] == 0)
	continue;
      
      K_red->el(i,ii) = A->el(i,j);
      ii +=1;
    }
  }

  double* res_c = new double[c];
  kk =0;
  for(j=0; j<kc; j++){
      if(index[j] == 0)
	  continue;
      res_c[kk] = beta[j];
      kk++;
  }


  ii = 0;
  double* theta_vec = new double[c];
  for(i=0;i<kc;i++){
    if(index[i] == 0)
      continue;
    theta_vec[ii] = Theta[i];
    ii +=1;
  }
 
  double epsi_tilde =0.;
  for(i=0;i<kc;i++){
    if(index[i] == 0)
      continue;
    epsi_tilde -=  Theta[i]*beta[i];
  }

  epsi_tilde += mu;

  
  getCoeffnew(n, K_red,   YTd, res_c, epsi_tilde , c, theta_vec);
   
  double* optVal; 
  
  optVal = optimize_qp();
 
  ii = 0;
  for(i=0;i<kc;i++){
      if(index[i] == 0){
	  h[i] = 0.0;
	  continue;
      }
      h[i] = optVal[ii];
      ii +=1;
  }

 

  for(i=0;i<kc;i++)
      old_index[i] = index[i];
 
  delete K_red;
  
  delete [] theta_vec;
  delete [] res_c;
}


 
void LogregLASSO::regress()
{
  int i;
  double* Evec = new double[kc];

  SUBlasso(Evec, Y);  
 
  double* vec3 = new double[n]; 

  multA_Beta_index(Evec,vec3);
  
  for(i=0; i<n;i++)
      Projected[i] = vec3[i];
  
  delete [] vec3;
   
  for(i=0; i<kc; i++)
      alpha[i] = Evec[i];
  delete [] Evec;
  
}