pr_loqo.c

MATLAB的SVM算法实现

  tau       = malloc(n*sizeof(double));
  sigma     = malloc(n*sizeof(double));

  gamma_z   = malloc(n*sizeof(double));
  gamma_s   = malloc(n*sizeof(double));

  hat_nu    = malloc(n*sizeof(double));
  hat_tau   = malloc(n*sizeof(double));

  delta_x   = malloc(n*sizeof(double));
  delta_y   = malloc(m*sizeof(double));
  delta_s   = malloc(n*sizeof(double));
  delta_z   = malloc(n*sizeof(double));
  delta_g   = malloc(n*sizeof(double));
  delta_t   = malloc(n*sizeof(double));

  d         = malloc(n*sizeof(double));

  /* pointers into the external variables */
  x = primal;			/* n */
  g = x + n;			/* n */
  t = g + n;			/* n */

  y = dual;			/* m */
  z = y + m;			/* n */
  s = z + n;			/* n */

  /* initial settings */
  b_plus_1 = 1;
  c_plus_1 = 0;
  for (i=0; i<n; i++) c_plus_1 += c[i];

  /* get diagonal terms */
  for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i]; 
  
  /* starting point */
  if (restart == 1) {
				/* x, y already preset */
    for (i=0; i<n; i++) {	/* compute g, t for primal feasibility */
      g[i] = max(ABS(x[i] - l[i]), bound);
      t[i] = max(ABS(u[i] - x[i]), bound); 
    }

    matrix_vector(n, h_x, x, h_dot_x); /* h_dot_x = h_x * x */

    for (i=0; i<n; i++) {	/* sigma is a dummy variable to calculate z, s */
      sigma[i] = c[i] + h_dot_x[i];
      for (j=0; j<m; j++)
	sigma[i] -= a[n*j + i] * y[j];

      if (sigma[i] > 0) {
	s[i] = bound;
	z[i] = sigma[i] + bound;
      }
      else {
	s[i] = bound - sigma[i];
	z[i] = bound;
      }
    }
  }
  else {			/* use default start settings */
    for (i=0; i<m; i++)
      for (j=i; j<m; j++)
	h_y[i*m + j] = (i==j) ? 1 : 0;
    
    for (i=0; i<n; i++) {
      c_x[i] = c[i];
      h_x[(n+1)*i] += 1;
    }
    
    for (i=0; i<m; i++)
      c_y[i] = b[i];
    
    /* and solve the system [-H_x A'; A H_y] [x, y] = [c_x; c_y] */
    solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,
		  PREDICTOR);
    
    /* initialize the other variables */
    for (i=0; i<n; i++) {
      g[i] = max(ABS(x[i] - l[i]), bound);
      z[i] = max(ABS(x[i]), bound);
      t[i] = max(ABS(u[i] - x[i]), bound); 
      s[i] = max(ABS(x[i]), bound); 
    }
  }

  for (i=0, mu=0; i<n; i++)
    mu += z[i] * g[i] + s[i] * t[i];
  mu = mu / (2*n);

  /* the main loop */
  if (verb >= STATUS) {
    printf("counter | pri_inf  | dual_inf  | pri_obj   | dual_obj  | ");
    printf("sigfig | alpha  | nu \n");
    printf("-------------------------------------------------------");
    printf("---------------------------\n");
  }
  
  while (status == STILL_RUNNING) {
    /* predictor */
    
    /* put back original diagonal values */
    for (i=0; i<n; i++) 
      h_x[(n+1) * i] = diag_h_x[i];

    matrix_vector(n, h_x, x, h_dot_x); /* compute h_dot_x = h_x * x */

    for (i=0; i<m; i++) {
      rho[i] = b[i];
      for (j=0; j<n; j++)
	rho[i] -= a[n*i + j] * x[j];
    }
    
    for (i=0; i<n; i++) {
      nu[i] = l[i] - x[i] + g[i];
      tau[i] = u[i] - x[i] - t[i];

      sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];
      for (j=0; j<m; j++)
	sigma[i] -= a[n*j + i] * y[j];

      gamma_z[i] = - z[i];
      gamma_s[i] = - s[i];
    }

    /* instrumentation */
    x_h_x = 0;
    primal_inf = 0;
    dual_inf = 0;
    
    for (i=0; i<n; i++) {
      x_h_x += h_dot_x[i] * x[i];
      primal_inf += sqr(tau[i]);
      primal_inf += sqr(nu[i]);
      dual_inf += sqr(sigma[i]);
    }
    for (i=0; i<m; i++) 
      primal_inf += sqr(rho[i]);
    primal_inf = sqrt(primal_inf)/b_plus_1;
    dual_inf = sqrt(dual_inf)/c_plus_1;

    primal_obj = 0.5 * x_h_x;
    dual_obj = -0.5 * x_h_x;
    for (i=0; i<n; i++) {
      primal_obj += c[i] * x[i];
      dual_obj += l[i] * z[i] - u[i] * s[i];
    }
    for (i=0; i<m; i++)
      dual_obj += b[i] * y[i];

    sigfig = log10(ABS(primal_obj) + 1) -
             log10(ABS(primal_obj - dual_obj));
    sigfig = max(sigfig, 0);
		   
    /* the diagnostics - after we computed our results we will
       analyze them */

    if (counter > counter_max) status = ITERATION_LIMIT;
    if (sigfig  > sigfig_max)  status = OPTIMAL_SOLUTION;
    if (primal_inf > 10e100)   status = PRIMAL_INFEASIBLE;
    if (dual_inf > 10e100)     status = DUAL_INFEASIBLE;
    if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE;
    if (ABS(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;
    if (ABS(dual_obj) > 10e100) status = DUAL_UNBOUNDED;

    /* write some nice routine to enforce the time limit if you
       _really_ want, however it's quite useless as you can compute
       the time from the maximum number of iterations as every
       iteration costs one cholesky decomposition plus a couple of 
       backsubstitutions */

    /* generate report */
    if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)))
      printf("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",
	     counter, primal_inf, dual_inf, primal_obj, dual_obj,
	     sigfig, alfa, mu);

    counter++;

    if (status == 0) {		/* we may keep on going, otherwise
				   it'll cost one loop extra plus a
				   messed up main diagonal of h_x */
      /* intermediate variables (the ones with hat) */
      for (i=0; i<n; i++) {
	hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
	hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
	/* diagonal terms */
	d[i] = z[i] / g[i] + s[i] / t[i];
      }
      
      /* initialization before the cholesky solver */
      for (i=0; i<n; i++) {
	h_x[(n+1)*i] = diag_h_x[i] + d[i];
	c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - 
	  s[i] * hat_tau[i] / t[i];
      }
      for (i=0; i<m; i++) {
	c_y[i] = rho[i];
	for (j=i; j<m; j++) 
	  h_y[m*i + j] = 0;
      }
      
      /* and do it */
      solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
		    PREDICTOR);
      
      for (i=0; i<n; i++) {
	/* backsubstitution */
	delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
	delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
	
	delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
	delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
	
	/* central path (corrector) */
	gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];
	gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];
	
	/* (some more intermediate variables) the hat variables */
	hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
	hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
	
	/* initialization before the cholesky */
	c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];
      }
      
      for (i=0; i<m; i++) {	/* comput c_y and rho */
	c_y[i] = rho[i];
	for (j=i; j<m; j++)
	  h_y[m*i + j] = 0;
      }
      
      /* and do it */
      solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
		    CORRECTOR);
      
      for (i=0; i<n; i++) {
	/* backsubstitution */
	delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
	delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
	
	delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
	delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
      }
      
      alfa = -1;
      for (i=0; i<n; i++) {
	alfa = min(alfa, delta_g[i]/g[i]);
	alfa = min(alfa, delta_t[i]/t[i]);
	alfa = min(alfa, delta_s[i]/s[i]);
	alfa = min(alfa, delta_z[i]/z[i]);
      }
      alfa = (margin - 1) / alfa;
      
      /* compute mu */
      for (i=0, mu=0; i<n; i++)
	mu += z[i] * g[i] + s[i] * t[i];
      mu = mu / (2*n);
      mu = mu * sqr((alfa - 1) / (alfa + 10));
      
      for (i=0; i<n; i++) {
	x[i] += alfa * delta_x[i];
	g[i] += alfa * delta_g[i];
	t[i] += alfa * delta_t[i];
	z[i] += alfa * delta_z[i];
	s[i] += alfa * delta_s[i];
      }
      
      for (i=0; i<m; i++) 
	y[i] += alfa * delta_y[i];
    }
  }
  if ((status == 1) && (verb >= STATUS)) {
    printf("----------------------------------------------------------------------------------\n");
    printf("optimization converged\n");
  }
  
  /* free memory */
  free(workspace);
  free(diag_h_x);
  free(h_y);
  free(c_x);
  free(c_y);
  free(h_dot_x);
  
  free(rho);
  free(nu);
  free(tau);
  free(sigma);
  free(gamma_z);
  free(gamma_s);
  
  free(hat_nu);
  free(hat_tau);
    
  free(delta_x);
  free(delta_y);
  free(delta_s);
  free(delta_z);
  free(delta_g);
  free(delta_t);
    
  free(d);

  /* and return to sender */
  return status;
}