regsmoimproved.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.
 */

/*
 *    RegSMOImproved.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.TechnicalInformation;
import weka.core.TechnicalInformationHandler;
import weka.core.Utils;
import weka.core.TechnicalInformation.Field;
import weka.core.TechnicalInformation.Type;

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

/**
 <!-- globalinfo-start -->
 * Learn SVM for regression using SMO with Shevade, Keerthi, et al. adaption of the stopping criterion.<br/>
 * <br/>
 * For more information see:<br/>
 * <br/>
 * S.K. Shevade, S.S. Keerthi, C. Bhattacharyya, K.R.K. Murthy: Improvements to the SMO Algorithm for SVM Regression. In: IEEE Transactions on Neural Networks, 1999.<br/>
 * <br/>
 * S.K. Shevade, S.S. Keerthi, C. Bhattacharyya, K.R.K. Murthy (1999). Improvements to the SMO Algorithm for SVM Regression. Control Division, Dept. of Mechanical Engineering.
 * <p/>
 <!-- globalinfo-end -->
 *
 <!-- technical-bibtex-start -->
 * BibTeX:
 * <pre>
 * &#64;inproceedings{Shevade1999,
 *    author = {S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy},
 *    booktitle = {IEEE Transactions on Neural Networks},
 *    title = {Improvements to the SMO Algorithm for SVM Regression},
 *    year = {1999},
 *    PS = {http://guppy.mpe.nus.edu.sg/\~mpessk/svm/ieee_smo_reg.ps.gz}
 * }
 * 
 * &#64;techreport{Shevade1999,
 *    address = {Control Division, Dept. of Mechanical Engineering},
 *    author = {S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy},
 *    institution = {National University of Singapore},
 *    number = {CD-99-16},
 *    title = {Improvements to the SMO Algorithm for SVM Regression},
 *    year = {1999},
 *    PS = {http://guppy.mpe.nus.edu.sg/\~mpessk/svm/smoreg_mod.ps.gz}
 * }
 * </pre>
 * <p/>
 <!-- technical-bibtex-end -->
 *
 <!-- options-start -->
 * Valid options are: <p/>
 * 
 * <pre> -T &lt;double&gt;
 *  The tolerance parameter for checking the stopping criterion.
 *  (default 0.001)</pre>
 * 
 * <pre> -V
 *  Use variant 1 of the algorithm when true, otherwise use variant 2.
 *  (default true)</pre>
 * 
 * <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 RegSMOImproved
  extends RegSMO
  implements TechnicalInformationHandler {
  
  /** for serialization */
  private static final long serialVersionUID = 471692841446029784L;
  
  public final static int I0 = 3;
  public final static int I0a = 1;
  public final static int I0b = 2;
  public final static int I1 = 4;
  public final static int I2 = 8;
  public final static int I3 = 16;
  
  /** The different sets used by the algorithm. */
  protected SMOset m_I0;
  
  /** Index set {i: 0 < m_alpha[i] < C || 0 < m_alphaStar[i] < C}} */
  protected int [] m_iSet;
  
  /** b.up and b.low boundaries used to determine stopping criterion */
  protected double m_bUp, m_bLow;
  
  /** index of the instance that gave us b.up and b.low */
  protected int m_iUp, m_iLow;
  
  /** tolerance parameter used for checking stopping criterion b.up < b.low + 2 tol */
  double m_fTolerance = 0.001;
  
  /** set true to use variant 1 of the paper, otherwise use variant 2 */
  boolean m_bUseVariant1 = true;
  
  /**
   * Returns a string describing the object
   * 
   * @return 		a description suitable for
   * 			displaying in the explorer/experimenter gui
   */
  public String globalInfo() {
    return 
        "Learn SVM for regression using SMO with Shevade, Keerthi, et al. " 
      + "adaption of the stopping criterion.\n\n"
      + "For more information see:\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;
    TechnicalInformation	additional;
    
    result = new TechnicalInformation(Type.INPROCEEDINGS);
    result.setValue(Field.AUTHOR, "S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy");
    result.setValue(Field.TITLE, "Improvements to the SMO Algorithm for SVM Regression");
    result.setValue(Field.BOOKTITLE, "IEEE Transactions on Neural Networks");
    result.setValue(Field.YEAR, "1999");
    result.setValue(Field.PS, "http://guppy.mpe.nus.edu.sg/~mpessk/svm/ieee_smo_reg.ps.gz");
    
    additional = result.add(Type.TECHREPORT);
    additional.setValue(Field.AUTHOR, "S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy");
    additional.setValue(Field.TITLE, "Improvements to the SMO Algorithm for SVM Regression");
    additional.setValue(Field.INSTITUTION, "National University of Singapore");
    additional.setValue(Field.ADDRESS, "Control Division, Dept. of Mechanical Engineering");
    additional.setValue(Field.NUMBER, "CD-99-16");
    additional.setValue(Field.YEAR, "1999");
    additional.setValue(Field.PS, "http://guppy.mpe.nus.edu.sg/~mpessk/svm/smoreg_mod.ps.gz");
    
    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 tolerance parameter for checking the stopping criterion.\n" 
	+ "\t(default 0.001)", 
	"T", 1, "-T <double>"));
    
    result.addElement(new Option(
	"\tUse variant 1 of the algorithm when true, otherwise use variant 2.\n" 
	+ "\t(default true)", 
	"V", 0, "-V"));
    
    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> -T &lt;double&gt;
   *  The tolerance parameter for checking the stopping criterion.
   *  (default 0.001)</pre>
   * 
   * <pre> -V
   *  Use variant 1 of the algorithm when true, otherwise use variant 2.
   *  (default true)</pre>
   * 
   * <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('T', options);
    if (tmpStr.length() != 0) {
      setTolerance(Double.parseDouble(tmpStr));
    } else {
      setTolerance(0.001);
    }
    
    setUseVariant1(Utils.getFlag('V', options));
    
    super.setOptions(options);
  }
  
  /**
   * Gets the current settings of the object.
   *
   * @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("-T");
    result.add("" + getTolerance());
    
    if (m_bUseVariant1)
      result.add("-V");

    return (String[]) result.toArray(new String[result.size()]);	  
  }
  
  /**
   * Returns the tip text for this property
   * 
   * @return 		a description suitable for
   * 			displaying in the explorer/experimenter gui
   */
  public String toleranceTipText() {
    return "tolerance parameter used for checking stopping criterion b.up < b.low + 2 tol";
  }
  
  /**
   * returns the current tolerance
   * 
   * @return	the tolerance
   */
  public double getTolerance() {
    return m_fTolerance;
  }
  
  /**
   * sets the tolerance
   * 
   * @param d	the new tolerance
   */
  public void setTolerance(double d) {
    m_fTolerance = d;
  }
  
  /**
   * Returns the tip text for this property
   * 
   * @return 		a description suitable for
   * 			displaying in the explorer/experimenter gui
   */
  public String useVariant1TipText() {
    return "set true to use variant 1 of the paper, otherwise use variant 2.";
  }
  
  /**
   * Whether variant 1 is used
   * 
   * @return		true if variant 1 is used
   */
  public boolean isUseVariant1() {
    return m_bUseVariant1;
  }
  
  /**
   * Sets whether to use variant 1
   * 
   * @param b		if true then variant 1 is used
   */
  public void setUseVariant1(boolean b) {
    m_bUseVariant1 = b;
  }
  
  /** 
   * takeStep method from Shevade et al.s paper.
   * parameters correspond to pseudocode from paper.
   * 
   * @param i1
   * @param i2
   * @param alpha2
   * @param alpha2Star
   * @param phi2
   * @return
   * @throws Exception
   */
  protected int takeStep(int i1, int i2, double alpha2, double alpha2Star, double phi2) throws Exception {
    //procedure takeStep(i1, i2)
    //
    //  if (i1 == i2) 
    //    return 0 
    if (i1 == i2) {
      return 0;
    }
    double C1 = m_C * m_data.instance(i1).weight();
    double C2 = m_C * m_data.instance(i2).weight();
    //  alpha1, alpha1' = Lagrange multipliers for i1 
    double alpha1 = m_alpha[i1];
    double alpha1Star = m_alphaStar[i1];
//  double y1 = m_target[i1];
    // TODO: verify we do not need to recompute m_error[i1] here
    // TODO: since m_error is only updated for indices in m_I0
    double phi1 = m_error[i1];
//  if ((m_iSet[i1] & I0)==0) {
//  phi1 = -SVMOutput(i1) - m_b + m_target[i1];
//  m_error[i1] = phi1;
//  }
    //  k11 = kernel(point[i1], point[i1]) 
    //  k12 = kernel(point[i1], point[i2]) 
    //  k22 = kernel(point[i2], point[i2]) 
    //  eta = -2*k12+k11+k22 
    //  gamma = alpha1-alpha1'+alpha2-alpha2'
    //
    double k11 = m_kernel.eval(i1, i1, m_data.instance(i1));
    double k12 = m_kernel.eval(i1, i2, m_data.instance(i1));
    double k22 = m_kernel.eval(i2, i2, m_data.instance(i2));
    double eta = -2 * k12 + k11 + k22;
    double gamma = alpha1 - alpha1Star + alpha2 - alpha2Star;
//  if (eta < 0) {
    // this may happen due to numeric instability
    // due to Mercer's condition, this should not happen, hence we give up
//  return 0;
//  }
    //  % We assume that eta > 0. Otherwise one has to repeat the complete 
    //  % reasoning similarly (i.e. compute objective functions at L and H 
    //  % and decide which one is largest
    //
    //  case1 = case2 = case3 = case4 = finished = 0 
    //  alpha1old = alpha1, 
    //  alpha1old' = alpha1' 
    //  alpha2old = alpha2, 
    //  alpha2old' = alpha2' 
    //  deltaphi = F1 - F2 
    //
    
    //  while !finished
    //    % This loop is passed at most three times 
    //    % Case variables needed to avoid attempting small changes twice 
    //    if (case1 == 0) &&
    //       (alpha1 > 0 || (alpha1' == 0 && deltaphi > 0)) && 
    //       (alpha2 > 0 || (alpha2' == 0 && deltaphi < 0))
    //        compute L, H (w.r.t. 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 
    //    elseif (case2 == 0) &&
    //           (alpha1 > 0 || (alpha1' == 0 && deltaphi > 2*epsilon)) && 
    //           (alpha2' > 0 || (alpha2 == 0 && deltaphi > 2*epsilon))
    //
    //        compute L, H (w.r.t. 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 
    //    elseif (case3 == 0) &&
    //           (alpha1' > 0 || (alpha1 == 0 && deltaphi < -2*epsilon)) && 
    //           (alpha2 > 0 || (alpha2' == 0 && deltaphi < -2*epsilon))
    //         compute L, H (w.r.t. 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 
    //    elseif (case4 == 0) &&
    //           (alpha1' > 0) || (alpha1 == 0 && deltaphi < 0)) && 
    //           (alpha2' > 0) || (alpha2 == 0 && deltaphi > 0))
    //         compute L, H (w.r.t. 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 
    //         case4 = 1 
    //    else
    //      finished = 1 
    //    endif 
    //    update deltaphi 
    //  endwhile 
    
    double alpha1old = alpha1;
    double alpha1Starold = alpha1Star;
    double alpha2old = alpha2;
    double alpha2Starold = alpha2Star;
    double deltaPhi = phi1 - phi2;
    
    if (findOptimalPointOnLine(i1, alpha1, alpha1Star, C1, i2, alpha2, alpha2Star, C2, gamma, eta, deltaPhi)) {
      
      alpha1 = m_alpha[i1];
      alpha1Star = m_alphaStar[i1];
      alpha2 = m_alpha[i2];
      alpha2Star = m_alphaStar[i2];
      
      //  if changes in alpha('), alpha2(') are larger than some eps
      //    Update f-cache[i] for i in I.0 using new Lagrange multipliers 
      //    Store the changes in alpha, alpha' array 
      //    Update I.0, I.1, I.2, I.3 
      //    Compute (i.low, b.low) and (i.up, b.up) by applying the conditions mentioned above, using only i1, i2 and indices in I.0 
      //    return 1 
      //  else
      //    return 0
      //endif endprocedure
      
      //		Update error cache using new Lagrange multipliers 
      double dAlpha1 = alpha1 - alpha1old - (alpha1Star - alpha1Starold);
      double dAlpha2 = alpha2 - alpha2old - (alpha2Star - alpha2Starold);
      for (int j = m_I0.getNext(-1); j != -1; j = m_I0.getNext(j)) {
	if ((j != i1) && (j != i2)) {
	  m_error[j] -= dAlpha1 * m_kernel.eval(i1, j, m_data.instance(i1)) 
	  + dAlpha2 * m_kernel.eval(i2, j, m_data.instance(i2));
	}
      }
      m_error[i1] -= dAlpha1 * k11 + dAlpha2 * k12;
      m_error[i2] -= dAlpha1 * k12 + dAlpha2 * k22;
      
      updateIndexSetFor(i1, C1);
      updateIndexSetFor(i2, C2);
      
      //    Compute (i.low, b.low) and (i.up, b.up) by applying the conditions mentioned above, using only i1, i2 and indices in I.0