penaltyconstrainedminimizer.java

Standord Classifier实现了一个基于Java的最大熵分类器。用于模式识别

package edu.stanford.nlp.optimization;

/**
 * A wrapper class which builds a constrained minimizer out of an
 * unconstrained one by finding the unconstrained minimum of a
 * sequence of penalized surfaces.  It assumes that the objective it
 * is given is a <code>DiffFunction</code>, even if the base
 * unconstrained minimzer just ignores the derivative information.
 *
 * The desired unconstrained minimizer is passed in on construction:
 *
 * <code>Minimizer m = new SomeUnconstrainedMinimizer();</code>
 * <code>ConstrainedMinimizer cm = new PenaltyConstrainedMinimizer(m);</code>
 *
 * @author <a href="mailto:[e-mail]">Dan Klein</a>
 * @version 1.0
 * @since 1.0
 * @see ConstrainedMinimizer
 */
public class PenaltyConstrainedMinimizer implements ConstrainedMinimizer {
  private Minimizer minimizer;

  private boolean silent = true;

  class PenalizedFunction implements DiffFunction {
    DiffFunction function;

    DiffFunction[] eqConstraints;
    double[] eqLagrange;
    double[] eqPenalties;

    DiffFunction[] ineqConstraints;
    double[] ineqLagrange;
    double[] ineqPenalties;

    public int domainDimension() {
      return function.domainDimension();
    }

    public double valueAt(double[] x) {
      double val = 0.0;
      val += function.valueAt(x);
      // equality constraints
      for (int i=0; i<eqConstraints.length; i++) {
	double c = eqConstraints[i].valueAt(x);
	val += eqLagrange[i]*c;
	val += 0.5*eqPenalties[i]*c*c;
      }
      // inequality constraints
      for (int i=0; i<ineqConstraints.length; i++) {
	double c = ineqConstraints[i].valueAt(x);
	//val += Math.exp(-1.0*ineqPenalties[i]*c);
	double tmp = fmax(c,2.0*c+ineqLagrange[i]/ineqPenalties[i]);
	val += ineqLagrange[i]*tmp;
	val += 0.5*ineqPenalties[i]*tmp*tmp;
	//val += (tmp2*tmp2 - ineqLagrange[i]*ineqLagrange[i])/2.0*ineqPenalties[i];
      }
      return val;
    }

    public double[] derivativeAt(double[] x) {
      double[] deriv = copyArray(function.derivativeAt(x));
      // equality constraints
      for (int i=0; i<eqConstraints.length; i++) {
	double[] cDeriv = eqConstraints[i].derivativeAt(x);
	double cVal = eqConstraints[i].valueAt(x);
	for (int d=0; d<domainDimension(); d++) {
	  deriv[d] += eqLagrange[i]*cDeriv[d];
	  deriv[d] += eqPenalties[i]*cVal*cDeriv[d];
	}
      }
      // inequality constraints
      for (int i=0; i<ineqConstraints.length; i++) {
	double[] cDeriv = ineqConstraints[i].derivativeAt(x);
	double cVal = ineqConstraints[i].valueAt(x);
	for (int d=0; d<domainDimension(); d++) {
	  //deriv[d] += ineqLagrange[i]*cDeriv[d];
	  //deriv[d] += -1.0*ineqPenalties[i]*Math.exp(-1.0*ineqPenalties[i]*cVal)*cDeriv[d];
	  //deriv[d] += (cVal > 0 ? 0 : ineqPenalties[i]*cVal*cDeriv[d]);
	  double tmp = cVal;
	  double tmpD = cDeriv[d];
	  double crit = cVal + ineqLagrange[i]/ineqPenalties[i];
	  if (crit > 0) {
	    tmp += crit;
	    tmpD += cDeriv[d];
	  }
	  deriv[d] += ineqLagrange[i]*tmpD;
	  deriv[d] += ineqPenalties[i]*tmp*tmpD;
	}
      }
      return deriv;
    }

    PenalizedFunction(Function function, Function[] eqConstraints, double[] eqLagrange, double[] eqPenalties, Function[] ineqConstraints, double[] ineqLagrange, double[] ineqPenalties) {
      this.function = (DiffFunction)function;
      this.eqConstraints = (DiffFunction[])eqConstraints;
      this.eqLagrange = eqLagrange;
      this.eqPenalties = eqPenalties;
      this.ineqConstraints = (DiffFunction[])ineqConstraints;
      this.ineqLagrange = ineqLagrange;
      this.ineqPenalties = ineqPenalties;
    }
  }

  double[] copyArray(double[] a) {
    double[] result = new double[a.length];
    for(int i=0;i<a.length;i++)
      result[i]=a[i];
    return result;
  }

  private double arrayMax(double[] x) {
    double max = Double.NEGATIVE_INFINITY;
    for (int i=0; i<x.length; i++) {
      if (max < x[i])
	max = x[i];
    }
    return max;
  }

  private double fabs(double x) {
    if (x<0)
      return -1.0*x;
    return x;
  }

  private double fmax(double x, double y) {
    if (x<y)
      return y;
    return x;
  }

  private String arrayToString(double[] x) {
    StringBuffer sb = new StringBuffer("(");
    for (int j=0; j<x.length; j++) {
      sb.append(x[j]);
      if (j != x.length-1)
	sb.append(", ");
    }
    sb.append(")");
    return sb.toString();
  }

  public double[] minimize(Function function, double functionTolerance, Function[] eqConstraints, double eqConstraintTolerance, Function[] ineqConstraints, double ineqConstraintTolerance, double[] initial) {

    int dimension = function.domainDimension();
    int numEqConstraints = eqConstraints.length;
    int numIneqConstraints = ineqConstraints.length;

    // check that the function is a DiffFunction
    if (! (function instanceof DiffFunction))
      throw new UnsupportedOperationException();

    // check the constraints
    for (int i=0; i<numEqConstraints; i++) {
      if (! (eqConstraints[i] instanceof DiffFunction))
	throw new UnsupportedOperationException();
    }
    for (int i=0; i<numIneqConstraints; i++) {
      if (! (ineqConstraints[i] instanceof DiffFunction))
	throw new UnsupportedOperationException();
    }

    // use a penalty method to solve for the lagrange multipliers and the x

    double[] eqLagrange = new double[numEqConstraints];
    double[] eqPenalties = new double[numEqConstraints];
    for (int i=0; i<numEqConstraints; i++) {
      eqLagrange[i] = 0;
      eqPenalties[i] = 1.0;
    }

    double[] ineqLagrange = new double[numIneqConstraints];
    double[] ineqPenalties = new double[numIneqConstraints];
    for (int i=0; i<numIneqConstraints; i++) {
      ineqLagrange[i] = 0;
      ineqPenalties[i] = 1.0;
    }

    double worstEqViolation = 1.0+2.0*eqConstraintTolerance;
    double worstIneqViolation = 1.0+2.0*ineqConstraintTolerance;
    double lastWorstEqViolation = 10*worstEqViolation;
    double lastWorstIneqViolation = 10*worstIneqViolation;

    int iter = 0;

    double[] x = copyArray(initial);

    // do a series of penalized minimizations
    boolean lowTolIter = false;
    while (iter < 100 && 
	   (worstEqViolation > eqConstraintTolerance ||
	    worstIneqViolation > ineqConstraintTolerance ||
	    !lowTolIter)) {

      if (worstEqViolation <= eqConstraintTolerance &&
	  worstIneqViolation <= ineqConstraintTolerance) {
	lowTolIter = true;
      } else {
	lowTolIter = false;
      }

      // build a penalized surface
      DiffFunction penalizedFunction = new PenalizedFunction(function, eqConstraints, eqLagrange, eqPenalties, ineqConstraints, ineqLagrange, ineqPenalties);

      // minimize it
      double thisTol = (lowTolIter ? functionTolerance : Math.sqrt(functionTolerance));
      x = minimizer.minimize(penalizedFunction, thisTol, x);
      penalizedFunction.derivativeAt(x);

      // update lagrange multipliers and penalties

      // EQUALITY CONSTRAINTS
      lastWorstEqViolation = worstEqViolation;
      worstEqViolation = Double.NEGATIVE_INFINITY;
      double[] eqViolations = new double[numEqConstraints];
      double[] eqValues = new double[numEqConstraints];
      for (int i=0; i<numEqConstraints; i++) {
	eqValues[i] = eqConstraints[i].valueAt(x);
	eqViolations[i] = fabs(eqValues[i]);
      }
      worstEqViolation = arrayMax(eqViolations);
      // update penalties and lagrange multipliers
      for (int i=0; i<numEqConstraints; i++) {
	// lagrange multipliers should take over the forces previously
	// exerted by the penalty terms
	eqLagrange[i] += eqPenalties[i]*eqValues[i];
	// penalties increase more or less based on violation severity
	if (eqViolations[i] >= worstEqViolation/2.0 &&
	    worstEqViolation >= eqConstraintTolerance)
	  eqPenalties[i] *= 1.0;//*= 1.2;
	else
	  eqPenalties[i] *= 1.0;//*= 2.0;
	// penalties also increase if constraint satisfaction is too slow
	if (worstEqViolation / (lastWorstEqViolation + 1e-100) > 0.25 &&
	    worstEqViolation >= eqConstraintTolerance)
	  eqPenalties[i] *= 5.0;
      }

      if (! silent)
	System.err.println("!e "+arrayMax(eqPenalties)+"/"+worstEqViolation+"!");

      // INEQUALTIY CONSTRAINTS
      lastWorstIneqViolation = worstIneqViolation;
      worstIneqViolation = Double.NEGATIVE_INFINITY;
      double[] ineqViolations = new double[numIneqConstraints];
      double[] ineqValues = new double[numIneqConstraints];
      for (int i=0; i<numIneqConstraints; i++) {
	// ineq violations can be < 0 constraint values OR
	//   >= 0 constraints with non-zero lagrange multipliers / penalties
	double cVal = ineqConstraints[i].valueAt(x);
	ineqValues[i] = cVal;
	ineqViolations[i] = ((cVal < 0 ? -1.0*cVal : 0)+fabs(ineqLagrange[i]*cVal));
      }
      worstIneqViolation = fmax(0.0,arrayMax(ineqViolations));
      // update penalties and lagrange multipliers
      for (int i=0; i<numIneqConstraints; i++) {
	// lagrange multipliers should take over the forces previously
	// exerted by the penalty terms ... sort of
	double crit =  ineqValues[i] + ineqLagrange[i]/ineqPenalties[i];
	if (crit > 0)
	  ineqLagrange[i] += ineqPenalties[i]*crit;
	ineqLagrange[i] += ineqPenalties[i]*ineqValues[i];
	if (ineqLagrange[i] > 0)
	  ineqLagrange[i] = 0;
	// penalties increase more or less based on violation severity
	if (ineqViolations[i] >= worstIneqViolation/2.0 &&
	    worstIneqViolation >= ineqConstraintTolerance)
	  ineqPenalties[i] *= 1.0;//+= 0.2;
	else
	  ineqPenalties[i] *= 1.0;//+= 1.0;
	// penalties also increase if constraint satisfaction is too slow
	if (worstIneqViolation / (lastWorstIneqViolation + 1e-100) > 0.25 &&
	    worstIneqViolation >= ineqConstraintTolerance)
	  ineqPenalties[i] *= 5.0;
      }

      if (! silent)
	System.err.println("!i "+arrayMax(ineqPenalties)+"/"+worstIneqViolation+"!");

    }

    return copyArray(x);
  }


  public double[] minimize(Function function, double functionTolerance, double[] initial) {
    return minimizer.minimize(function, functionTolerance, initial);
  }

  public PenaltyConstrainedMinimizer(Minimizer minimizer) {
    this.minimizer = minimizer;
  }
}