regsmo.java

Weka

/*
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 *    RegSMO.java
 *    Copyright (C) 2006 University of Waikato, Hamilton, New Zealand
 *
 */

package weka.classifiers.functions.supportVector;

import weka.core.Instances;
import weka.core.Option;
import weka.core.Utils;
import weka.core.TechnicalInformation;
import weka.core.TechnicalInformationHandler;
import weka.core.TechnicalInformation.Field;
import weka.core.TechnicalInformation.Type;

import java.util.Enumeration;
import java.util.Vector;

/**
 <!-- globalinfo-start -->
 * Implementation of SMO for support vector regression as described in :<br/>
 * <br/>
 * A.J. Smola, B. Schoelkopf (1998). A tutorial on support vector regression.
 * <p/>
 <!-- globalinfo-end -->
 *
 <!-- technical-bibtex-start -->
 * BibTeX:
 * <pre>
 * &#64;misc{Smola1998,
 *    author = {A.J. Smola and B. Schoelkopf},
 *    note = {NeuroCOLT2 Technical Report NC2-TR-1998-030},
 *    title = {A tutorial on support vector regression},
 *    year = {1998}
 * }
 * </pre>
 * <p/>
 <!-- technical-bibtex-end -->
 *
 <!-- options-start -->
 * Valid options are: <p/>
 * 
 * <pre> -P &lt;double&gt;
 *  The epsilon for round-off error.
 *  (default 1.0e-12)</pre>
 * 
 * <pre> -L &lt;double&gt;
 *  The epsilon parameter in epsilon-insensitive loss function.
 *  (default 1.0e-3)</pre>
 * 
 * <pre> -W &lt;double&gt;
 *  The random number seed.
 *  (default 1)</pre>
 * 
 <!-- options-end -->
 *
 * @author  Remco Bouckaert ([e-mail],[e-mail])
 * @version $Revision: 1.3 $
 */
public class RegSMO 
  extends RegOptimizer
  implements TechnicalInformationHandler {
  
  /** for serialization */
  private static final long serialVersionUID = -7504070793279598638L;

  /** tolerance parameter, smaller changes on alpha in inner loop will be ignored **/
  protected double m_eps = 1.0e-12;
  
  /** Precision constant for updating sets */
  protected final static double m_Del = 1e-10; //1000 * Double.MIN_VALUE;
  
  /** error cache containing m_error[i] = SVMOutput(i) - m_target[i] - m_b <br/>
   * note, we don't need m_b in the cache, since if we do, we need to maintain 
   * it when m_b is updated */
  double[] m_error;
  
  /** alpha value for first candidate **/  
  protected double m_alpha1;
  
  /** alpha* value for first candidate **/  
  protected double m_alpha1Star;

  /** alpha value for second candidate **/  
  protected double m_alpha2;
  
  /** alpha* value for second candidate **/
  protected double m_alpha2Star;
  
  /**
   * default constructor
   */
  public RegSMO() {
    super();
  }
  
  /**
   * Returns a string describing classifier
   * 
   * @return 		a description suitable for
   * 			displaying in the explorer/experimenter gui
   */
  public String globalInfo() {
    return
        "Implementation of SMO for support vector regression as described "
      + "in :\n\n"
      + getTechnicalInformation().toString();
  }

  /**
   * Returns an instance of a TechnicalInformation object, containing 
   * detailed information about the technical background of this class,
   * e.g., paper reference or book this class is based on.
   * 
   * @return 		the technical information about this class
   */
  public TechnicalInformation getTechnicalInformation() {
    TechnicalInformation 	result;
    
    result = new TechnicalInformation(Type.MISC);
    result.setValue(Field.AUTHOR, "A.J. Smola and B. Schoelkopf");
    result.setValue(Field.TITLE, "A tutorial on support vector regression");
    result.setValue(Field.NOTE, "NeuroCOLT2 Technical Report NC2-TR-1998-030");
    result.setValue(Field.YEAR, "1998");
    
    return result;
  }
  
  /**
   * Returns an enumeration describing the available options
   * 
   * @return an enumeration of all the available options
   */
  public Enumeration listOptions() {
    Vector result = new Vector();
    
    result.addElement(new Option(
	"\tThe epsilon for round-off error.\n" 
	+ "\t(default 1.0e-12)", 
	"P", 1, "-P <double>"));
    
    Enumeration enm = super.listOptions();
    while (enm.hasMoreElements()) {
      result.addElement(enm.nextElement());
    }
    
    return result.elements();
  }
  
  /**
   * Parses a given list of options. <p/>
   *
   <!-- options-start -->
   * Valid options are: <p/>
   * 
   * <pre> -P &lt;double&gt;
   *  The epsilon for round-off error.
   *  (default 1.0e-12)</pre>
   * 
   * <pre> -L &lt;double&gt;
   *  The epsilon parameter in epsilon-insensitive loss function.
   *  (default 1.0e-3)</pre>
   * 
   * <pre> -W &lt;double&gt;
   *  The random number seed.
   *  (default 1)</pre>
   * 
   <!-- options-end -->
   *
   * @param options the list of options as an array of strings
   * @throws Exception if an option is not supported 
   */
  public void setOptions(String[] options) throws Exception {
    String	tmpStr;
    
    tmpStr = Utils.getOption('P', options);
    if (tmpStr.length() != 0) {
      setEpsilon(Double.parseDouble(tmpStr));
    } else {
      setEpsilon(1.0e-12);
    }
    
    super.setOptions(options);
  }
  
  /**
   * Gets the current settings of the classifier.
   *
   * @return an array of strings suitable for passing to setOptions
   */
  public String[] getOptions() {
    int       	i;
    Vector    	result;
    String[]  	options;

    result = new Vector();

    options = super.getOptions();
    for (i = 0; i < options.length; i++)
      result.add(options[i]);
    
    result.add("-P");
    result.add("" + getEpsilon());

    return (String[]) result.toArray(new String[result.size()]);	  
  }
  
  /**
   * Returns the tip text for this property
   * 
   * @return 		tip text for this property suitable for
   * 			displaying in the explorer/experimenter gui
   */
  public String epsilonTipText() {
    return "The epsilon for round-off error (shouldn't be changed).";
  }
  
  /**
   * Get the value of epsilon.
   * 
   * @return 		Value of epsilon.
   */
  public double getEpsilon() {
    return m_eps;
  }
  
  /**
   * Set the value of epsilon.
   * 
   * @param v  		Value to assign to epsilon.
   */
  public void setEpsilon(double v) {
    m_eps = v;
  }
  
  /** initialize various variables before starting the actual optimizer 
   * 
   * @param data 	data set used for learning
   * @throws Exception	if something goes wrong
   */
  protected void init(Instances data) throws Exception {
    super.init(data);
    
    //init error cache
    m_error = new double[m_nInstances];
    for (int i = 0; i < m_nInstances; i++) {
      m_error[i] = -m_target[i];
    }
  }
  
  /** 
   * wrap up various variables to save memeory and do some housekeeping after optimization
   * has finished.
   *
   * @throws Exception 	if something goes wrong
   */
  protected void wrapUp() throws Exception {
    m_error = null;
    super.wrapUp();
  }

  /** 
   * Finds optimal point on line constrained by first (i1) and second (i2) 
   * candidate. Parameters correspond to pseudocode (see technicalinformation)
   * 
   * @param i1
   * @param alpha1
   * @param alpha1Star
   * @param C1
   * @param i2
   * @param alpha2
   * @param alpha2Star
   * @param C2
   * @param gamma
   * @param eta
   * @param deltaPhi
   * @return
   */
  protected boolean findOptimalPointOnLine(int i1, double alpha1, double alpha1Star, double C1, 
      int i2, double alpha2, double alpha2Star, double C2, 
      double gamma, double eta, double deltaPhi) {
    if (eta <= 0) {
      // this may happen due to numeric instability
      // due to Mercer's condition, this should not happen, hence we give up
      return false;
    }
    
    boolean case1 = false;
    boolean case2 = false;
    boolean case3 = false;
    boolean case4 = false;
    boolean finished = false;
    
    //		while !finished 
    //		% this loop is passed at most three times 
    //		% case variables needed to avoid attempting small changes twice 
    while (!finished) {
      //			if (case1 == 0) && 
      //				(alpha1 > 0 || (alpha1* == 0 && deltaPhi > 0)) && 
      //				(alpha2 > 0 || (alpha2* == 0 && deltaPhi < 0)) 
      //				compute L, H (wrt. alpha1, alpha2) 
      //				if L < H 
      //					a2 = alpha2 ? - deltaPhi/eta 
      //					a2 = min(a2, H) 
      //					a2 = max(L, a2) 
      //					a1 = alpha1 ? - (a2 ? alpha2) 
      //					update alpha1, alpha2 if change is larger than some eps 
      //				else 
      //					finished = 1 
      //				endif 
      //				case1 = 1; 
      
      if ((case1 == false) && 
	  (alpha1 > 0 || (alpha1Star == 0 && deltaPhi > 0)) && 
	  (alpha2 > 0 || (alpha2Star == 0 && deltaPhi < 0))) {
	// compute L, H (wrt. alpha1, alpha2) 
	double L = Math.max(0, gamma - C1);
	double H = Math.min(C2, gamma);
	if (L < H) {
	  double a2 = alpha2 - deltaPhi / eta;
	  a2 = Math.min(a2, H);
	  a2 = Math.max(L, a2);
	  // To prevent precision problems
	  if (a2 > C2 - m_Del * C2) {
	    a2 = C2;
	  } else if (a2 <= m_Del * C2) {
	    a2 = 0;
	  }
	  double a1 = alpha1 - (a2 - alpha2);
	  if (a1 > C1 - m_Del * C1) {
	    a1 = C1;
	  } else if (a1 <= m_Del * C1) {
	    a1 = 0;
	  }
	  // update alpha1, alpha2 if change is larger than some eps
	  if (Math.abs(alpha1 - a1) > m_eps) {
	    deltaPhi += eta * (a2 - alpha2);
	    alpha1 = a1;
	    alpha2 = a2;
	  }
	} else {
	  finished = true;
	}
	case1 = true;
      }
      
      //			elseif (case2 == 0) && 
      //				(alpha1 > 0 || (alpha1* == 0 && deltaPhi > 2 epsilon)) && 
      //				(alpha2* > 0 || (alpha2 == 0 && deltaPhi > 2 epsilon)) 
      //				compute L, H (wrt. alpha1, alpha2*) 
      //				if L < H 
      //					a2 = alpha2* + (deltaPhi ?- 2 epsilon)/eta 
      //					a2 = min(a2, H) 
      //					a2 = max(L, a2) 
      //					a1 = alpha1 + (a2 ? alpha2*) 
      //					update alpha1, alpha2* if change is larger than some eps 
      //				else 
      //					finished = 1 
      //				endif 
      //				case2 = 1; 
      
      else if (
	  (case2 == false)
	  && (alpha1 > 0 || (alpha1Star == 0 && deltaPhi > 2 * m_epsilon))
	  && (alpha2Star > 0 || (alpha2 == 0 && deltaPhi > 2 * m_epsilon))) {
	// compute L, H (wrt. alpha1, alpha2*) 
	double L = Math.max(0, -gamma);
	double H = Math.min(C2, -gamma + C1);
	if (L < H) {
	  double a2 = alpha2Star + (deltaPhi - 2 * m_epsilon) / eta;
	  a2 = Math.min(a2, H);
	  a2 = Math.max(L, a2);
	  // To prevent precision problems
	  if (a2 > C2 - m_Del * C2) {
	    a2 = C2;
	  } else if (a2 <= m_Del * C2) {
	    a2 = 0;
	  }
	  double a1 = alpha1 + (a2 - alpha2Star);
	  if (a1 > C1 - m_Del * C1) {
	    a1 = C1;
	  } else if (a1 <= m_Del * C1) {
	    a1 = 0;
	  }
	  // update alpha1, alpha2* if change is larger than some eps 
	  if (Math.abs(alpha1 - a1) > m_eps) {
	    deltaPhi += eta * (-a2 + alpha2Star);
	    alpha1 = a1;
	    alpha2Star = a2;
	  }
	} else {
	  finished = true;
	}
	case2 = true;
      }
      
      //			elseif (case3 == 0) && 
      //				(alpha1* > 0 || (alpha1 == 0 && deltaPhi < -2 epsilon)) && 
      //				(alpha2 > 0 || (alpha2* == 0 && deltaPhi < -2 epsilon)) 
      //				compute L, H (wrt. alpha1*, alpha2) 
      //				if L < H 
      //					a2 = alpha2 ?- (deltaPhi ?+ 2 epsilon)/eta 
      //					a2 = min(a2, H) 
      //					a2 = max(L, a2) 
      //					a1 = alpha1* + (a2 ? alpha2) 
      //					update alpha1*, alpha2 if change is larger than some eps 
      //				else 
      //					finished = 1 
      //				endif 
      //				case3 = 1; 
      
      else if (
	  (case3 == false)
	  && (alpha1Star > 0 || (alpha1 == 0 && deltaPhi < - 2 * m_epsilon))
	  && (alpha2 > 0 || (alpha2Star == 0 && deltaPhi < - 2 * m_epsilon))) {
	// compute L, H (wrt. alpha1*, alpha2)
	double L = Math.max(0, gamma);
	double H = Math.min(C2, C1 + gamma);
	if (L < H) {
	  // note Smola's psuedocode has a minus, where there should be a plus in the following line, Keerthi's is correct