loqo_solver.cpp

支持向量分类算法在linux操作系统下的是实现

      ierr = PLA_Obj_local_buffer(x_loc_view, (void **) &addr_buf);
      CheckError(ierr);
      alpha_new = (LOQOfloat *) addr_buf;
      ierr = PLA_Obj_local_buffer(yy_loc, (void **) &addr_buf);
      CheckError(ierr);
      dual = (LOQOfloat *) addr_buf;
      ierr = PLA_Obj_local_buffer(dist_loc_view, (void **) &addr_buf);
      dist_ = (LOQOfloat *) addr_buf;

      ++iter;
      // Setup problem for parallel LOQO solver
      time = MPI_Wtime();
      info("setup_problem...");
      setup_problem(work_set, not_work_set);
      info("done.\n"); info_flush();
      info("setup_up...");
      setup_up(work_set);
      info("done.\n"); info_flush();
      info("setup_low...");
      setup_low();
      info("done.\n"); info_flush();
      //print_problem();
      time = MPI_Wtime() - time;
      problem_setup_time += time;

      // Compute obj before optimization
      info("Computing obj before..."); info_flush();
      obj_before = 0.0;
      int start_i = 0;
      for(int i=0; i<n; ++i, ++start_i)
	{
	  if(work_set[i] >= l_low_loc)
	    break;
	}
      //      for(int i=n_low_loc; i<n_up_loc; ++i)
      for(int i=start_i; work_set[i] < l_up_loc && i<n; ++i)
	{
	  //obj_before += alpha[work_set[i]]*c[i-n_low_loc];
	  obj_before += alpha[work_set[i]]*c[i-start_i];
 	  for(int j=0; j<n; ++j)
//  	    obj_before += 0.5*alpha[work_set[i]]*Q_bb[local_n*j+(i-n_low_loc)]*
// 	      alpha[work_set[j]];
 	    obj_before += 0.5*alpha[work_set[i]]*Q_bb[num_rows*j+(i-start_i)]*
	      alpha[work_set[j]];
	}
      // Get local changes for obj_before from other processors
      double obj_before_loc=0.0;
      double obj_before_other=0.0;
      for(int i=0; i<size; ++i)
	{
	  if(i == rank)
	    ierr = MPI_Bcast(&obj_before, 1, MPI_DOUBLE, i, comm);
	  else
	    {
	      ierr = MPI_Bcast(&obj_before_loc, 1, MPI_DOUBLE, i, comm);
	      obj_before_other += obj_before_loc;
	    }
	}
      obj_before += obj_before_other;
      printf("obj before = %g\n", obj_before); fflush(stdout);
      info("done.\n"); info_flush();

      // Run solver
      time = MPI_Wtime();
      status = solve_inner(work_set);
      time = MPI_Wtime() - time;
      inner_solver_time += time;

      // Update gradient.
      time = MPI_Wtime();
      int *nz = new int[n];
      for(int i=0; i<n; ++i)
	{
	  delta_alpha[i] = alpha_new[i] - alpha[work_set[i]];
	  if(fabs(delta_alpha[i]) > TOL_ZERO)
	    nz[i] = 1; 
	  else
	    nz[i] = 0;
	}
      info("Updating G_b..."); info_flush();
      // Compute G_b
      for(int i=start_i; work_set[i] < l_up_loc && i<n; ++i)
      // for(int i=n_low_loc; i<n_up_loc; ++i)
	{
 	  for(int j=0; j<n; ++j)
	    {
 	      if(nz[j])
// 		G[work_set[i]] += Q_bb[j*local_n+(i-n_low_loc)] * 
// 		  delta_alpha[j];
		G[work_set[i]] += Q_bb[j*num_rows+(i-start_i)] * 
		  delta_alpha[j];
	    }
	}
      info("done.\n"); info_flush();

      // Update local part of the parallel cache status
      for(int i=0; i<l; ++i)
	{
	  if(Q.is_cached(i))
	    p_cache_status[rank*size + i] = CACHED;
	  else
	    p_cache_status[rank*size + i] = NOT_CACHED;
	}
      // Synchronize parallel cache status
      for(int k=0; k<size; ++k)
	{
	  ierr = MPI_Bcast(&p_cache_status[k*size], l, MPI_CHAR, k, comm); 
	  CheckError(ierr);	  
	}
//       printf("p_cache_status %d =",rank);
//       for(int k=0; k<size; ++k)
// 	{
// 	  for(int i=0; i<l; ++i)
// 	    {
// 	      printf(" %d", p_cache_status[k*size+i]);
// 	    }
// 	  printf("\n");
// 	}
//       printf("\n");

      // Smart parallel cache handling
      int *idx_cached = new int[n];
      int *idx_not_cached = new int[n];
      int count_cached = 0; int count_not_cached = 0;
      int next_k = 0; bool found = false;
      int next_not_cached = 0;
      for(int i=0; i<n; ++i)
	{
	  if(nz[i])
	    {
	      for(int k=next_k; !found && k<size; ++k)
		{
		  if(p_cache_status[k*size + work_set[i]] == CACHED)
		    {
		      if(k == rank) // Do we have it?
			{
			  idx_cached[count_cached] = i;
			  ++count_cached;
			}
		      found = true;
		      next_k = k == size-1 ? 0 : k+1;
		    }
		}
	      for(int k=0; !found && k<next_k; ++k)
		{
		  if(p_cache_status[k*size + work_set[i]] == CACHED)
		    {
		      if(k == rank) // Do we have it?
			{
			  idx_cached[count_cached] = i;
			  ++count_cached;
			}
		      found = true;
		      next_k = k == size-1 ? 0 : k+1;
		    } 
		}
	      if(!found) // not in any cache
		{
		  if(rank == next_not_cached) // Do we have to compute it?
		    {
		      idx_not_cached[count_not_cached] = i;
		      ++count_not_cached;
		    }
		  next_not_cached = next_not_cached == size-1 ? 0 : 
		    next_not_cached + 1;
		}
	      found = false;
	    } // if(nz[i])
	}

      info("Updating G_n..."); info_flush();
      // Compute G_n
      for(int j=0; j<lmn; ++j)
	G_n[j] = 0;
//       printf("idx_cached %d =", rank);
//       for(int i=0; i<count_cached; ++i)
// 	printf(" %d", idx_cached[i]);
//       printf("\n");
//       printf("idx_not_cached %d =", rank);
//       for(int i=0; i<count_not_cached; ++i)
// 	printf(" %d", idx_not_cached[i]);
//       printf("\n");
      // First update the cached part...
      for(int i=0; i<count_cached; ++i)
	{
	  const Qfloat *Q_i = Q.get_Q_subset(work_set[idx_cached[i]],
					      not_work_set,lmn);
	  for(int j=0; j<lmn; ++j)
	    G_n[j] += Q_i[not_work_set[j]] * delta_alpha[idx_cached[i]];
	}
      // now update the non-cached part
      for(int i=0; i<count_not_cached; ++i)
	{
	  const Qfloat *Q_i = Q.get_Q_subset(work_set[idx_not_cached[i]],
					      not_work_set,lmn);
	  for(int j=0; j<lmn; ++j)
	    G_n[j] += Q_i[not_work_set[j]] * delta_alpha[idx_not_cached[i]];
	}
//       for(int i=start_i; work_set[i] < l_up_loc && i<n; ++i)
// 	 //      for(int i=n_low_loc; i<n_up_loc; ++i)
// 	{
// 	  if(nz[i]){
// 	    const Qfloat *Q_i = Q.get_Q_subset(work_set[i],not_work_set,lmn);
// 	    for(int j=0; j<lmn; ++j)
// 	      G_n[j] += Q_i[not_work_set[j]] * delta_alpha[i];

// 	  }
// 	}
      
      info("done.\n"); info_flush();
      delete[] idx_cached;
      delete[] idx_not_cached;
      delete[] nz;
      time = MPI_Wtime() - time;
      gradient_updating_time += time;

      // Synchronize gradient with other processors
      info("Synchronizing gradient..."); info_flush();
      sync_gradient(work_set, not_work_set);
      info("done.\n"); info_flush();

//       info("synced G proc[%d]=",rank);
//       for(int i=0; i<l; ++i)
// 	info(" %g", G[i]);
//       info("\n");
//       info_flush();
      // Update alpha
      for(int i=0; i<n; ++i)
	{
	  alpha[work_set[i]] = alpha_new[i];
	  update_alpha_status(work_set[i]);
	}
    } // while(1)
  // Calculate rho
  si->rho = calculate_rho();

  // Calculate objective value
  {
    double v = 0;
    int i;
    for(i=0;i<l;i++)
      v += alpha[i] * (G[i] + b[i]);

    si->obj = v/2;
  }

  // Put back the solution
  {
//     printf("alpha = ");
    for(int i=0;i<l;i++)
      {
	alpha_[i] = alpha[i];
// 	printf(" %g\n",alpha[i]);
      }
//     printf("\n");
  }

  si->upper_bound_p = Cp;
  si->upper_bound_n = Cn;

  info("\noptimization finished, #iter = %d\n",iter);

  total_time = MPI_Wtime() - total_time;

  // print timing statistics
  if(rank == 0)
    {
      info("Total opt. time = %lf\n", total_time); info_flush();
      info("Problem setup time = %lf (%lf%%)\n", problem_setup_time,
	   problem_setup_time/total_time*100); info_flush();
      info("Inner solver time = %lf (%lf%%)\n", inner_solver_time,
	   inner_solver_time/total_time*100); info_flush();
      info("Gradient updating time = %lf (%lf%%)\n", gradient_updating_time,
	   gradient_updating_time/total_time*100); info_flush();
    }

  // Clean up
  delete[] b;
  delete[] y;
  delete[] alpha;
  delete[] alpha_status;
  delete[] delta_alpha;
  delete[] work_status;
  delete[] work_count;
  delete[] work_set;
  delete[] not_work_set;
  delete[] G;
  delete[] Q_bb;
  delete[] c;
  delete[] up;
  delete[] low;
  delete[] a;
  delete[] d;
  delete[] l_low;
  delete[] l_up;
  delete[] n_low;
  delete[] n_up;
  delete[] lmn_low;
  delete[] lmn_up;
  delete[] G_n;
  delete[] p_cache_status;
  ierr = PLA_Obj_free(&Q_bb_global); CheckError(ierr);
  ierr = PLA_Obj_free(&c_global); CheckError(ierr);
  ierr = PLA_Obj_free(&up_global); CheckError(ierr);
  ierr = PLA_Obj_free(&low_global); CheckError(ierr);
  ierr = PLA_Obj_free(&a_global); CheckError(ierr);
  ierr = PLA_Obj_free(&d_global); CheckError(ierr);
  ierr = PLA_Obj_free(&x_loc); CheckError(ierr);
  ierr = PLA_Obj_free(&x); CheckError(ierr);
  ierr = PLA_Obj_free(&g); CheckError(ierr);
  ierr = PLA_Obj_free(&t); CheckError(ierr);
  ierr = PLA_Obj_free(&yy); CheckError(ierr)
  ierr = PLA_Obj_free(&yy_loc); CheckError(ierr);
  ierr = PLA_Obj_free(&z); CheckError(ierr);
  ierr = PLA_Obj_free(&s); CheckError(ierr);
}

void Solver_Parallel_LOQO::sync_gradient(int *work_set, int *not_work_set)
{
  int start_i = 0;
  for(int i=0; i<n; ++i, ++start_i)
    {
      if(work_set[i] >= l_low_loc)
	break;
    }
  // Synchronize G_b
  //  double *G_send = new double[local_n];
  double *G_send = new double[num_rows];
  double *G_recv = new double[n];
  int count=0;
  for(int i=start_i; work_set[i] < l_up_loc && i<n; ++i)
    //for(int i=n_low_loc; i<n_up_loc; ++i)
    {
      G_send[count] = G[work_set[i]];
      ++count;
    }
//   printf("G_send = ");
//   for(int i=start_i; work_set[i] < l_up_loc && i<n; ++i)
//     printf(" %g", G_send[i]);
//   printf("\n");
//   printf("count = %d\n", count);
//   printf("num_rows = %d\n", num_rows);
  int *sendcounts = new int[size];
  int *sdispls = new int[size];
  for(int i=0; i<size; ++i)
    {
      //      sendcounts[i] = local_n;
      sendcounts[i] = num_rows;
      sdispls[i] = 0;
    }
  int *recvcounts = new int[size];
  int *rdispls = new int[size];

  for(int k=0; k<size; ++k)
    {
      int start_i_other = 0;
      for(int i=0; i<n; ++i, ++start_i_other)
	if(work_set[i] >= l_low[k])
	  break;
      int num_rows_other = 0;
      for(int i=start_i_other; work_set[i] < l_up[k] && i<n; ++i)
	++num_rows_other;
      recvcounts[k] = num_rows_other;
      if(k == 0)
	rdispls[k] = 0;
      else
	rdispls[k] = rdispls[k-1] + recvcounts[k-1];
    }
//   recvcounts[0] = n_up[0] - n_low[0];
//   rdispls[0] = 0;
//   for(int i=1; i<size; ++i)
//     {
//       recvcounts[i] = n_up[i] - n_low[i];
//       rdispls[i] = rdispls[i-1] + recvcounts[i-1];
//     }
  ierr = MPI_Alltoallv(G_send, sendcounts, sdispls, MPI_DOUBLE,
		       G_recv, recvcounts, rdispls, MPI_DOUBLE, comm);
  CheckError(ierr);
  // Update local G_b
  for(int k=0; k<size; ++k) 
    {
      if(k != rank)
	{
	  int start_i_other = 0;
	  for(int i=0; i<n; ++i, ++start_i_other)
	    if(work_set[i] >= l_low[k])
	      break;
	  // G_b
	  for(int i=start_i_other; work_set[i] < l_up[k] && i<n; ++i)
	    G[work_set[i]] = G_recv[rdispls[k]+(i-start_i_other)];
// 	  for(int j=n_low[i]; j<n_up[i]; ++j)
// 	    G[work_set[j]] = G_recv[rdispls[i]+(j-n_low[i])];
	}
    }
  delete[] G_send;
  delete[] G_recv;
  delete[] sendcounts;
  delete[] sdispls;
  delete[] recvcounts;
  delete[] rdispls;

  double *G_buf = new double[lmn];
  // Synchronize G_n
  for(int i=0; i<size; ++i)
    {
      if(rank == i)
	{
	  for(int j=0; j<lmn; ++j)
	    G_buf[j] = G_n[j];
	}
      ierr = MPI_Bcast(G_buf, lmn, MPI_DOUBLE, i, comm); CheckError(ierr);
      // Accumulate contributions
      for(int j=0; j<lmn; ++j)
	G[not_work_set[j]] += G_buf[j];
    }
  delete[] G_buf;
}

int Solver_LOQO::select_working_set(int *work_set, int *not_work_set)
{
 printf("selecting working set...");
  // reset work status
  n = n_old;
  for(int i=0; i<l; ++i)
      work_status[i] = WORK_N;

  double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }
  double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }
  for(int i=0; i<l; ++i)
    {
      if(!is_upper_bound(i))
	{
	  if(y[i] == +1)
	    {
	      if(-G[i] > Gmax1)
		Gmax1 = -G[i];
	    }
	  else
	    {
	      if(-G[i] > Gmax2)
		Gmax2 = -G[i];
	    }
	}
      if(!is_lower_bound(i))
	{
	  if(y[i] == +1)
	    {
	      if(G[i] > Gmax2)
		Gmax2 = G[i];
	    }
	  else
	    {
	      if(G[i] > Gmax1)
		Gmax1 = G[i];
	    }
	}
    }
  // check for optimality
  printf("Gmax1 + Gmax2 = %g < %g\n", Gmax1+Gmax2,eps);
  if(Gmax1 + Gmax2 < eps)
      return 1;

  // Compute yG
  double *yG = new double[l];
  int *pidx = new int[l];
  for(int i=0; i<l; ++i)
    {
      if(y[i] == +1)
	yG[i] = G[i];
      else
	yG[i] = -G[i];
      pidx[i] = i;
    }
  quick_sort(yG, pidx, 0, l-1);
//   printf("yG = ");
//   for(int i=0; i<l; ++i)
//     printf(" %g",yG[i]);
//   printf("\n");
//   printf("pidx = ");
//   for(int i=0; i<l; ++i)
//     printf(" %d",pidx[i]);
//   printf("\n");
  int top=l-1; int bot=0; int count=0;
  // Select a full set initially
  int nselect = iter == 0 ? n : q;
  while(top > bot && count < nselect)
    {
      while(!(is_free(pidx[top]) 
	      || (is_upper_bound(pidx[top]) && y[pidx[top]] == +1)
	      || (is_lower_bound(pidx[top]) && y[pidx[top]] == -1)
	      ))
	--top;
      while(!(is_free(pidx[bot])
	      || (is_upper_bound(pidx[bot]) && y[pidx[bot]] == -1)
	      || (is_lower_bound(pidx[bot]) && y[pidx[bot]] == +1)
	      ))
	++bot;
      if(top > bot)