www.pudn.com > SVM.rar > pr_loqo.c

/*
 * File:        pr_loqo.c
 * Purpose:     solves quadratic programming problem for pattern recognition
 *              for support vectors
 *
 * Author:      Alex J. Smola
 * Created:     10/14/97
 * Updated:     11/08/97
 * Updated:     13/08/98 (removed exit(1) as it crashes svm lite when the margin
 *                        in a not sufficiently conservative manner)
 *
 * 
 * Copyright (c) 1997  GMD Berlin - All rights reserved
 * THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE of GMD Berlin
 * The copyright notice above does not evidence any
 * actual or intended publication of this work.
 *
 * Unauthorized commercial use of this software is not allowed
 */

#include 
#include 
#include 
#include 
#include "pr_loqo.h"

#define	max(A, B)	((A) > (B) ? (A) : (B))
#define	min(A, B)	((A) < (B) ? (A) : (B))
#define sqr(A)          ((A) * (A))
#define	ABS(A)  	((A) > 0 ? (A) : (-(A)))

#define PREDICTOR 1
#define CORRECTOR 2

/*****************************************************************
  replace this by any other function that will exit gracefully
  in a larger system
  ***************************************************************/

void nrerror(char error_text[])
{
  printf("ERROR: terminating optimizer - %s\n", error_text);
  /* exit(1); */
}

/*****************************************************************
   taken from numerical recipes and modified to accept pointers
   moreover numerical recipes code seems to be buggy (at least the
   ones on the web)

   cholesky solver and backsubstitution
   leaves upper right triangle intact (rows first order)
   ***************************************************************/

void choldc(double a[], int n, double p[])
{
  void nrerror(char error_text[]);
  int i, j, k;
  double sum;

  for (i = 0; i < n; i++){
    for (j = i; j < n; j++) {
      sum=a[n*i + j];
      for (k=i-1; k>=0; k--) sum -= a[n*i + k]*a[n*j + k];
      if (i == j) {
	if (sum <= 0.0) {
	  nrerror("choldc failed, matrix not positive definite");
	  sum = 0.0;
	}
	p[i]=sqrt(sum);
      } else a[n*j + i] = sum/p[i];
    }
  }
}

void cholsb(double a[], int n, double p[], double b[], double x[])
{
  int i, k;
  double sum;

  for (i=0; i=0; k--) sum -= a[n*i + k]*x[k];
    x[i]=sum/p[i];
  }

  for (i=n-1; i>=0; i--) {
    sum=x[i];
    for (k=i+1; k=0; k--) sum -= a[n*i + k]*x[k];
    x[i]=sum/p[i];
  }
}

void chol_backward(double a[], int n, double p[], double b[], double x[])
{
  int i, k;
  double sum;

  for (i=n-1; i>=0; i--) {
    sum=b[i];
    for (k=i+1; k 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= 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 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= 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;
}