hs046.c

优化算法loqo的算法源代码。Purpose: solves quadratic programming problem for pattern recognition for support vec

/*********************************************************************/
/***    Copyright (c) Robert J. Vanderbei, 1994                    ***/
/***    All Rights Reserved                                        ***/
/*********************************************************************/

/* This C-source file illustrates how to solve the Hock and Schittkowski 
   problem number 46 using the LOQO subroutine library.  
   The AMPL model for this problem is:

	minimize obj: 
	  (x[1]-x[2])^2 + (x[3]-1)^2 + (x[4]-1)^4 + (x[5]-1)^6
	  ;

	subject to constr1: x[1]^2*x[4] + sin(x[4] - x[5]) = 1;
	subject to constr2: x[2] + x[3]^4*x[4]^2 = 2;

   In LOQO, all numbering starts from 0.  Therefore, it is convenient
   to record here the objective and constraints using 0-based indexing:

	minimize obj: 
	  (x[0]-x[1])^2 + (x[2]-1)^2 + (x[3]-1)^4 + (x[4]-1)^6
	  ;

	subject to constr0: x[0]^2*x[3] + sin(x[3] - x[4]) = 1;
	subject to constr1: x[1] + x[2]^4*x[3]^2 = 2;
   */

#include <string.h>
#include <math.h>
#include <stdlib.h>

#ifdef __cplusplus
extern "C" {
#endif

#include "loqo.h"
#include "myalloc.h"

void nlinit( void *vlp );
int nlupdate( void *vlp );

double objval( double *x );
void objgrad( double *c, double *x );
void hessian( double *Q, double *x, double *y );
void conval( double *h, double *x );
void congrad( double *A, double *At, double *x );

main()
{
	int 	status;
	char	fname[80];	/* solution file name */
	LOQO	*lp;

	int i,j,k;

	int bad_opns;
	char *opnstr;

	lp = openlp();
	if (lp == NULL) { my_exit(11,"cannot open another LP problem"); }

	/* The LOQO data structure pointed to by lp contains all the
	   information about a problem.  It consists of many components.
	   They are listed in loqo.h.  Every component has some sort
	   of reasonable default and therefore only those for which the
	   default doesn't match the problem at hand must be set explicitly.
	   The most obvious examples of components that don't match the defaults
	   are the number of variables and the number of constraints.
	   So, we start be setting these values and some others.
	   */

	lp->n = 5;	/* number of variables   (indexed 0,1,...,n-1) */
	lp->m = 2;      /* number of constraints (indexed 0,1,...,m-1) */
	lp->nz = 6; 	/* number of nonzeros in linearized constraint matrix */
	lp->qnz = 13;   /* number of nonzeros in Hessian (details below) */

	/* Most of the data defining the quadratic approximation to a 
	   given nonlinear problem will be changed at each iteration.
	   Later we will define subroutines to do that.  Here, however,
	   we must allocate storage for the arrays and we must set up the
	   sparse representation of the matrices.  We begin by allocating
	   storage.
	   */

	/* The matrix containing the linearized constraints is called A.
	   It is stored sparsely in 3 arrays:  A, iA, kA.
	   Array A contains the values of the nonzeros list one column after
	   another.
	   Array iA contains the row index of each corresponding value in A.
	   Array kA contains a list of indices in the A/iA array indicating
	   where each new column starts.  Hence, kA[0]=0, kA[n]=nz, and
	   kA[j+1]-kA[j] is the number of nonzero elements in column j.
	   */
	MALLOC(  lp->A,  lp->nz, double );
	MALLOC(  lp->iA, lp->nz,  int );
	MALLOC(  lp->kA, lp->n+1, int );

	/* The Hessian is called Q.  
	   It is stored sparsely in 3 arrays: Q, iQ, kQ.
	   */
	MALLOC(  lp->Q,  lp->qnz, double );
	MALLOC(  lp->iQ, lp->qnz,  int );
	MALLOC(  lp->kQ, lp->n+1, int );

	/* The vector of constant terms in the linearized constraints is
	   moved to the right-hand side of the equation and is called b.
	   It is stored as a dense (i.e., ordinary) array.
	   Since the rhs is assumed to be constant, we can initialize its
	   values here.
	   */
	MALLOC( lp->b, lp->m, double );
	lp->b[0] = 1;
	lp->b[1] = 2;

	/* The vector of coefficients on the linear terms of the quadratic
	   approximation to the objective are stored in an array c.
	   It is stored as a dense (i.e., ordinary) array.
	   */
	MALLOC( lp->c, lp->n, double );

	/* The vector of lower bounds on the variables is called l.  If
	   l==NULL, then solvelp() will set it to a zero vector.  Since the
	   variables in this problem are free, we reset this vector as follows.
	   */
	MALLOC( lp->l, lp->n, double );
	for (j=0; j<lp->n; j++) { lp->l[j] = -HUGE_VAL; }

	/* The vector of upper bounds on the variables is called u.  If
	   u==NULL, then solvelp() will set it to a vector of HUGE_VAL's.
	   This default behavior is correct for the problem at hand and so
	   nothing needs to be done.  Nonetheless, we set it here so one
	   can see how it is done.
	   */
	MALLOC( lp->u, lp->n, double );
	for (j=0; j<lp->n; j++) { lp->u[j] =  HUGE_VAL; }

	/* Each constraint is assumed to be an inequality constraint with
	   a range, r, on the inequality:
		b <= a_i(x) <= b+r .
	   An equality constraint is specified by setting r[i] to 0 (this
	   is the default).  
	   An inequality constraint is specified by setting it to HUGE_VAL.
	   Since the default behavior is easy to forget, it is a good idea
	   to set r explicitly as we do here.
	   */
	MALLOC( lp->r, lp->m, double );
	for (i=0; i<lp->m; i++) { lp->r[i] = 0; }

	/* Now, we define the sparse structure of A. 
	   The sparsity pattern is as follows:
		  (var  1  2  3  4  5)
		   col  0  1  2  3  4
		row
		 0      *        *  *
		 1         *  *  *
	   */
	lp->kA[0] = 0;
			lp->iA[0] = 0;
	lp->kA[1] = 1;
			lp->iA[1] = 1;
	lp->kA[2] = 2;
			lp->iA[2] = 1;
	lp->kA[3] = 3;
			lp->iA[3] = 0;
			lp->iA[4] = 1;
	lp->kA[4] = 5;
			lp->iA[5] = 0;
	lp->kA[5] = 6;

	/* Next, we define the sparse structure of Q.  Since Q involves
	   quadratic information about the objective function and each
	   nonlinear constraint, it must be the union of these sparsity 
	   patterns.  The sparsity pattern for the Hessian of the objective
	   function is:
		       (var  1  2  3  4  5)
		        col  0  1  2  3  4
		var  row
		 1    0      *  *         
		 2    1      *  *           
		 3    2            *      
		 4    3               *   
		 5    4                  *

	   For the first constraint, the pattern is as follows:
		       (var  1  2  3  4  5)
		        col  0  1  2  3  4
		var  row
		 1    0      *        *   
		 2    1                     
		 3    2                   
		 4    3      *        *  *
		 5    4               *  *

	   For the second constraint, the pattern is as follows:
		       (var  1  2  3  4  5)
		        col  0  1  2  3  4
		var  row
		 1    0                   
		 2    1                     
		 3    2            *  *   
		 4    3            *  *   
		 5    4                   

	   The union of these three patterns is this:
		       (var  1  2  3  4  5)
		        col  0  1  2  3  4
		var  row
		 1    0      *  *     *   
		 2    1      *  *           
		 3    2            *  *   
		 4    3      *     *  *  *
		 5    4               *  *

	   We store this pattern in iQ,kQ as follows.
	   */
	lp->kQ[0] = 0;
			lp->iQ[0]  = 0;
			lp->iQ[1]  = 1;
			lp->iQ[2]  = 3;
	lp->kQ[1] = 3;
			lp->iQ[3]  = 0;
			lp->iQ[4]  = 1;
	lp->kQ[2] = 5;
			lp->iQ[5]  = 2;
			lp->iQ[6]  = 3;
	lp->kQ[3] = 7;
			lp->iQ[7]  = 0;
			lp->iQ[8]  = 2;
			lp->iQ[9]  = 3;
			lp->iQ[10] = 4;
	lp->kQ[4] = 11;
			lp->iQ[11] = 3;
			lp->iQ[12] = 4;
	lp->kQ[5] = 13;

	/* Specific data for A, c, and Q will be computed in nlupdate(),
	   defined later.  For now, we fill these arrays with zeros.
	   */
	for (k=0; k<lp->nz;  k++) { lp->A[k] = 0; }
	for (j=0; j<lp->n;   j++) { lp->c[j] = 0; }
	for (k=0; k<lp->qnz; k++) { lp->Q[k] = 0; }

	/* The function that updates the quadratic approximation
	   to the nonlinear problem is defined later in this file.  Here, we
	   simply store a pointer to it in the LOQO data structure.
	   We also need to have a once called routine nlinit() saved for
	   future use.
	   */
	nlsetup(lp);


	/* solvelp() is happy to compute its own default starting point but 
	   in NLP it is customary for the user to supply one.  Here is 
	   a starting point.
	   */
	MALLOC( lp->x, lp->n, double );
	lp->x[0] = sqrt(2.)/2;
	lp->x[1] = 1.75;
	lp->x[2] = 0.5;
	lp->x[3] = 2;
	lp->x[4] = 2;

	/* Finally, before calling the solver we can set a number of parameters
	   that control the algorithm.  A complete list of them can be
	   found in loqo.h.  There are three ways to change these parameters	
	   from their default values:
	   (1) Hardcode the changes in this file: 
	   One of the parameters is called verbose.  It is an integer 
	   variable.  The higher the value, the more output produced during the
	   solution process.  The default value is zero (no output).  We can
	   set it higher by uncommenting the following line:
		lp->verbose=2;
	   (2) Set the options using an environment variable called 
	   loqo_options:
	   If you have an environment variable called loqo_options which is 
	   equal to a string containing a list of options and their new values,
	   they will be read if you uncomment the following lines:
		opnstr = getenv("loqo_options");
		if (opnstr != NULL) {
		        if (!rd_opns(opnstr)) {
		                printf("loqo_options not set correctly\n");
		                return 0;
		        }
		} 
		set_opns(lp);
	   (3) Set the options in a specification file:
	   A file called LOQO.SPC has been provided with this distribution to
	   serve as a template.  Uncomment the following lines to read and use 
	   the file:
		bad_opns = rd_specs( "LOQO.SPC" );
		if (bad_opns) {
		        printf("Error in LOQO.SPC, line %d \n", bad_opns);
		        return 0;
		}
		set_opns( lp );
	*/

	if (lp->verbose>1) {
	    printf("%s\n%s\n%s\n%s\n%s\n%s\n",
		"\t+-------------------------------------------------+",
		"\t|                                                 |",
		"\t|       Hock and Schittkowski                     |",
		"\t|       Test Problem Number 46                    |",
		"\t|                                                 |",
		"\t+-------------------------------------------------+");
	    fflush(stdout);
	}

	status = solvelp( lp );

	inv_clo();  /* frees memory associated with matrix factorization */

	if (lp->verbose>1) {
	    printf("x: \n");
	    for (j=0; j<lp->n; j++) { 
		printf("%2d %12.6f \n", j+1, lp->x[j]);
	    }
	    printf("total time in seconds = %lf \n", cputimer() );
	}

	closelp(lp);	/* frees memory used to store the problem */

	return(0);
}

/* Function objval(x) computes and returns obj(x).
   */
double objval( double *x )
{
	return pow(x[0]-x[1],2) + pow(x[2]-1,2) 
	     + pow(x[3]-1,4) + pow(x[4]-1,6);
}

/* Function objgrad(c,x) fills in c with grad(obj(x)).
   */
void objgrad( double *c, double *x )
{
	c[0] =  2*(x[0]-x[1]);
	c[1] = -2*(x[0]-x[1]);
	c[2] =  2*(x[2]-1);
	c[3] =  4*pow(x[3]-1,3);
	c[4] =  6*pow(x[4]-1,5);
}

/* Function hessian(Q,x,y) fills in Q with 

	Hessian(obj(x)) - sum_i y_i * Hessian(a_i(x))
   */
void hessian( double *Q, double *x, double *y )
{
	int k;

	/* Initialize Q[] to zero 
	   */
	for (k=0; k<13; k++) { Q[k] = 0; }

	/* Now feed in the nonzeros associated with the Hessian of the
	   objective function.  
	   Recall the sparsity pattern for Q:
		       (var  1  2  3  4  5)
		        col  0  1  2  3  4
		var  row
		 1    0      *  *     *   
		 2    1      *  *           
		 3    2            *  *   
		 4    3      *     *  *  *
		 5    4               *  *
	    */
	Q[ 0] +=  2;
	Q[ 1] += -2;
	Q[ 3] += -2;
	Q[ 4] +=  2;
	Q[ 5] +=  2;
	Q[ 9] +=  4*3*pow(x[3]-1,2);
	Q[12] +=  6*5*pow(x[4]-1,4);

	/* Now add in y[0] times the Hessian of the first constraint 
	   */
	Q[ 0] -=  y[0]*( 2*x[3] );
	Q[ 2] -=  y[0]*( 2*x[0] );
	Q[ 7] -=  y[0]*( 2*x[0] );
	Q[ 9] -=  y[0]*( -sin(x[3]-x[4]) );
	Q[10] -=  y[0]*(  sin(x[3]-x[4]) );
	Q[11] -=  y[0]*(  sin(x[3]-x[4]) );
	Q[12] -=  y[0]*( -sin(x[3]-x[4]) );

	/* Now add in y[1] times the Hessian of the second constraint 
	   */
	Q[ 5] -=  y[1]*( 4*3*pow(x[2],2)*pow(x[3],2) );
	Q[ 6] -=  y[1]*( 4*pow(x[2],3)*2*x[3] );
	Q[ 8] -=  y[1]*( 4*pow(x[2],3)*2*x[3] );
	Q[ 9] -=  y[1]*( pow(x[2],4)*2 );
}

/* Function conval(h, x) fills in vector h with a_i(x), i=1..m.
   */
void conval( double *h, double *x )
{
	h[0] = pow(x[0],2)*x[3] + sin(x[3]-x[4]);
	h[1] = x[1] + pow(x[2],4)*pow(x[3],2);
}

/* Function congrad(A, At, x) fills in A and At with grad(a_i(x)), i=1..m.
   */
void congrad( double *A, double *At, double *x )
{
	/*  To fill in the values of A[], recall the sparsity pattern for A:
	  (var  1  2  3  4  5)
	   col  0  1  2  3  4
	row 
	 0      *        *  *
	 1         *  *  *
	 */
	A[0] = 2*x[0]*x[3];
	A[1] = 1;
	A[2] = 4*pow(x[2],3)*pow(x[3],2);
	A[3] = pow(x[0],2) + cos(x[3]-x[4]);
	A[4] = 2*pow(x[2],4)*x[3];
	A[5] = -cos(x[3]-x[4]);

	/* At this point, both A and its transpose, At, exist.  At[] must be 
	   updated too.  To do it, read the elements of A[] rowwise.
	   For big problems, libloqo.a has a function atnum(...) that
	   will recompute the entries of At[].  It is probably more efficient
	   to do it directly as shown below.
	   */
	At[0] = A[0];
	At[1] = A[3];
	At[2] = A[5];
	At[3] = A[1];
	At[4] = A[2];
	At[5] = A[4];
}

#ifdef __cplusplus
}
#endif